From 29e6840dce094fa3c9550797ddb9c9812f092e1f Mon Sep 17 00:00:00 2001 From: ClothildeCorsi Date: Fri, 25 Jun 2021 19:04:05 +1000 Subject: [PATCH 1/9] Implementation of a Beam-Down system with 2D-crossed CPC (Compound Parabolic Concentrator) and secondary reflector (hyperboloid) with vertical axis and equivalent eccentricity in xOz and yOz planes. --- example/gen_bd.py | 102 ++++++ solsticepy/__init__.py | 1 + solsticepy/cal_layout.py | 269 ++++----------- solsticepy/cal_sun.py | 8 +- solsticepy/design_bd.py | 677 +++++++++++++++++++++++++++++++++++++ solsticepy/design_cpc.py | 408 ++++++++++++++++++++++ solsticepy/design_crs.py | 189 +++++------ solsticepy/design_dish.py | 10 +- solsticepy/gen_yaml.py | 299 ++++++++-------- solsticepy/input.py | 139 ++++---- solsticepy/master.py | 80 ++--- solsticepy/output_motab.py | 182 ++-------- solsticepy/process_raw.py | 159 ++++----- 13 files changed, 1702 insertions(+), 821 deletions(-) create mode 100644 example/gen_bd.py create mode 100644 solsticepy/design_bd.py create mode 100644 solsticepy/design_cpc.py diff --git a/example/gen_bd.py b/example/gen_bd.py new file mode 100644 index 0000000..6a787ed --- /dev/null +++ b/example/gen_bd.py @@ -0,0 +1,102 @@ +# +# This is an example script to simulate a central receiver system (CRS) +# using Solstice via solsticepy +# +import solsticepy +from solsticepy.design_bd import * +import numpy as np +import os + +#================================================== +# INPUT PARAMETERS + +# whether run a new case (True) or load pre-existing input.yaml inputs (False): +new_case=True + +# define a unique case folder for the user +# ========= + +snum = 0 +suffix = "" +while 1: + import datetime,time + dt = datetime.datetime.now() + ds = dt.strftime("%a-%H-%M") + casefolder = os.path.join(os.getcwd(),'case-%s-%s%s'%(os.path.basename(__file__),ds,suffix)) + if os.path.exists(casefolder): + snum+=1 + suffix = "-%d"%(snum,) + if snum > 200: + raise RuntimeError("Some problem with creating casefolder") + else: + # good, we have a new case dir + break + +pm=Parameters() + +# Variables +pm.H_tower=65. # tower height or vertical distance to aiming point (located at the center of xOy plan) +theta_deg=20. # acceptance half angle of the CPC in degree +rim_angle = 45. # rim angle of the heliostat field in the xOz plan in degree +foci_ratio = 0.6 # ratio of the foci distance and apex distance to the origin [1/2,1] + +# fixed parameters +# ========= + +pm.Q_in_rcv=40e6 +num_hst = 300 +pm.nd=5 +pm.nh=5 +pm.H_rcv=10. +pm.W_rcv=10. +pm.dependent_par() +pm.saveparam(casefolder) +print(pm.fb) +print(pm.H_tower) +pm.field_type = 'surround' +pm.Zhelio = 0. +pm.slope_error = 0. +pm.R1 = 30. + + +# enter the parameters for the beam-down components +receiver='beam_down' +# Polygon receiver +rec_w=1.2 +rec_l=10. +rec_z=0. +# 2D-crossed cpc with n faces +n_CPC_faces=4 +rec_grid=200 +n_Z=30 +# Secondary refector 'hyperboloid' +refl_sec = 0.95 +slope_error = 2.e-3 # radian + +# parameters recalculated (pre-optimized before optimization) +secref_vert = None # np.array([[-15,25],[-15,-25],[15,-25],[15,25]]) +cpc_h=None + +# create the environment and scene +# ========= +tablefile=casefolder+'/OELT_Solstice.motab' +weafile='./demo_TMY3_weather.motab' + +bd=BD(latitude=pm.lat, casedir=casefolder) + +bd.receiversystem(receiver=receiver, rec_abs=float(pm.alpha_rcv), rec_w=float(rec_w), rec_l=float(rec_l), rec_z=float(rec_z), rec_grid=int(rec_grid), cpc_nfaces=int(n_CPC_faces), cpc_theta_deg=float(theta_deg), cpc_h=cpc_h, cpc_nZ=float(n_Z), field_rim_angle=float(rim_angle), aim_z=float(pm.H_tower), secref_fratio=foci_ratio, refl_sec=float(refl_sec), slope_error=float(slope_error), secref_vert = secref_vert) + +bd.heliostatfield(field=pm.field_type, hst_rho=pm.rho_helio, slope=pm.slope_error, hst_w=pm.W_helio, hst_h=pm.H_helio, tower_h=pm.H_tower, tower_r=pm.R_tower, hst_z=pm.Z_helio, num_hst=num_hst, R1=pm.R1, fb=pm.fb, dsep=pm.dsep, x_max=150., y_max=150.) + + + +bd.yaml(sunshape=pm.sunshape,csr=pm.crs,half_angle_deg=pm.half_angle_deg,std_dev=pm.std_dev) + +oelt, A_land=bd.field_design_annual(dni_des=900., num_rays=int(1e5), nd=pm.nd, nh=pm.nh, weafile=weafile, method=1, Q_in_des=pm.Q_in_rcv, n_helios=None, zipfiles=False, gen_vtk=True, plot=False) + + +if (A_land==0): + tablefile=None +else: + A_helio=pm.H_helio*pm.W_helio + output_matadata_motab(table=oelt, field_type=pm.field_type, aiming='single', n_helios=num_hst, A_helio=A_helio, eff_design=bd.eff_des, H_rcv=pm.H_rcv, W_rcv=pm.W_rcv, H_tower=pm.H_tower, Q_in_rcv=pm.Q_in_rcv, A_land=A_land, savedir=tablefile) diff --git a/solsticepy/__init__.py b/solsticepy/__init__.py index 67e58a1..8b6bef5 100644 --- a/solsticepy/__init__.py +++ b/solsticepy/__init__.py @@ -8,3 +8,4 @@ from .gen_yaml import * from .process_raw import * from .master import * +from .design_cpc import * diff --git a/solsticepy/cal_layout.py b/solsticepy/cal_layout.py index 7de3c46..9b18627 100644 --- a/solsticepy/cal_layout.py +++ b/solsticepy/cal_layout.py @@ -1,13 +1,12 @@ import numpy as np import sys -import os import matplotlib -#matplotlib.use("agg") +matplotlib.use("agg") import matplotlib.pyplot as plt from .cal_sun import * from .gen_vtk import gen_vtk -def radial_stagger(latitude, num_hst, width, height, hst_z, towerheight, R1, fb, dsep=0., field='polar', num_aperture=0, gamma=0., rec_w=0., rec_z=[], savedir='.', verbose=False, plot=False, plt_aiming=None): +def radial_stagger(latitude, num_hst, width, height, hst_z, towerheight, R1, fb, dsep=0., field='polar', savedir='.', plot=False, r_field_max=-1.): '''Generate a radial-stagger heliostat field, ref. Collado and Guallar, 2012, Campo: Generation of regular heliostat field. ``Arguments`` @@ -20,13 +19,10 @@ def radial_stagger(latitude, num_hst, width, height, hst_z, towerheight, R1, fb, * R1 (float) : distance from the first row to the bottom of the tower, i.e. (0, 0, 0) * fb (float) : the field layout growing factor, in (0, 1) * dsep (float) : separation distance (m) - * field (str) : 'polar-half' or 'surround-half' or 'polar' or 'surround' field or 'multi-aperture', the 'half' option is for simulation a symmetric field - * num_aperture(int): number of apertures, for a multi-aperture configuration - * gamma (float) : the anangular range of the multi-aperture configration (deg) - * rec_z (list) : a list of the elevation heights of the apertures + * field (str) : 'polar-half' or 'surround-half' or 'polar' or 'surround' field, the 'half' option is for simulation a symmetric field * savedir (str) : directory of saving the pos_and_aiming.csv - * verbose(bool) : write results to disk or not * plot (bool) : True - plot the layout by Matplotlib + * r_field_max : Maximum field radius (m), if r_max < -1, it is not considered ``Returns`` @@ -54,8 +50,8 @@ def radial_stagger(latitude, num_hst, width, height, hst_z, towerheight, R1, fb, ''' - # heliostat diagonal distantce - DH=np.sqrt(height**2+width**2) + # heliostat diagonal distance + DH=np.sqrt(height**2+width**2) # distance between contiguous helistat center on the X and Y plane DM=DH+dsep @@ -81,7 +77,16 @@ def radial_stagger(latitude, num_hst, width, height, hst_z, towerheight, R1, fb, sys.stderr.write('DM '+repr(DM)+'\n') sys.stderr.write('dRm'+repr(delta_Rmin)+'\n') - while num 0.: + max_helio = np.inf + r_field_max *= 1.5 + else: + max_helio = num_hst*3 + r_field_max = np.inf + + r_field=0. + while num1.5*np.pi)+(azimuth<0.5*np.pi) + idx= (azimuth>-np.pi/2.)*(azimuth(1.5*np.pi+i*np.pi/40.))+(azimuth<(np.pi/2.-i*np.pi/40.)) + idx=(azimuth>-np.pi/2.+i*np.pi/40.)* (azimuth=5) * nh (int): number of columns of the lookup table (hours in a day, i.e. 24h) - * verbose (bool): write results to disk or not ``Returns`` @@ -409,7 +407,7 @@ def annual_angles(self, latitude, casefolder=None, nd=5, nh=5, verbose=False): #azimuth=case_list[1:,-2].astype(float) #zenith=case_list[1:,-1].astype(float) - if casefolder!=None and verbose: + if casefolder!='NOTSAVE': np.savetxt(casefolder+'/table_view.csv', table, fmt='%s', delimiter=',') np.savetxt(casefolder+'/annual_simulation_list.csv', case_list, fmt='%s', delimiter=',') diff --git a/solsticepy/design_bd.py b/solsticepy/design_bd.py new file mode 100644 index 0000000..da88431 --- /dev/null +++ b/solsticepy/design_bd.py @@ -0,0 +1,677 @@ +import os +import sys +import time +import numpy as np +import matplotlib.pyplot as plt +from uncertainties import ufloat +from scipy.interpolate import interp1d,interp2d +import matplotlib.cm as cm +from scipy.optimize import curve_fit + +from .process_raw import * +from .cal_layout import radial_stagger +from .cal_field import * +from .cal_sun import * +from .gen_yaml import gen_yaml, Sun +from .gen_vtk import * +from .input import Parameters +from .output_motab import output_matadata_motab, output_motab +from .master import * + + +class BD: + ''' + The Beam-Down (BD) model includes five parts: + the sun, the field, the secondary reflector, the CPC and the receiver. + ''' + + def __init__(self, latitude, casedir): + ''' + Arguements: + casedir : str, the directory of the case + ''' + self.casedir=casedir + + if not os.path.exists(casedir): + os.makedirs(casedir) + self.latitude=latitude + self.sun=SunPosition() + self.master=Master(casedir) + + def cpcradius(self, rec_w=2., rec_l=2., cpc_nfaces=4): + ''' + Ouput + rec_radius_ref: receiver radius in the direction of theta (CPC half acceptance angle), it is used for the design criteria + rec_radius: receiver radius in the opposite direction (equal to rec_radius_theta if nber of CPC faces is different to 4, because symmetry applies) + ''' + if cpc_nfaces is 4: + rec_radius_ref = np.minimum(rec_w, rec_l)/2. + rec_radius = np.maximum(rec_w, rec_l)/2. + else: + rec_radius = np.sqrt(rec_w*rec_w+rec_l*rec_l)/2. + rec_radius_ref = rec_radius + + return rec_radius_ref, rec_radius + + def hyperboloidparameters(self, rec_z, cpc_theta_deg, cpc_h, field_rim_angle, aim_z): + ''' + ''' + #intersection point of line passing by aiming point and furthest heliostat & line passing by secondary real foci (aperture of cpc) forming an angle theta with the vertical (cpc acceptance angle) + theta = cpc_theta_deg*np.pi/180. + phy = field_rim_angle*np.pi/180. + x_inter = ((rec_z+cpc_h) - aim_z) / (1/np.tan(phy) - 1/np.tan(theta)) + z_inter = x_inter/np.tan(phy)+aim_z + # translate the intersection point to place the mid distance of teh 2 focci at the origin + foci_dist = (aim_z - (rec_z+cpc_h))/2. + z_inter -= foci_dist + + assert z_inter>0., 'there is no corresponding hyperbola to fit the design criteria' + + foci_dist_sqrt = foci_dist**2 + z_inter_sqrt = z_inter**2 + x_inter_sqrt = x_inter**2 + # c0, c1, c2 are the hyperbola equation parameters with the distance to the vertex as the unknown parameter + c2 = 1. + c1 = -(z_inter_sqrt + x_inter_sqrt + foci_dist_sqrt) + c0 = z_inter_sqrt*foci_dist_sqrt + delta = c1**2 - 4*c0*c2 + vertex_dist_sqrt = (-c1-np.sqrt(delta))/2. + check = z_inter_sqrt/vertex_dist_sqrt - x_inter_sqrt/(foci_dist_sqrt-vertex_dist_sqrt) + print('equation resolution should be 1, and it is: ', check) + vertex_dist = np.sqrt(vertex_dist_sqrt) + + return vertex_dist, x_inter, foci_dist + + def intersectionlinehyperbol(self, a_hyper, b_hyper, m_line, z0_line): + ''' + intersection point between hyperbola and line: + z^2/a^2 - y^2/b^2 = 1 + z = m*y + z0 + ''' + ### Previous method if the coefficient of the line is outside the tangent of the hyperbola, there will be no intersection dummy! + # hyper_tangent = a_hyper/b_hyper + # print('tangent: ', hyper_tangent) + # print('m_line', m_line) + # if abs(m_line) <= hyper_tangent: + # m_line = hyper_tangent * 1.05 + + b_hyper_sqrt = b_hyper**2 + a_hyper_sqrt = a_hyper**2 + m_line_sqrt = m_line**2 + aCoef_poly = b_hyper_sqrt - a_hyper_sqrt/m_line_sqrt + bCoef_poly = 2*z0_line*a_hyper_sqrt/m_line_sqrt + cCoef_poly = -a_hyper_sqrt*((z0_line**2)/m_line_sqrt+b_hyper_sqrt) + delta = bCoef_poly**2 - 4*aCoef_poly*cCoef_poly + root1 = (-bCoef_poly+np.sqrt(delta))/(2*aCoef_poly) + root2 = (-bCoef_poly-np.sqrt(delta))/(2*aCoef_poly) + print('root1: ', root1) + print('root2: ', root2) + z_inter = max(root1,root2) + y_inter = (z_inter-z0_line)/m_line + check = (z_inter**2)/a_hyper_sqrt - (y_inter**2)/b_hyper_sqrt + print('equation resolution should be 1, and it is: ', check) + + return y_inter + + def zhyperboloid(self, x_hyper, y_hyper, vertex_dist, foci_dist): + ''' + calculate the z coordinate of the hyperboloid: + (x^2+y^2)/a^2 - (z+z0-g/2)^2/b^2 + 1 = 0 + with the vetex of the hyperboloid located in the xOy plan. + ''' + g_para = 2*foci_dist + f_para = (vertex_dist+foci_dist)/g_para + a_para_sqrt = (g_para**2)*(f_para-f_para**2) + b_para = g_para*(f_para-1/2.) + z0_hyper = abs(b_para) + g_para/2. + + z_hyper = b_para*np.sqrt((x_hyper**2+y_hyper**2)/a_para_sqrt+1) - z0_hyper+g_para/2. + return z_hyper + + def receiversystem(self, receiver, rec_abs=1., rec_w=1.2, rec_l=10., rec_z=0., rec_grid=200, cpc_nfaces=4, cpc_theta_deg=20., cpc_h=None, + cpc_nZ=20., field_rim_angle=30., aim_z=62., secref_fratio=None, refl_sec=0.95, slope_error=0.0, secref_vert=np.array([[-15,25],[-15,-25],[15,-25],[15,25]])): + + ''' + Variables: + - cpc_theta_deg + - field_rim_angle + - aim_z (tower height) + - rec_z (optional) + Arguments: + (1) receiver : str, type of the receiver; 'beam-down', + # Arguments for flat receiver + (2) rec_abs : float, receiver surface absorptivity, e.g. 0.9 + (3) rec_w : float, width of the receiver (m) + (4) rec_l : float, length of the receiver (m) + (5) rec_z : float, z (vertical) location of the receiver (m) + (6) rec_grid : number of slices for solstice flux map + # Arguments for Compound Parabolic Concentrator (CPC) + (7) cpc_nfaces : int, number of faces of the CPC + (8) cpc_theta_deg : float, acceptance angle of CPC (deg) + (9) cpc_h : float, vertical distance of the maximum and minimum point of the parabola + in local coordinate system of the parabola in the xOy plan + WARNING: cpc_h is different from the total height of the CPC + (10) cpc_nZ : int, number of number of incrementation for the clipping polygon for the construction of each CPC face + # Arguments for Secondary Reflector + (11) field_rim_angle : float, rim angle of the field formed wih vertical axis in the zOx plan (deg) + (12) aim_z : float, z (vertical) location of the heliostats' aiming point (m) + (13) secref_fratio : flaot, ratio of the foci distance and apex distance to the origin [1/2,1] + (14) r_f : float, secondary mirror and CPC reflectivity, e.g. 0.9 + (15) slope_error : float, slope error of secondary mirror and CPC refelctivity (?) + (16) secref_vert : array, clipping polygon of the secondary reflector + ''' + self.receiver=receiver + self.rec_abs=rec_abs + + assert cpc_nfaces > 2, f'The number of faces for the CPC should be minimum 3, and it is {cpc_nfaces=}' + if secref_fratio != None: + assert 0.5<=secref_fratio<1, f'The ratio of the hyperbole distances to foci and apex should be between 0.5 and 1, and it is {secref_fratio=}' + rec_rad_ref, rec_rad = self.cpcradius(rec_w, rec_l, cpc_nfaces) + + # Calculate the maximum height of the CPC if not defined + if cpc_h is None: + cpc_theta = cpc_theta = cpc_theta_deg*np.pi/180. + cpc_h = rec_rad_ref*(1+1/np.sin(cpc_theta))/np.tan(cpc_theta) + print('height of cpc: ', cpc_h) + + assert aim_z > (cpc_h+rec_z), 'The imaginary foci of the hyperbol is lower than its real foci' + + # Calculate the characteristics of the secondary reflector + if secref_fratio is None: + vertex_dist, x_inter, foci_dist = self.hyperboloidparameters(rec_z, cpc_theta_deg, cpc_h, field_rim_angle, aim_z) + else: + foci_dist = (aim_z-(cpc_h+rec_z))/2. + vertex_dist = secref_fratio*foci_dist + x_inter = None + secref_z = rec_z + cpc_h + vertex_dist + foci_dist + print('hyperboloid ratio: ', vertex_dist/foci_dist) + + # Calculate the clipping polygon of the secondary reflector + if secref_vert is None: + b_hyper = np.sqrt(foci_dist**2-vertex_dist**2) + if x_inter is None: + m_line = -1/np.tan(field_rim_angle*np.pi/180.) + x_inter = self.intersectionlinehyperbol(vertex_dist, b_hyper, m_line, foci_dist) + + if rec_rad_ref != rec_rad: + m_line = cpc_h/rec_rad + # if the coefficient of the line is outside the tangent of the hyperbola, there will be no intersection + hyper_tangent = vertex_dist/b_hyper + if abs(m_line) <= hyper_tangent: + y_inter = rec_rad + (x_inter-rec_rad_ref) + else: + y_inter = self.intersectionlinehyperbol(vertex_dist, b_hyper, m_line, -foci_dist) + else: + y_inter = x_inter + print('x_inter: ', x_inter) + print('THEY SHOULD BE EQUAL:', y_inter) + secref_vert = np.array([[-x_inter,y_inter],[-x_inter,-y_inter],[x_inter,-y_inter],[x_inter,y_inter]]) + + # calculate the field maxium radius along x and y axis + self.x_max = aim_z*np.tan(field_rim_angle*np.pi/180.) + y_inter = max(secref_vert[:,1]) + print('THEY SHOULD BE EQUAL:', y_inter) + z_hyper = secref_z + self.zhyperboloid( 0.0, y_inter, vertex_dist, foci_dist) + phy = np.arctan(y_inter/(aim_z-z_hyper)) + self.y_max = aim_z*np.tan(phy) + + self.rec_param=np.array([rec_w, rec_l, rec_z, rec_grid, cpc_nfaces, cpc_theta_deg, cpc_h, cpc_nZ, aim_z, secref_z, refl_sec, slope_error, secref_vert]) + + def heliostatfield(self, field, hst_rho, slope, hst_w, hst_h, tower_h, tower_r=0.01, hst_z=0., num_hst=0., R1=0., fb=0., dsep=0., x_max=-1.0, y_max=-1.0): + + ''' + Arguments: + (1) field : str, + -- 'polar', 'polar-half', 'surround' or 'surround-half' for desiging a new field + -- the directory of the layout file + the layout file is a 'csv' file, (n+2, 7) + - n is the total number of heliostats + - the 1st row is the number of each column + - the 2nd row is the unit + - the first three columns are X, Y, Z coordinates of each heliostat + - the fourth column is the focal length + - the last three columns are the aiming points + (2) hst_w : float, width of a heliostat (m) + (3) hst_h : float, height of a heliostat (m) + (4) hst_z : float, the installation height of the heliostat + (5) hst_rho : float, reflectivity of heliostat + (6) slope : float, slope error(radians) + (7) R1 : float, layout parameter, the distance of the first row of heliostat + (8) dsep : float, layout parameter, the separation distance of heliostats (m) + (9) tower_h : float, tower height (m) + (10)tower_r : float, radius of tower (m) + (11)num_hst : int, number of heliostats used in the field design + (12)Q_in_rcv : int, required heat of the receiver in the field design + ''' + + if field[-3:]=='csv': + print('KNOWN FIELD') + layout=np.loadtxt(field, delimiter=',', skiprows=2) + + else: + # design a new field + savefolder=self.casedir+'/des_point' + if not os.path.exists(savefolder): + os.makedirs(savefolder) + + if (x_max>0) and (y_max>0): + x_max = self.x_max + y_max = self.y_max + r_max = np.sqrt(x_max**2 + y_max**2) + pos_and_aiming=radial_stagger(latitude=self.latitude, num_hst=num_hst, width=hst_w, height=hst_h, hst_z=hst_z, towerheight=tower_h, R1=R1, fb=fb, dsep=0., field=field, savedir=savefolder, plot=False, r_field_max= r_max) + else: + pos_and_aiming=radial_stagger(latitude=self.latitude, num_hst=num_hst, width=hst_w, height=hst_h, hst_z=hst_z, towerheight=tower_h, R1=R1, fb=fb, dsep=0., field=field, savedir=savefolder, plot=False) + + layout=pos_and_aiming[2:, :] + + self.hst_w=hst_w + self.hst_h=hst_h + self.hst_rho=hst_rho + self.slope=slope + self.tower_h=tower_h + self.tower_r=tower_r + + self.hst_pos=layout[:,:3].astype(float) + self.hst_foc=layout[:,3].astype(float) + self.hst_aims=layout[:,4:7].astype(float) + self.hst_row=layout[:,-1].astype(float) + + Xmax = max(self.hst_pos[:,0]) + if (x_max>R1) and (Xmax>x_max*1.0): + select_hst=np.array([]) + for i in range(len(self.hst_pos)): + Xhelio = self.hst_pos[i,0] + if np.abs(Xhelio)R1) and (Ymax>y_max*1.0): + select_hst=np.array([]) + for i in range(len(self.hst_pos)): + Yhelio = self.hst_pos[i,1] + if np.abs(Yhelio)0: + hemisphere='North' + else: + hemisphere='South' + gen_yaml(sun, self.hst_pos, self.hst_foc, self.hst_aims, self.hst_w + , self.hst_h, self.hst_rho, self.slope, self.receiver, self.rec_param + , self.rec_abs, outfile_yaml=outfile_yaml, outfile_recv=outfile_recv + , hemisphere='North', tower_h=self.tower_h, tower_r=self.tower_r + , spectral=False , medium=att_factor, one_heliostat=False) + + + def field_design_annual(self, dni_des, num_rays, nd, nh, weafile, method, Q_in_des=None, n_helios=None, zipfiles=False, gen_vtk=False, plot=False): + ''' + Design a field according to the ranked annual performance of heliostats + (DNI weighted) + + ''' + print('') + print('Start field design') + dni_weight=self.dni_TMY(weafile, nd, nh) + + AZI, ZENITH,table,case_list=self.sun.annual_angles(self.latitude, casefolder=self.casedir, nd=nd, nh=nh) + case_list=case_list[1:] + SOLSTICE_AZI, SOLSTICE_ELE=self.sun.convert_convention('solstice', AZI, ZENITH) + + oelt=table + run=np.r_[0] + nhst=len(self.hst_pos) + ANNUAL=np.zeros(nhst) + hst_annual={} + + for i in range(len(case_list)): + c=int(case_list[i,0].astype(float)) + if c not in run: + azimuth=SOLSTICE_AZI[c-1] + elevation= SOLSTICE_ELE[c-1] + + if np.sin(elevation*np.pi/180.)>=1.e-5: + dni=1618.*np.exp(-0.606/(np.sin(elevation*np.pi/180.)**0.491)) + else: + dni=0. + + sys.stderr.write("\n"+green('Sun position: %s \n'%c)) + print('azimuth: %.2f'% azimuth, ', elevation: %.2f'%elevation) + + onesunfolder=os.path.join(self.casedir,'sunpos_%s'%(c)) + + # run solstice + if elevation<1.: + efficiency_total=0 + performance_hst=np.zeros((nhst, 9)) + efficiency_hst=np.zeros(nhst) + else: + efficiency_total, performance_hst=self.master.run(azimuth, elevation, num_rays, self.hst_rho, dni, folder=onesunfolder, gen_vtk=False, printresult=False, system='beamdown') + + efficiency_hst=performance_hst[:,-1]/performance_hst[:,0] + + + hst_annual[c]=performance_hst + sys.stderr.write(yellow("Total efficiency: {:f}\n".format(efficiency_total))) + + cc=0 + for a in range(len(table[3:])): + for b in range(len(table[0,3:])): + val=re.findall(r'\d+', table[a+3,b+3]) + if val==[]: + table[a+3,b+3]=0 + else: + if c==float(val[0]): + #table[a+3,b+3]=efficiency_total # this line cause wired problem + if cc==0: # avoid adding the symetrical point + ANNUAL+=dni_weight[a,b]*efficiency_hst + cc+=1 + run=np.append(run,c) + np.savetxt(self.casedir+'/annual_hst.csv',ANNUAL, fmt='%.2f', delimiter=',') + + designfolder=self.casedir+'/des_point' + day=self.sun.days(21, 'Mar') + dec=self.sun.declination(day) + hra=0. # solar noon + zen=self.sun.zenith(self.latitude, dec, hra) + azi=self.sun.azimuth(self.latitude, zen, dec, hra) + azi_des, ele_des=self.sun.convert_convention('solstice', azi, zen) + + sys.stderr.write("\n"+green('Design Point: \n')) + efficiency_total, performance_hst_des=self.master.run(azi_des, ele_des, num_rays, self.hst_rho, dni_des, folder=designfolder, gen_vtk=True, printresult=False, system='beamdown') + + Qin=performance_hst_des[:,-1] + Qsolar=performance_hst_des[0,0] + + ID=ANNUAL.argsort() + ID=ID[::-1] + + select_hst=np.array([]) + rmax=np.max(self.hst_row) + if method==1: + print('') + print('Method 1') + power=0. + self.Q_in_rcv=Q_in_des + for i in range(len(ID)): + if power0.1*rmax: + select_hst=np.append(select_hst, idx) + power+=Qin[idx] + + else: + print('') + print('Method 2') + num_hst=0 + power=0. + for i in range(len(ID)): + if num_hst= 0.0001: + # Calculate the local y coordinate + z_para += delta_Z + y_para = np.sqrt(4*focal_length*z_para) + # Calculate the radius of the CPC at the height z_para + _, y_CPC, z_CPC = self.thetagammarotations(y_para, z_para, theta, 0.0) + z_CPC += z_trans + if z_CPC <= self.h_CPC + 0.000001: + # Calculate the local/global x coordinate + x_para = (y_CPC+y_trans)*increment+margin + poly_vertex[i]=[x_para, y_para] + poly_vertex[2*cpc_nZ+1-i]=[-x_para, y_para] + i += 1 + else: + z_para -= delta_Z + delta_Z/=2 + + return poly_vertex + + def dependantparameters(self, rec_param): + ''' + Return the clipping polygon (vertex) of the CPC face + ''' + rec_w=rec_param[0] + rec_l=rec_param[1] + rec_z=rec_param[2] + cpc_nfaces=int(rec_param[4]) + cpc_theta_deg=rec_param[5] + cpc_h=rec_param[6] + + cpc_nfaces = int(cpc_nfaces) + if cpc_nfaces==4: + self.rec_x = rec_w/2. + self.rec_y = rec_l/2. + self.rec_radius = np.minimum(self.rec_x,self.rec_y) + else: + self.rec_x = np.sqrt(rec_w*rec_w+rec_l*rec_l)/2. + self.rec_y = self.rec_x + self.rec_radius = self.rec_x + + + # Parameters of the parabola + self.theta = cpc_theta_deg*np.pi/180. + self.focal_length = self.rec_radius*(np.sin(self.theta) + 1) + + # Lowest and highest point of the parabola curve of the CPC face given in local coordinates of the parabola + p_A = 1. + p_B = 4*self.focal_length*np.tan(self.theta) + p_C = -4*self.focal_length*self.focal_length + self.y_min_para, self.z_min_para = self.intersectionpoint(p_A, p_B, p_C, self.focal_length) + ''' OLD + if cpc_h<=0.: + p_B = -4.0*self.focal_length/np.tan(2.*self.theta) + self.y_max_para, self.z_max_para = self.intersectionpoint(p_A, p_B, p_C, self.focal_length) + else: + self.z_max_para = self.z_min_para+cpc_h + self.y_max_para = np.sqrt(4*self.focal_length*self.z_max_para) + # height of the CPC + _, _, z_trans = self.xyztranslations(self.theta, 0.0, 0.0, self.rec_radius, 0.0) + _, _, z_CPC = self.thetagammarotations(self.y_max_para, self.z_max_para, self.theta, 0.0) + self.h_CPC = z_trans + z_CPC + print('calculated height: ', self.h_CPC) + h_CPC = self.rec_radius*(1+1/np.sin(self.theta))/np.tan(self.theta) + print('theoretical critical height: ', h_CPC) + ''' + h_CPC = self.rec_radius*(1+1/np.sin(self.theta))/np.tan(self.theta) + p_B = -4.0*self.focal_length/np.tan(2.*self.theta) + _, self.z_max_para = self.intersectionpoint(p_A, p_B, p_C, self.focal_length) + if cpc_h<=0.: + self.h_CPC = h_CPC + else: + self.h_CPC = cpc_h + + print('theoretical critical height: ', h_CPC) + print('selected height: ', self.h_CPC) + + def beamdowncomponents(self, rec_param): + ''' + Generate YAML files for the Solstice simulation + + Arguments: + # Receiver (flat) + (1) rec_w : float, width of the receiver (m) + (2) rec_l : float, length of the receiver (m) + (3) rec_z : float, z (vertical) location of the receiver (m) + (4) rec_grid : number of slices for solstice flux map + # Compound Parabolic Concentrator (CPC) + (5) cpc_nfaces : int, number of faces of the CPC + (6) cpc_theta_deg : float, acceptance angle of CPC (deg) + (7) cpc_h : float, vertical distance of the maximum and minimum point of the parabola + in local coordinate system of the parabola in the xOy plan + WARNING: cpc_h is different from the total height of the CPC + (8) cpc_nZ : int, number of number of incrementation for the clipping polygon for the construction of each CPC face + # Secondary Reflector (hyperboloid) + (9) aim_z : float, z (vertical) location of the heliostats' aiming point (m) + (10) secref_z : flaot, z (vertical) location of the apex of the hyperboloid secondary reflector (m) + (11) r_f : float, secondary mirror and CPC reflectivity, e.g. 0.9 + (12) slope_error : float, slope error of secondary mirror and CPC refelctivity (?) + (13) secref_vert : array, (x,y) clipping polygon of the secondary reflector + ''' + + rec_z=rec_param[2] + rec_grid=int(rec_param[3]) # it assumes equal number of slices in x and y directions + cpc_nfaces=int(rec_param[4]) + cpc_theta_deg=rec_param[5] + cpc_nZ=int(rec_param[7]) + aim_z=rec_param[8] + secref_z=rec_param[9] + r_f=rec_param[10] # front + slope_error=rec_param[11] + secref_vert=rec_param[12] + self.dependantparameters(rec_param[:7]) + + iyaml='\n' + + # + # Materials + # + # CREATE a specular material for the secondary reflector + iyaml+='- material: &%s\n' % 'material_sec_mirror' + iyaml+=' front:\n' + iyaml+=' mirror: {reflectivity: %6.4f, slope_error: %15.8e }\n' % (r_f, slope_error) + iyaml+=' back:\n' + iyaml+=' matte: {reflectivity: 0.0 }\n' + iyaml+='\n' + + # + # Entities + # + # Receiver Polygon + rec_vert=self.regularpolygon(self.rec_x, self.rec_y, cpc_nfaces, 1) + iyaml+='- entity:\n' + iyaml+=' name: %s\n' % 'receiver' + iyaml+=' primary: 0\n' + iyaml+=' transform: { translation: %s, rotation: %s }\n' % ([0, 0, rec_z], [0, 0, 0]) + iyaml+=' geometry:\n' + iyaml+=' - material: *%s\n' % 'material_target' + iyaml+=' plane:\n' + iyaml+=' slices: %d\n' % rec_grid + iyaml+=' clip: \n' + iyaml+=' - operation: AND \n' + iyaml+=' vertices: \n' + for i in range(cpc_nfaces): + iyaml+=' - %s\n' % ([rec_vert[i][0],rec_vert[i][1]]) + iyaml+='\n' + # + # CREATE a virtual target entity at the receiver level + pts = np.array([ [-1.0, 1.0], [-1.0, -1.0], [1.0, -1.0], [1.0, 1.0] ]) + pts = pts*(self.rec_radius * 50.0) + slices = 4 + iyaml+='\n- entity:\n' + iyaml+=' name: virtual_target\n' + iyaml+=' primary: 0\n' + iyaml+=' transform: { translation: %s, rotation: %s }\n' % ([0, 0, -5], [0., 0, 0]) + iyaml+=' geometry: \n' + iyaml+=' - material: *%s\n' % 'material_virtual' + iyaml+=' plane: \n' + iyaml+=' slices: %d\n' % slices + iyaml+=' clip: \n' + iyaml+=' - operation: AND \n' + iyaml+=' vertices: \n' + for i in range(len(pts)): + iyaml+=' - %s\n' % ([pts[i][0],pts[i][1]]) + iyaml+='\n' + # + # Secondary Reflector + # focal image: distance between hyperbola vertex and aiming point of heliostats + # focal real: distance between hyperbola vertex and aiming point of secondary reflector (aperture of the CPC) + focal_image = aim_z - secref_z + focal_real = abs(rec_z + self.h_CPC - secref_z) # aim at the aperture of the CPC + iyaml+='- entity:\n' + iyaml+=' name: %s\n' % 'secondary_reflector' + iyaml+=' primary: 0\n' + iyaml+=' transform: { translation: %s, rotation: %s }\n' % ([0, 0, secref_z], [0., 0, 0]) + iyaml+=' geometry:\n' + iyaml+=' - material: *%s\n' % 'material_sec_mirror' + iyaml+=' hyperbol:\n' + iyaml+=' focals: &hyperbol_focals { real: %s, image: %s }\n' % (focal_real, focal_image) + iyaml+=' clip: \n' + iyaml+=' - operation: AND \n' + iyaml+=' vertices: \n' + for i in range(len(secref_vert)): + iyaml+=' - %s\n' % ([secref_vert[i][0],secref_vert[i][1]]) + iyaml+='\n' + # + # CREATE a virtual target entity above the secondary reflector + pts = pts*5 + iyaml+='\n- entity:\n' + iyaml+=' name: virtual_sec_ref\n' + iyaml+=' primary: 0\n' + iyaml+=' transform: { translation: %s, rotation: %s }\n' % ([0, 0, aim_z+50], [0., 180., 0]) + iyaml+=' geometry: \n' + iyaml+=' - material: *%s\n' % 'material_virtual' + iyaml+=' plane: \n' + iyaml+=' slices: %d\n' % slices + iyaml+=' clip: \n' + iyaml+=' - operation: AND \n' + iyaml+=' vertices: \n' + for i in range(len(pts)): + iyaml+=' - %s\n' % ([pts[i][0],pts[i][1]]) + iyaml+='\n' + # + # CREATE the receiver objects in receiver yaml + # receiver at the CPC exit and 2 virtuals surfaces + rcv = '\n' + rcv+='- name: receiver \n' + rcv+=' side: %s \n' % 'FRONT_AND_BACK' + rcv+=' per_primitive: %s \n' % 'INCOMING_AND_ABSORBED' + rcv+='- name: virtual_sec_ref \n' + rcv+=' side: %s \n' % 'FRONT' + rcv+=' per_primitive: %s \n' % 'INCOMING' + rcv+='- name: virtual_target \n' + rcv+=' side: %s \n' % 'FRONT' + rcv+=' per_primitive: %s \n' % 'INCOMING' + # + # CREATE the facets of the CPC + # For each face, the x,y,z translations and the clipping polygon are calculated + # The clipping polygon is individual to a face if nCPC = 4 and the receiver is not a regular polygon (xRec!=yRec) + cpc_face_pos=self.regularpolygon(self.rec_x, self.rec_y, cpc_nfaces, 0) + cpc_face_polygon=self.cpcfaceclipping(cpc_nZ, self.theta, self.focal_length, cpc_nfaces, self.rec_radius) + rec_vert = np.insert(rec_vert,0,rec_vert[cpc_nfaces-1],axis=0) + increment = np.zeros(cpc_nfaces) + if self.rec_x != self.rec_y: + for i in range(cpc_nfaces): + increment[i] = np.linalg.norm(rec_vert[i]-rec_vert[i+1])/2 - self.rec_radius + # implement individual facets of the CPC + for i in range(cpc_nfaces): + gamma = i*2.*np.pi/cpc_nfaces + iyaml+='- entity:\n' + iyaml+=' name: %s\n' % ('CPC_face_'+str(i)) + x_trans, y_trans, z_trans = self.xyztranslations(self.theta, gamma, cpc_face_pos[i][0], cpc_face_pos[i][1], rec_z) + iyaml+=' transform: { translation: %s }\n' % ([x_trans, y_trans, z_trans]) + iyaml+=' children: \n' + iyaml+=' - name: RotationGamma\n' + iyaml+=' transform: { rotation: %s } \n' % ([0, 0, gamma*180.0/np.pi]) + iyaml+=' children: \n' + iyaml+=' - name: RotationTheta\n' + iyaml+=' primary: 0\n' + iyaml+=' transform: {rotation: %s } \n' % ([cpc_theta_deg, 0, 0.]) + iyaml+=' geometry:\n' + iyaml+=' - material: *%s\n' % 'material_sec_mirror' + iyaml+=' parabolic-cylinder:\n' + iyaml+=' focal: %s\n' % self.focal_length + iyaml+=' clip: \n' + iyaml+=' - operation: AND \n' + iyaml+=' vertices: \n' + midlength = int(len(cpc_face_polygon)/2) + for j in range(midlength): + iyaml+=' - %s\n' % ([cpc_face_polygon[j][0]+increment[i],cpc_face_polygon[j][1]]) + for j in range(midlength,midlength*2): + iyaml+=' - %s\n' % ([cpc_face_polygon[j][0]-increment[i],cpc_face_polygon[j][1]]) + iyaml+='\n' + + + return iyaml, rcv + +if __name__=='__main__': + start=time.time() +# casedir='./test-CPC-design' + +# enter the parameters + rec_w=11. + rec_l=4. + rec_z=0. + rec_grid=20. + n_CPC_faces=4 + theta_deg=20. + cpc_h=0. + n_Z=20 + aim_z = 30. + sec_z = 22. + r_f= 1. # front + slope_error = 0. + secref_vert = np.array([[-15,25],[-15,-10],[15,-10],[15,25]]) + + cpc=CPC() + cpc.beamdowncomponents(rec_w, rec_l, rec_z, rec_grid, n_CPC_faces, theta_deg, cpc_h, n_Z, aim_z, sec_z, r_f, slope_error, secref_vert) + + eta=1. + print('total efficiency:', eta) + + end=time.time() + print('total time %.2f'%((end-start)/60.), 'min') diff --git a/solsticepy/design_crs.py b/solsticepy/design_crs.py index c1560d6..3a7d3da 100644 --- a/solsticepy/design_crs.py +++ b/solsticepy/design_crs.py @@ -21,15 +21,15 @@ class CRS: ''' - The Central Receiver System (CRS) model includes three parts: + The Central Receiver System (CRS) model includes three parts: the sun, the field and the receiver. ''' def __init__(self, latitude, casedir, nproc=None, verbose=False): ''' Arguements: - casedir : str, the directory of the case - nproc (int): number of processors, e.g. nproc=1 will run in serial mode, + casedir : str, the directory of the case + nproc (int): number of processors, e.g. nproc=1 will run in serial mode, nproc=4 will run with 4 processors in parallel nproc=None will run with any number of processors that are available verbose : bool, write results to files or not @@ -47,18 +47,18 @@ def receiversystem(self, receiver, rec_w=0., rec_h=0., rec_x=0., rec_y=0., rec_z ''' Arguements: - (1) receiver : str, type of the receiver, i.e. 'flat', 'cylinder', 'multi_aperture' or directory of the 'stl', - (2) rec_w : float, width of a flat receiver or radius of a cylindrical receiver (m) + (1) receiver : str, type of the receiver, i.e. 'flat', 'cylinder', 'multi_aperture' or directory of the 'stl', + (2) rec_w : float, width of a flat receiver or radius of a cylindrical receiver (m) (3) rec_h : float, height of the receiver (m) (4) rec_x : float, x location of the receiver (m) (5) rec_y : float, y location of the receiver (m) - (6) rec_z : float, z location of the receiver (m), or a list (for multi-aperture) of elevation heights of apertures + (6) rec_z : float, z location of the receiver (m), or a list (for multi-aperture) of elevation heights of apertures (7) rec_tilt : float, tilt angle of the receiver (deg), 0 is where the receiver face to the horizontal (8) rec_grid_w : int, number of elements in the horizontal(x)/circumferential direction (9) rec_grid_h : int, number of elements in the vertical(z) direction (10) rec_abs : float, receiver surface absorptivity, e.g. 0.9 (11) num_aperture : int, number of apertures if it is a multi-aperture receiver - (12) gamma : float, the anangular range of the multi-aperture configration (deg) + (12) gamma : float, the anangular range of the multi-aperture configration (deg) ''' self.receiver=receiver self.rec_abs=rec_abs @@ -71,16 +71,16 @@ def receiversystem(self, receiver, rec_w=0., rec_h=0., rec_x=0., rec_y=0., rec_z self.rec_param=[rec_w, rec_h, receiver, rec_x, rec_y, rec_z, rec_tilt] elif receiver=='flat': self.rec_param=[rec_w, rec_h, rec_grid_w, rec_grid_h, rec_x, rec_y, rec_z, rec_tilt] - elif receiver=='cylinder': - self.rec_param=[rec_w, rec_h, rec_grid_w, rec_grid_h, rec_x, rec_y, rec_z, rec_tilt] - elif receiver=='multi-aperture': + elif receiver=='cylinder': + self.rec_param=[rec_w, rec_h, rec_grid_w, rec_grid_h, rec_x, rec_y, rec_z, rec_tilt] + elif receiver=='multi-aperture': # rec_w and rec_h is the size of one aperture # rec_grid_w and rec_gird_h is the number of elements of one aperture # rec_x, rec_y, rec_z is the location of the central point of the multi-aperture receiver # rec_tilt is the tilt angle of each aperture # num_aperture is the number of apertures - # gamma is the anangular range of the multi-aperture configration (deg) - self.rec_param=[rec_w, rec_h, rec_grid_w, rec_grid_h, rec_z, rec_tilt, num_aperture, gamma] + # gamma is the anangular range of the multi-aperture configration (deg) + self.rec_param=[rec_w, rec_h, rec_grid_w, rec_grid_h, rec_z, rec_tilt, num_aperture, gamma] self.rec_w=rec_w @@ -90,40 +90,40 @@ def heliostatfield(self, field, hst_rho, slope, hst_w, hst_h, tower_h, tower_r=0 ''' Arguements: (1) field : str, - -- 'polar', 'polar-half', 'surround' or 'surround-half' or 'multi-aperture' for desiging a new field + -- 'polar', 'polar-half', 'surround' or 'surround-half' or 'multi-aperture' for desiging a new field -- or the directory of the layout file the layout file is a 'csv' file, (n+2, 7) - - n is the total number of heliostats + - n is the total number of heliostats - the 1st row is the number of each column - - the 2nd row is the unit + - the 2nd row is the unit - the first three columns are X, Y, Z coordinates of each heliostat - the fourth column is the focal length - the last three columns are the aiming points - (2) hst_w : float, width of a heliostat (m) + (2) hst_w : float, width of a heliostat (m) (3) hst_h : float, height of a heliostat (m) (4) hst_z : float, the installation height of the heliostat - (5) hst_rho : float, reflectivity of heliostat + (5) hst_rho : float, reflectivity of heliostat (6) slope : float, slope error(radians) - (7) R1 : float, layout parameter, the distance of the first row of heliostat + (7) R1 : float, layout parameter, the distance of the first row of heliostat (8) dsep : float, layout parameter, the separation distance of heliostats (m) (9) tower_h : float, tower height (m) (10)tower_r : float, radius of tower (m) - (11)num_hst : int, number of heliostats used in the field design + (11)num_hst : int, number of heliostats used in the field design ''' - + if field[-3:]=='csv': print('KNOWN FIELD') layout=np.loadtxt(field, delimiter=',', skiprows=2) else: - # design a new field + # design a new field savefolder=self.casedir+'/des_point' if not os.path.exists(savefolder): os.makedirs(savefolder) - pos_and_aiming, self.Nzones, self.Nrows=radial_stagger(latitude=self.latitude, num_hst=num_hst, width=hst_w, height=hst_h, hst_z=hst_z, towerheight=tower_h, R1=R1, fb=fb, dsep=0., field=field, num_aperture=self.num_aperture, gamma=self.gamma, rec_w=self.rec_w, rec_z=self.rec_z, savedir=savefolder, verbose=self.verb ) - + pos_and_aiming, self.Nzones, self.Nrows=radial_stagger(latitude=self.latitude, num_hst=num_hst, width=hst_w, height=hst_h, hst_z=hst_z, towerheight=tower_h, R1=R1, fb=fb, dsep=0., field=field, num_aperture=self.num_aperture, gamma=self.gamma, rec_w=self.rec_w, rec_z=self.rec_z, savedir=savefolder, verbose=self.verb ) + layout=pos_and_aiming[2:, :] @@ -131,11 +131,11 @@ def heliostatfield(self, field, hst_rho, slope, hst_w, hst_h, tower_h, tower_r=0 self.hst_h=hst_h self.hst_rho=hst_rho self.slope=slope - self.tower_h=tower_h + self.tower_h=tower_h self.tower_r=tower_r self.hst_pos=layout[:,:3].astype(float) - self.hst_foc=layout[:,3].astype(float) + self.hst_foc=layout[:,3].astype(float) self.hst_aims=layout[:,4:7].astype(float) self.hst_aim_idx=layout[:,7].astype(float) @@ -169,12 +169,12 @@ def yaml(self, dni=1000,sunshape=None,csr=0.01,half_angle_deg=0.2664,std_dev=0.2 def field_design_annual(self, dni_des, num_rays, nd, nh, weafile, method, Q_in_des=None, n_helios=None, zipfiles=False, gen_vtk=False, plot=False): ''' - Design a field according to the ranked annual performance of heliostats + Design a field according to the ranked annual performance of heliostats (DNI weighted) - ''' + ''' print('') - print('Start field design') + print('Start field design') system=self.receiver AZI, ZENITH,table,case_list=self.sun.annual_angles(self.latitude, casefolder=self.casedir,nd=nd, nh=nh, verbose=self.verb) @@ -182,14 +182,14 @@ def field_design_annual(self, dni_des, num_rays, nd, nh, weafile, method, Q_in_ SOLSTICE_AZI, SOLSTICE_ELE=self.sun.convert_convention('solstice', AZI, ZENITH) run=np.r_[0] - nhst=len(self.hst_pos) + nhst=len(self.hst_pos) #oelt=table - ANNUAL=np.zeros(nhst) - annual_solar=0. + ANNUAL=np.zeros(nhst) + annual_solar=0. hst_annual={} - for i in range(len(case_list)): + for i in range(len(case_list)): c=int(case_list[i,0].astype(float)) if c not in run: # the morning positions @@ -209,32 +209,32 @@ def field_design_annual(self, dni_des, num_rays, nd, nh, weafile, method, Q_in_ # run solstice if elevation<1.: efficiency_total=0 - performance_hst=np.zeros((nhst, 9)) + performance_hst=np.zeros((nhst, 9)) efficiency_hst=np.zeros(nhst) else: efficiency_total, performance_hst=self.master.run(azimuth, elevation, num_rays, self.hst_rho, dni, folder=onesunfolder, gen_vtk=gen_vtk, printresult=False, verbose=self.verb, system=system) - + #res=np.loadtxt(onesunfolder+'/result-formatted.csv', dtype=str, delimiter=',') #res_hst=np.loadtxt(onesunfolder+'/heliostats-raw.csv', dtype=str, delimiter=',') #efficiency_total=res[-2,1].astype(float)/res[1,1].astype(float) #performance_hst=res_hst[1:,-9:].astype(float) - + efficiency_hst=performance_hst[:,-1]/performance_hst[:,0] hst_annual[c]=performance_hst sys.stderr.write(yellow("Total efficiency: {:f}\n".format(efficiency_total))) - run=np.append(run,c) + run=np.append(run,c) cc=0 for a in range(len(table[3:])): for b in range(len(table[0,3:])): val=re.findall(r'\d+', table[a+3,b+3]) if str(c) in val: - if cc==0: + if cc==0: # i.e. morning positions ANNUAL+=dni*efficiency_hst annual_solar+=dni - cc+=1 + cc+=1 else: # the symetrical points (i.e. afternoon) eff_symetrical=np.array([]) @@ -251,37 +251,37 @@ def field_design_annual(self, dni_des, num_rays, nd, nh, weafile, method, Q_in_ eff_row=eff_row[::-1] else: eff_row[1:]=eff_row[1:][::-1] - + eff_symetrical=np.append(eff_symetrical, eff_row) #print(np.shape(eff_symetrical)) #check=np.append(self.hst_zone, (self.hst_row, self.hst_num_idx, efficiency_hst, eff_symetrical)) #print(np.shape(check)) #check=check.reshape(5,int(len(check)/5)) - #np.savetxt('./check.csv', check.T, fmt='%.5f', delimiter=',') - ANNUAL+=dni*eff_symetrical - annual_solar+=dni - - ANNUAL/=annual_solar - if self.verb: + #np.savetxt('./check.csv', check.T, fmt='%.5f', delimiter=',') + ANNUAL+=dni*eff_symetrical + annual_solar+=dni + + ANNUAL/=annual_solar + if self.verb: np.savetxt(self.casedir+'/annual_hst.csv',ANNUAL, fmt='%.2f', delimiter=',') - + designfolder=self.casedir+'/des_point' day=self.sun.days(21, 'Mar') dec=self.sun.declination(day) hra=0. # solar noon zen=self.sun.zenith(self.latitude, dec, hra) - azi=self.sun.azimuth(self.latitude, zen, dec, hra) - azi_des, ele_des=self.sun.convert_convention('solstice', azi, zen) + azi=self.sun.azimuth(self.latitude, zen, dec, hra) + azi_des, ele_des=self.sun.convert_convention('solstice', azi, zen) - sys.stderr.write("\n"+green('Design Point: \n')) + sys.stderr.write("\n"+green('Design Point: \n')) efficiency_total, performance_hst_des=self.master.run(azi_des, ele_des, num_rays, self.hst_rho, dni_des, folder=designfolder, gen_vtk=gen_vtk, printresult=False, verbose=self.verb, system=system) - + #res=np.loadtxt(designfolder+'/result-formatted.csv', dtype=str, delimiter=',') #res_hst=np.loadtxt(designfolder+'/heliostats-raw.csv', dtype=str, delimiter=',') #efficiency_total=res[-2,1].astype(float)/res[1,1].astype(float) #performance_hst_des=res_hst[1:,-9:].astype(float) - + Qin=performance_hst_des[:,-1] Qsolar=performance_hst_des[0,0] @@ -298,7 +298,7 @@ def field_design_annual(self, dni_des, num_rays, nd, nh, weafile, method, Q_in_ if method==1: hst_aim_idx=self.hst_aim_idx[ID] - print('') + print('') print('Method 1') self.Q_in_rcv=Q_in_des if self.receiver=='multi-aperture': @@ -327,7 +327,7 @@ def field_design_annual(self, dni_des, num_rays, nd, nh, weafile, method, Q_in_ else: # initial selection - # for single-aperture receiver + # for single-aperture receiver # or multi-aperture receiver configuration that selects heliostats based on the total required heat assert isinstance(Q_in_des, float), "Q_in_des should be float, which is the total required incident power to the receiver" @@ -339,12 +339,12 @@ def field_design_annual(self, dni_des, num_rays, nd, nh, weafile, method, Q_in_ ap_idx=int(hst_aim_idx[i]) self.Q_in_rcv_i[ap_idx]+=Qin[idx] - + else: select_hst=np.array([]) - print('') - print('Method 2') - #TODO the Method 2 does not include multi-aperture option + print('') + print('Method 2') + #TODO the Method 2 does not include multi-aperture option num_hst=0 power=0. for i in range(len(ID)): @@ -356,7 +356,7 @@ def field_design_annual(self, dni_des, num_rays, nd, nh, weafile, method, Q_in_ power+=Qin[idx] self.Q_in_rcv=power - + select_hst=select_hst.astype(int) self.hst_pos= self.hst_pos[select_hst,:] @@ -390,9 +390,9 @@ def field_design_annual(self, dni_des, num_rays, nd, nh, weafile, method, Q_in_ np.savetxt(self.casedir+'/pos_and_aiming.csv', design_pos_and_aim, fmt='%s', delimiter=',') np.savetxt(self.casedir+'/selected_hst.csv', select_hst, fmt='%.0f', delimiter=',') - annual_solar=0. + annual_solar=0. annual_field=0. - + oelt={} QTOT=np.zeros(np.shape(table)) QIN=np.zeros(np.shape(table)) @@ -410,9 +410,9 @@ def field_design_annual(self, dni_des, num_rays, nd, nh, weafile, method, Q_in_ print('Aperture %s'%ap) print('num helios', np.sum(idx_apt_i)) - for i in range(len(case_list)): + for i in range(len(case_list)): c=int(case_list[i,0].astype(float)) - if c not in run: + if c not in run: #sundir=designfolder+'/sunpos_%s'%c res_hst=hst_annual[c] Qtot=res_hst[select_hst,0] @@ -443,8 +443,8 @@ def field_design_annual(self, dni_des, num_rays, nd, nh, weafile, method, Q_in_ QIN[a+3,b+3]+=np.sum(Qin[idx_apt_i]) annual_solar+=dni annual_field+=dni*eff - - + + oelt[ap][2, 3:]=table[2, 3:].astype(float) oelt[ap][3:,2]=table[3:,2].astype(float) @@ -454,7 +454,7 @@ def field_design_annual(self, dni_des, num_rays, nd, nh, weafile, method, Q_in_ if self.num_aperture==1: return oelt[0], A_land else: - oelt[self.num_aperture]= np.divide(QIN, QTOT, out=np.zeros(QIN.shape, dtype=float), where=QTOT!=0) + oelt[self.num_aperture]= np.divide(QIN, QTOT, out=np.zeros(QIN.shape, dtype=float), where=QTOT!=0) oelt[self.num_aperture][2, 3:]=table[2, 3:].astype(float) oelt[self.num_aperture][3:,2]=table[3:,2].astype(float) @@ -465,17 +465,17 @@ def field_design_annual(self, dni_des, num_rays, nd, nh, weafile, method, Q_in_ if self.verb: for ap in range(self.num_aperture+1): if ap==self.num_aperture: - np.savetxt(self.casedir+'/lookup_table_total.csv', oelt[ap], fmt='%s', delimiter=',') + np.savetxt(self.casedir+'/lookup_table_total.csv', oelt[ap], fmt='%s', delimiter=',') else: - np.savetxt(self.casedir+'/lookup_table_%s.csv'%ap, oelt[ap], fmt='%s', delimiter=',') - return oelt, A_land + np.savetxt(self.casedir+'/lookup_table_%s.csv'%ap, oelt[ap], fmt='%s', delimiter=',') + return oelt, A_land def annual_oelt(self, dni_des, num_rays, nd, nh, zipfiles=False, gen_vtk=False, plot=False): ''' Annual performance of a known field - ''' - self.n_helios=len(self.hst_pos) + ''' + self.n_helios=len(self.hst_pos) oelt, ANNUAL=self.master.run_annual(nd=nd, nh=nh, latitude=self.latitude, num_rays=num_rays, num_hst=self.n_helios,rho_mirror=self.hst_rho, dni=dni_des, verbose=self.verb) Xmax=max(self.hst_pos[:,0]) @@ -494,15 +494,15 @@ def annual_oelt(self, dni_des, num_rays, nd, nh, zipfiles=False, gen_vtk=False, dec=self.sun.declination(day) hra=0. # solar noon zen=self.sun.zenith(self.latitude, dec, hra) - azi=self.sun.azimuth(self.latitude, zen, dec, hra) - azi_des, ele_des=self.sun.convert_convention('solstice', azi, zen) + azi=self.sun.azimuth(self.latitude, zen, dec, hra) + azi_des, ele_des=self.sun.convert_convention('solstice', azi, zen) - sys.stderr.write("\n"+green('Design Point: \n')) + sys.stderr.write("\n"+green('Design Point: \n')) efficiency_total, performance_hst_des=self.master.run(azi_des, ele_des, num_rays, self.hst_rho, dni_des, folder=designfolder, gen_vtk=False, printresult=False, verbose=self.verb) self.eff_des=efficiency_total.n return oelt, A_land - + def dni_TMY(self, weafile, nd, nh, plot=False): ''' @@ -512,19 +512,19 @@ def dni_TMY(self, weafile, nd, nh, plot=False): # col0 - time (s) # col1 - GHI # col2 -DNI (W/m2) - # col3 -DHI + # col3 -DHI ... ''' with open(weafile) as f: content=f.read().splitlines() f.close() - + lines=len(content) seconds=np.array([]) dni=np.array([]) for i in range(2, lines-1): - l=content[i].split(",") - seconds=np.append(seconds, l[0]) + l=content[i].split(",") + seconds=np.append(seconds, l[0]) dni=np.append(dni, l[2]) seconds=seconds.astype(float)+1800. days=seconds/3600/24 @@ -543,36 +543,36 @@ def dni_TMY(self, weafile, nd, nh, plot=False): hra_lim=180.+dh/2. dec_lim=23.45+dd/2. - hra_bin=np.linspace(-hra_lim, hra_lim, nh+1) + hra_bin=np.linspace(-hra_lim, hra_lim, nh+1) dec_bin=np.linspace(-dec_lim, dec_lim, nd+1) bins=np.array([hra_bin, dec_bin]) dni_weight, xbins, ybins=np.histogram2d(wea_hra, wea_dec, bins, weights=wea_dni) dni_n, xbins, ybins=np.histogram2d(wea_hra, wea_dec, bins) - + dni_avg=np.divide(dni_weight, dni_n, out=np.zeros_like(dni_weight), where=dni_n!=0) if plot: X=np.linspace(-180., 180. , nh) - Y=np.linspace(-23.45, 23.45, nd) + Y=np.linspace(-23.45, 23.45, nd) plt.pcolormesh(hra_bin, dec_bin, dni_avg.T) cb=plt.colorbar() plt.xlabel('Solar hour angle') plt.ylabel('Delination angle') plt.savefig(self.casedir+'/DNI_weather.png', bbox_inches='tight') - plt.close() + plt.close() - # check the symmetricity of the weather dni data + # check the symmetricity of the weather dni data data=dni_avg.T plt.pcolormesh(hra_bin[:int(nh/2)], dec_bin, data[:,:int(nh/2), ]-np.fliplr(data[:,-int(nh/2):])) cb=plt.colorbar() plt.xlabel('Solar hour angle') plt.ylabel('Delination angle') plt.savefig(self.casedir+'/DNI_weather_diff.png', bbox_inches='tight') - plt.close() + plt.close() np.savetxt(self.casedir+'/DNI_weather.csv', data, fmt='%.4f', delimiter=',') - + return dni_weight.T, dni_avg.T def get_attenuation_factor(self): @@ -604,12 +604,12 @@ def generateVTK(self,eta_hst, savevtk): hst_tot =eta_hst[:,0] rec_abs=eta_hst[:,1] - hst_eff=rec_abs/hst_tot + hst_eff=rec_abs/hst_tot savevtk=savevtk+'/results-field.vtk' - TOT=np.repeat(hst_tot, element) - ABS=np.repeat(rec_abs, element) - EFF=np.repeat(hst_eff, element) + TOT=np.repeat(hst_tot, element) + ABS=np.repeat(rec_abs, element) + EFF=np.repeat(hst_eff, element) DATA={'tot': TOT, 'rec_abs': ABS, 'efficiency':EFF} gen_vtk(savedir=savevtk, points=COORD.T, indices=TRI, norms=NORMS, colormap=True, DATA=DATA) @@ -618,7 +618,7 @@ def generateVTK(self,eta_hst, savevtk): start=time.time() casedir='./test-crs-design' tablefile=casedir+'/OELT_Solstice.motab' - if os.path.exists(tablefile): + if os.path.exists(tablefile): print('') print('Load exsiting OELT') @@ -635,7 +635,7 @@ def generateVTK(self,eta_hst, savevtk): pm.saveparam(casedir) print(pm.fb) print(pm.H_tower) - crs=CRS(latitude=pm.lat, casedir=casedir) + crs=CRS(latitude=pm.lat, casedir=casedir) weafile='/home/yewang/solartherm-master/SolarTherm/Data/Weather/gen3p3_Daggett_TMY3.motab' crs.heliostatfield(field=pm.field_type, hst_rho=pm.rho_helio, slope=pm.slope_error, hst_w=pm.W_helio, hst_h=pm.H_helio, tower_h=pm.H_tower, tower_r=pm.R_tower, hst_z=pm.Z_helio, num_hst=pm.n_helios, R1=pm.R1, fb=pm.fb, dsep=pm.dsep) @@ -646,15 +646,12 @@ def generateVTK(self,eta_hst, savevtk): oelt, A_land=crs.field_design_annual(dni_des=900., num_rays=int(1e6), nd=pm.nd, nh=pm.nh, weafile=weafile, method=1, Q_in_des=pm.Q_in_rcv, n_helios=None, zipfiles=False, gen_vtk=False, plot=False) - if (A_land==0): + if (A_land==0): tablefile=None - else: + else: A_helio=pm.H_helio*pm.W_helio output_matadata_motab(table=oelt, field_type=pm.field_type, aiming='single', n_helios=self.n_helios, A_helio=A_helio, eff_design=self.eff_des, H_rcv=pm.H_rcv, W_rcv=pm.W_rcv, H_tower=pm.H_tower, Q_in_rcv=pm.Q_in_rcv, A_land=A_land, savedir=tablefile) end=time.time() print('total time %.2f'%((end-start)/60.), 'min') - - - diff --git a/solsticepy/design_dish.py b/solsticepy/design_dish.py index d105e34..506946f 100644 --- a/solsticepy/design_dish.py +++ b/solsticepy/design_dish.py @@ -21,18 +21,16 @@ class Dish: the sun, the parabolic concentrator and the receiver. ''' - def __init__(self, casedir, nproc=None, verbose=False): + def __init__(self, casedir): ''' Arguements: casedir : str, the directory of the case ''' self.casedir=casedir - self.nproc=nproc - self.verbose=verbose if not os.path.exists(casedir): os.makedirs(casedir) - self.master=Master(casedir, nproc) + self.master=Master(casedir) def yaml(self, dish_radius, dish_foc, rho_refl, slope_error, rec_r, rec_x, rec_y, rec_z, rec_grid_r, rec_abs, dni=1000, sunshape=None, csr=0.01, half_angle_deg=0.2664, std_dev=0.2): ''' @@ -171,12 +169,12 @@ def get_opt_eff(self, dni_des, num_rays, zipfiles=False, gen_vtk=False, plot=Fal ''' Assume the optical efficiency of a dish system is constant at different sun positions - ''' + ''' azi_des=0 ele_des=90 sys.stderr.write("\n"+green('Optical efficiency: \n')) - efficiency_total=self.master.run(azi_des, ele_des, num_rays, self.rho_refl, dni_des, folder=self.casedir, gen_vtk=gen_vtk, printresult=True, verbose=self.verbose, system='dish') + efficiency_total=self.master.run(azi_des, ele_des, num_rays, self.rho_refl, dni_des, folder=self.casedir, gen_vtk=gen_vtk, printresult=True, system='dish') eff_des=efficiency_total.n return eff_des diff --git a/solsticepy/gen_yaml.py b/solsticepy/gen_yaml.py index 6a793ad..296ecaf 100755 --- a/solsticepy/gen_yaml.py +++ b/solsticepy/gen_yaml.py @@ -1,12 +1,11 @@ from __future__ import print_function import numpy as np -import matplotlib.pyplot as plt #for python 2: #from builtins import super from .data_spectral import SolarSpectrum, MirrorRhoSpectrum -from .cal_layout import multi_aperture_pos +from .design_cpc import * import sys class Sun: @@ -48,9 +47,9 @@ def yaml(self,spectrum = None): s += ", spectrum = %s" % (spectrum,) if self.sunshape is not None: if self.sunshape=='pillbox': - s += ", pillbox: {half_angle: %6.4f}" % (self.half_angle_deg,) + s += ", pillbox: {half_angle: %6.4f}" % (self.half_angle_deg,) elif self.sunshape=='buie': - s += ", buie: {csr: %6.4f}" % (self.csr,) + s += ", buie: {csr: %6.4f}" % (self.csr,) elif self.sunshape=='gaussian': s += ", gaussian: {std_dev: %6.4f}" % (self.std_dev,) s += "}" @@ -80,23 +79,23 @@ def gen_yaml(sun, hst_pos, hst_foc, hst_aims,hst_w, hst_h * `tower_h` (float): tower height (m) * `tower_r` (float): tower radius (a cylindrical shape tower) (m) 3. the receiver - * `receiver` (str): ``'flat'``, ``'cylinder'``, or ``'stl' or 'multi-aperture'`` (first of the 'receiver' parameters) + * `receiver` (str): ``'flat'``, ``'cylinder'``, or ``'stl' or 'multi_cavity'`` (first of the 'receiver' parameters) * `rec_abs` (float): receiver absorptivity * `rec_param` (numpy array or str): each element contains the geometrical parameter of the corresponding receiver. 4. others * `spectral` (bool): True - simulate the spectral dependent performance (first of the 'other' parameters) * `medium` (float): if the atmosphere is surrounded by non-participant medium, medium=0; otherwise it is the extinction coefficient in m-1 * `one_heliosat` (boolean): if `True`, implements ray tracing from just one heliostat. - + Returns: nothing (requested files are created and written) - Note that the parameters are in groups that relate to the `sun`, the `field` and the `receiver` then `others`. + Note that the parameters are in groups that relate to the `sun`, the `field` and the `receiver` then `others`. Note also the type for `rec_param` should be as follows. - * if ``receiver == 'flat'``: np.array([width, height, grid_w, grid_h,, x, y, z, tilt angle (deg))] - * if ``receiver == 'cylinder'``: np.array([radius, height, grid_circ, grid_h, x, y, z, tilt angle (deg)]) + * if ``receiver == 'flat'``: np.array([width, height, slices, x, y, z, tilt angle (deg))] + * if ``receiver == 'cylinder'``: np.array([radius, height, slices]) * if ``receiver == 'stl'``: the directory of the stl file - * if ``receiver == 'multi-aperture'``: np.array([width, height, grid_w, grid_h,, x, y, z, tilt angle (deg),num_aperture, gamma (deg) ]) + * if ``receiver == 'multi_cavity'``: the directory of the stl file """ # FIXME Parameters should be named according to what they are, eg # the parameter should be called 'csr', not 'sunsize', to avoid confusion. @@ -107,9 +106,9 @@ def gen_yaml(sun, hst_pos, hst_foc, hst_aims,hst_w, hst_h iyaml='' # the input yaml file - # + # ### Section (1) - # set the spectral data: + # set the spectral data: # solar radiative intensity, refractive indexes, extinction coefficients, reflectivities #------------------------------ if spectral: @@ -117,9 +116,9 @@ def gen_yaml(sun, hst_pos, hst_foc, hst_aims,hst_w, hst_h # CREATE the spectrum for the sun iyaml+='- spectrum: &solar_spectrum \n' for i in range(0,len(I_sun)-1): - iyaml+=' - {wavelength: %s, data: %s }\n' % (I_sun[i][0],I_sun[i][1]) + iyaml+=' - {wavelength: %s, data: %s }\n' % (I_sun[i][0],I_sun[i][1]) i = len(I_sun)-1 - iyaml+=' - {wavelength: %s, data: %s }\n' % (I_sun[i][0],I_sun[i][1]) + iyaml+=' - {wavelength: %s, data: %s }\n' % (I_sun[i][0],I_sun[i][1]) iyaml+='\n' # CREATE the spectrum for the reflectivity (mirror) @@ -131,14 +130,14 @@ def gen_yaml(sun, hst_pos, hst_foc, hst_aims,hst_w, hst_h mirror_ref.append([4,0.9]) iyaml+='- spectrum: &%s \n' % 'ref_mirror' for i in range(0,len(mirror_ref)-1): - iyaml+=' - {wavelength: %15.8e, data: %15.8e }\n' % (float(mirror_ref[i][0]),float(mirror_ref[i][1])) + iyaml+=' - {wavelength: %15.8e, data: %15.8e }\n' % (float(mirror_ref[i][0]),float(mirror_ref[i][1])) i = len(mirror_ref)-1 - iyaml+=' - {wavelength: %15.8e, data: %15.8e }\n' % (float(mirror_ref[i][0]),float(mirror_ref[i][1])) + iyaml+=' - {wavelength: %15.8e, data: %15.8e }\n' % (float(mirror_ref[i][0]),float(mirror_ref[i][1])) iyaml+='\n' - # + # ### Section (2) - # set the medium types: + # set the medium types: # air, glass, vacuum, etc. gathering spectral data #------------------------------ # @@ -150,15 +149,15 @@ def gen_yaml(sun, hst_pos, hst_foc, hst_aims,hst_w, hst_h spectrum = "*solar_spectrum" else: spectrum = None - + iyaml += "- sun: %s\n" % (sun.yaml(spectrum),) if medium>1e-99: - iyaml+='- atmosphere: {extinction: %s}\n'%medium + iyaml+='- atmosphere: {extinction: %s}\n'%medium iyaml+='\n' - - # + + # ### Section (3) # set the materials # (gathering media) @@ -170,7 +169,7 @@ def gen_yaml(sun, hst_pos, hst_foc, hst_aims,hst_w, hst_h r_b = 0. # and back reflectivity iyaml+='- material: &%s\n' % 'material_black' iyaml+=' front:\n' - iyaml+=' matte: {reflectivity: %6.4f }\n' % r_f + iyaml+=' matte: {reflectivity: %6.4f }\n' % r_f iyaml+=' back:\n' iyaml+=' matte: {reflectivity: %6.4f }\n' % r_b iyaml+='\n' @@ -181,12 +180,12 @@ def gen_yaml(sun, hst_pos, hst_foc, hst_aims,hst_w, hst_h iyaml+='- material: &%s\n' % 'material_mirror' iyaml+=' front:\n' if spectral: - iyaml+=' mirror: {reflectivity: *%s, slope_error: %15.8e }\n' % ('ref_mirror', slope_error ) + iyaml+=' mirror: {reflectivity: *%s, slope_error: %15.8e }\n' % ('ref_mirror', slope_error ) else: - iyaml+=' mirror: {reflectivity: %6.4f, slope_error: %15.8e }\n' % (r_f, slope_error) + iyaml+=' mirror: {reflectivity: %6.4f, slope_error: %15.8e }\n' % (r_f, slope_error) iyaml+=' back:\n' - iyaml+=' matte: {reflectivity: %6.4f }\n' % r_b + iyaml+=' matte: {reflectivity: %6.4f }\n' % r_b iyaml+='\n' # # CREATE a material for the target @@ -194,7 +193,7 @@ def gen_yaml(sun, hst_pos, hst_foc, hst_aims,hst_w, hst_h r_b = 1.-rec_abs # and back reflectivity iyaml+='- material: &%s\n' % 'material_target' iyaml+=' front:\n' - iyaml+=' matte: {reflectivity: %6.4f }\n' % r_f + iyaml+=' matte: {reflectivity: %6.4f }\n' % r_f iyaml+=' back:\n' iyaml+=' matte: {reflectivity: %6.4f }\n' % r_b iyaml+='\n' @@ -205,7 +204,7 @@ def gen_yaml(sun, hst_pos, hst_foc, hst_aims,hst_w, hst_h iyaml+='\n' - # + # ### Section (4) # set the geometries # (gathering shapes and materials) @@ -216,10 +215,10 @@ def gen_yaml(sun, hst_pos, hst_foc, hst_aims,hst_w, hst_h # (cylindrical shape) # slices = 10 # slices for the envelop circle - iyaml+='- geometry: &%s\n' % 'tower_g' - iyaml+=' - material: *%s\n' % 'material_black' - #iyaml+=' transform: { translation: %s, rotation: %s }\n' % ([0, 0, h_tow*0.5], [0, 90, 0]) - iyaml+=' cylinder: {height: %7.3f, radius: %7.3f, slices: %d }\n' % (tower_h, tower_r, slices) + iyaml+='- geometry: &%s\n' % 'tower_g' + iyaml+=' - material: *%s\n' % 'material_black' + #iyaml+=' transform: { translation: %s, rotation: %s }\n' % ([0, 0, h_tow*0.5], [0, 90, 0]) + iyaml+=' cylinder: {height: %7.3f, radius: %7.3f, slices: %d }\n' % (tower_h, tower_r, slices) iyaml+='\n' # # Receiver Geometry @@ -235,9 +234,13 @@ def gen_yaml(sun, hst_pos, hst_foc, hst_aims,hst_w, hst_h elif receiver=='stl': rec_entt, rcv=STL_receiver(rec_param, hemisphere) - elif receiver=='multi-aperture': - geom, rec_entt, rcv =multi_aperture_receiver(rec_param, hemisphere) - iyaml+=geom + elif receiver=='multi_cavity': + geom, rec_entt, rcv =multi_cavity_receiver(rec_param, hemisphere) + + elif receiver=='beam_down': + sys.stderr.write("You chose a beam-down system\n") + cpc=CPC() + rec_entt, rcv = cpc.beamdowncomponents(rec_param) # # Heliostats Geometry # @@ -245,7 +248,7 @@ def gen_yaml(sun, hst_pos, hst_foc, hst_aims,hst_w, hst_h hst_x=np.r_[hst_pos[0]] hst_y=np.r_[hst_pos[1]] hst_z=np.r_[hst_pos[2]] - aim_x=np.r_[hst_aims[0]] + aim_x=np.r_[hst_aims[0]] aim_y=np.r_[hst_aims[1]] aim_z=np.r_[hst_aims[2]] num_hst=1 @@ -263,50 +266,50 @@ def gen_yaml(sun, hst_pos, hst_foc, hst_aims,hst_w, hst_h # CREATE a reflective facet (mirror) for i in range(0,num_hst): name_hst_g = 'hst_g_'+str(i) - iyaml+='- geometry: &%s\n' % name_hst_g - iyaml+=' - material: *%s\n' % 'material_mirror' + iyaml+='- geometry: &%s\n' % name_hst_g + iyaml+=' - material: *%s\n' % 'material_mirror' #iyaml+=' transform: { translation: %s, rotation: %s }\n' % ([hst_x[i], hst_y[i], hst_z[i]], [0, 0, 0]) ) iyaml+=' parabol: \n' iyaml+=' focal: %s\n' % hst_foc[i] - iyaml+=' clip: \n' + iyaml+=' clip: \n' iyaml+=' - operation: AND \n' iyaml+=' vertices: %s\n' % pts_hst - iyaml+=' slices: %d\n' % slices + iyaml+=' slices: %d\n' % slices # CREATE the pylon "pylon_g" geometry cylindrical shape h_pyl = 0.001 # pylon height r_pyl = 0.2 # pylon radius slices = 4 # slices for the envelop circle - iyaml+='- geometry: &%s\n' % 'pylon_g' - iyaml+=' - material: *%s\n' % 'material_black' - iyaml+=' transform: { translation: %s, rotation: %s }\n' % ([0, 0, -h_pyl*3], [0, 90, 0]) - iyaml+=' cylinder: {height: %7.3f, radius: %7.3f, slices: %d }\n' % (h_pyl,r_pyl,slices) - # + iyaml+='- geometry: &%s\n' % 'pylon_g' + iyaml+=' - material: *%s\n' % 'material_black' + iyaml+=' transform: { translation: %s, rotation: %s }\n' % ([0, 0, -h_pyl*3], [0, 90, 0]) + iyaml+=' cylinder: {height: %7.3f, radius: %7.3f, slices: %d }\n' % (h_pyl,r_pyl,slices) + # - # + # ### Section (5) # set the templates # (programming objects gathering geometries or pivot and geometries) #------------------------------ # CREATE the heliostat templates - for i in range(0,num_hst): + for i in range(0,num_hst): name_hst_t = 'hst_t_'+str(i) - iyaml+='- template: &%s\n' % name_hst_t + iyaml+='- template: &%s\n' % name_hst_t name_hst_n = 'hst_'+ str(i) - iyaml+=' name: %s\n' % name_hst_n - iyaml+=' primary: 0\n' + iyaml+=' name: %s\n' % name_hst_n + iyaml+=' primary: 0\n' iyaml+=' geometry: *pylon_g\n' - iyaml+=' children: \n' + iyaml+=' children: \n' iyaml+=' - name: pivot\n' - iyaml+=' zx_pivot: {target: {position: %s}} \n' % ([aim_x[i],aim_y[i],aim_z[i]]) + iyaml+=' zx_pivot: {target: {position: %s}} \n' % ([aim_x[i],aim_y[i],aim_z[i]]) iyaml+=' children: \n' iyaml+=' - name: reflect_surface\n' iyaml+=' primary: 1\n' - iyaml+=' transform: {rotation: [-90,0,0]} \n' + iyaml+=' transform: {rotation: [-90,0,0]} \n' name_hst_g = 'hst_g_'+str(i) - iyaml+=' geometry: *%s\n' % name_hst_g + iyaml+=' geometry: *%s\n' % name_hst_g - # + # ### Section (6) # set the entities # (gather templates to be created and active in the scene) @@ -318,9 +321,9 @@ def gen_yaml(sun, hst_pos, hst_foc, hst_aims,hst_w, hst_h # tower entities iyaml+='\n- entity:\n' iyaml+=' name: tower_e\n' - iyaml+=' primary: 0\n' - iyaml+=' transform: { translation: %s, rotation: %s }\n' % ([0, -tower_r, tower_h*0.5], [0, 0, 0]) - iyaml+=' geometry: *%s\n' % 'tower_g' + iyaml+=' primary: 0\n' + iyaml+=' transform: { translation: %s, rotation: %s }\n' % ([0, -tower_r, tower_h*0.5], [0, 0, 0]) + iyaml+=' geometry: *%s\n' % 'tower_g' # # heliostat entities from the template for i in range(0,num_hst): @@ -328,23 +331,23 @@ def gen_yaml(sun, hst_pos, hst_foc, hst_aims,hst_w, hst_h name_hst_t = 'hst_t_'+str(i) iyaml+='\n- entity:\n' iyaml+=' name: %s\n' % name_e - iyaml+=' transform: { translation: %s, rotation: %s }\n' % ([hst_x[i], hst_y[i], hst_z[i]], [0, 0, 0]) - iyaml+=' children: [ *%s ]\n' % name_hst_t + iyaml+=' transform: { translation: %s, rotation: %s }\n' % ([hst_x[i], hst_y[i], hst_z[i]], [0, 0, 0]) + iyaml+=' children: [ *%s ]\n' % name_hst_t with open(outfile_yaml,'w') as f: f.write(iyaml) with open(outfile_recv,'w') as f: - f.write(rcv) + f.write(rcv) def flat_receiver(rec_param, hemisphere='North'): """ hemisphere : 'North' or 'South' hemisphere of the earth where the field located if North: the field is in the positive y direction - if South: the field is in the negtive y direction + if South: the field is in the negative y direction this will influence: - (1) the setting of the receiver tilt angle, + (1) the setting of the receiver tilt angle, if the front surface always facing to the field is desirable (2) the position of the virtual target """ @@ -366,10 +369,10 @@ def flat_receiver(rec_param, hemisphere='North'): geom+='- geometry: &%s\n' % 'target_g' geom+=' - material: *%s\n' % 'material_target' geom+=' plane: \n' - geom+=' clip: \n' + geom+=' clip: \n' geom+=' - operation: AND \n' geom+=' vertices: %s\n' % pts - geom+=' slices: %d\n' % slices + geom+=' slices: %d\n' % slices geom+='\n' # CREATE a receiver entity from "target_g" geometry (primary = 0) @@ -378,9 +381,9 @@ def flat_receiver(rec_param, hemisphere='North'): entt+=' name: target_e\n' entt+=' primary: 0\n' if hemisphere=='North': - entt+=' transform: { translation: %s, rotation: %s }\n' % ([x, y, z], [-90.-tilt, 0, 0]) + entt+=' transform: { translation: %s, rotation: %s }\n' % ([x, y, z], [-90.-tilt, 0, 0]) else: - entt+=' transform: { translation: %s, rotation: %s }\n' % ([x, y, z], [90.+tilt, 0, 0]) + entt+=' transform: { translation: %s, rotation: %s }\n' % ([x, y, z], [90.+tilt, 0, 0]) entt+=' geometry: *%s\n' % 'target_g' # CREATE a virtual target entity from "target_g" geometry (primary = 0) @@ -393,16 +396,16 @@ def flat_receiver(rec_param, hemisphere='North'): entt+=' transform: { translation: %s, rotation: %s }\n' % ([x, y-5., z], [-90.-tilt, 0, 0]) else: entt+=' transform: { translation: %s, rotation: %s }\n' % ([x, y+5., z], [90.+tilt, 0, 0]) - entt+=' geometry: \n' - entt+=' - material: *%s\n' % 'material_virtual' + entt+=' geometry: \n' + entt+=' - material: *%s\n' % 'material_virtual' entt+=' plane: \n' - entt+=' clip: \n' + entt+=' clip: \n' entt+=' - operation: AND \n' entt+=' vertices: %s\n' % pts - entt+=' slices: %d\n' % slices + entt+=' slices: %d\n' % slices rcv='' - rcv+='- name: target_e \n' + rcv+='- name: target_e \n' rcv+=' side: %s \n' % 'FRONT_AND_BACK' rcv+=' per_primitive: %s \n' % 'INCOMING_AND_ABSORBED' rcv+='- name: virtual_target_e \n' @@ -419,7 +422,7 @@ def cylindrical_receiver(rec_param, hemisphere='North'): if North: the field is in the positive y direction if South: the field is in the negtive y direction this will influence: - (1) the setting of the receiver tilt angle, + (1) the setting of the receiver tilt angle, if the front surface always facing to the field is desirable (2) the position of the virtual target ''' @@ -435,10 +438,10 @@ def cylindrical_receiver(rec_param, hemisphere='North'): geom+='- geometry: &%s\n' % 'target_g' geom+=' - material: *%s\n' % 'material_target' geom+=' cylinder: \n' - geom+=' height: %s\n'%rec_h - geom+=' radius: %s\n'%rec_r - geom+=' slices: %d\n' % slices - geom+=' stacks: %d\n' % stacks + geom+=' height: %s\n'%rec_h + geom+=' radius: %s\n'%rec_r + geom+=' slices: %d\n' % slices + geom+=' stacks: %d\n' % stacks geom+='\n' # CREATE a receiver entity from "target_g" geometry (primary = 0) @@ -447,7 +450,7 @@ def cylindrical_receiver(rec_param, hemisphere='North'): entt+=' name: target_e\n' entt+=' primary: 0\n' - entt+=' transform: { translation: %s, rotation: %s }\n' % ([x, y, z], [0., 0., 0.]) + entt+=' transform: { translation: %s, rotation: %s }\n' % ([x, y, z], [0., 0., 0.]) entt+=' geometry: *%s\n' % 'target_g' @@ -461,16 +464,16 @@ def cylindrical_receiver(rec_param, hemisphere='North'): entt+=' transform: { translation: %s, rotation: %s }\n' % ([x, y, z+rec_h/2.+1], [-180., 0, 0]) - entt+=' geometry: \n' - entt+=' - material: *%s\n' % 'material_virtual' + entt+=' geometry: \n' + entt+=' - material: *%s\n' % 'material_virtual' entt+=' plane: \n' - entt+=' clip: \n' + entt+=' clip: \n' entt+=' - operation: AND \n' entt+=' vertices: %s\n' % pts - entt+=' slices: %d\n' % slices + entt+=' slices: %d\n' % slices rcv='' - rcv+='- name: target_e \n' + rcv+='- name: target_e \n' rcv+=' side: %s \n' % 'FRONT_AND_BACK' rcv+=' per_primitive: %s \n' % 'INCOMING_AND_ABSORBED' rcv+='- name: virtual_target_e \n' @@ -479,7 +482,7 @@ def cylindrical_receiver(rec_param, hemisphere='North'): return geom, entt, rcv - + def STL_receiver(rec_param, hemisphere='North'): ''' @@ -487,7 +490,7 @@ def STL_receiver(rec_param, hemisphere='North'): if North: the field is in the positive y direction if South: the field is in the negtive y direction this will influence: - (1) the setting of the receiver tilt angle, + (1) the setting of the receiver tilt angle, if the front surface always facing to the field is desirable (2) the position of the virtual target ''' @@ -500,17 +503,17 @@ def STL_receiver(rec_param, hemisphere='North'): z=rec_param[5].astype(float) tilt=rec_param[6].astype(float) # need to figure out the initial mesh orientation - # CREATE a receiver entity from a STL file + # CREATE a receiver entity from a STL file entt='' entt+='\n- entity:\n' entt+=' name: STL_receiver_e\n' entt+=' primary: 0\n' if hemisphere=='North': - entt+=' transform: { translation: %s, rotation: %s }\n' % ([x, y, z], [-90.-tilt, 0, 0]) + entt+=' transform: { translation: %s, rotation: %s }\n' % ([x, y, z], [-90.-tilt, 0, 0]) else: # if it is the mesh model of the bladed receiver at CSIRO - entt+=' transform: { translation: %s, rotation: %s }\n' % ([x, y, z], [180.+tilt, 0, 0]) + entt+=' transform: { translation: %s, rotation: %s }\n' % ([x, y, z], [180.+tilt, 0, 0]) entt+=' geometry:\n' entt+=' - material: *material_target\n' entt+=' transform: {translation: [0, 0, 0], rotation: [0, 0, 0]}\n' @@ -527,16 +530,16 @@ def STL_receiver(rec_param, hemisphere='North'): entt+=' transform: { translation: %s, rotation: %s }\n' % ([x, y-5., z], [-90.-tilt, 0, 0]) else: entt+=' transform: { translation: %s, rotation: %s }\n' % ([x, y+5., z], [90.+tilt, 0, 0]) - entt+=' geometry: \n' - entt+=' - material: *%s\n' % 'material_virtual' + entt+=' geometry: \n' + entt+=' - material: *%s\n' % 'material_virtual' entt+=' plane: \n' - entt+=' clip: \n' + entt+=' clip: \n' entt+=' - operation: AND \n' entt+=' vertices: %s\n' % pts - entt+=' slices: %d\n' % slices + entt+=' slices: %d\n' % slices rcv='' - rcv+='- name: STL_receiver_e \n' + rcv+='- name: STL_receiver_e \n' rcv+=' side: %s \n' % 'FRONT_AND_BACK' rcv+=' per_primitive: %s \n' % 'INCOMING_AND_ABSORBED' rcv+='- name: virtual_target_e \n' @@ -545,98 +548,82 @@ def STL_receiver(rec_param, hemisphere='North'): return entt, rcv -def multi_aperture_receiver(rec_param, hemisphere='North', plot=False): +def multi_cavity_receiver(rec_param, hemisphere='North'): """ hemisphere : 'North' or 'South' hemisphere of the earth where the field located if North: the field is in the positive y direction if South: the field is in the negtive y direction this will influence: - (1) the setting of the receiver tilt angle, + (1) the setting of the receiver tilt angle, if the front surface always facing to the field is desirable (2) the position of the virtual target """ - # rec_w and rec_h is the size of one aperture - # rec_grid_w and rec_gird_h is the number of elements of one aperture - # rec_z is a list of the elevation height of the center of apertures - # rec_tilt is the tilt angle of each aperture (the default is facing to the horizon) - # num_aperture is the number of apertures - # gamma is the angular range of the multi-aperture configration - + # x, y, z is the front center of the multi-cavity receiver + # n_c is the number of cavities + # rec_w, rec_h is the size of one cavity + # alpha is the angle between each two cavities + # num_cavity is the number of cavities rec_w=rec_param[0] rec_h=rec_param[1] - rec_grid_w=rec_param[2] - rec_grid_h=rec_param[3] - - rec_z=rec_param[4] - rec_tilt=rec_param[5] + slices=rec_param[2] + x=rec_param[3] + y=rec_param[4] + z=rec_param[5] + tilt=rec_param[6] # receiver tilt angle: # 0 is vertical # the standby posiion of a plane in solstice is normal points to the +z axis # rotation anagle, positive is anti-clockwise - num_aperture=int(rec_param[6]) - gamma=rec_param[7] # angular range of the multi-aperture configration (deg) - geom='' - entt='' - vir_z=0. - for i in range(num_aperture): - - pts=[ [-rec_w[i]*0.5, -rec_h[i]*0.5], [-rec_w[i]*0.5, rec_h[i]*0.5], [rec_w[i]*0.5, rec_h[i]*0.5], [rec_w[i]*0.5,-rec_h[i]*0.5] ] - - geom+='- geometry: &%s\n' % 'target_g_%.0f\n'%(i) - geom+=' - material: *%s\n' % 'material_target' - geom+=' plane: \n' - geom+=' clip: \n' - geom+=' - operation: AND \n' - geom+=' vertices: %s\n' % pts - geom+=' slices: %d\n' % rec_grid_w - geom+='\n' - - ang_pos, xc, yc=multi_aperture_pos(rec_w, gamma, num_aperture, i) - - zc=rec_z[i] - vir_z+=zc - - # CREATE a receiver entity from "target_g" geometry (primary = 0) + pts=[ [-rec_w*0.5, -rec_h*0.5], [-rec_w*0.5, rec_h*0.5], [rec_w*0.5, rec_h*0.5], [rec_w*0.5,-rec_h*0.5] ] - entt+='\n- entity:\n' - entt+=' name: target_e_%.0f\n'%(i) - entt+=' primary: 0\n' - if hemisphere=='North': - entt+=' transform: { translation: %s, rotation: %s }\n' % ([xc, yc, zc], [-90.-rec_tilt, 90.-ang_pos,0]) - else: - entt+=' transform: { translation: %s, rotation: %s }\n' % ([-xc, -yc, zc], [90.+rec_tilt, 90.-ang_pos,0]) - entt+=' geometry: *%s\n' % 'target_g_%.0f\n'%(i) + geom+='- geometry: &%s\n' % 'target_g' + geom+=' - material: *%s\n' % 'material_target' + geom+=' plane: \n' + geom+=' clip: \n' + geom+=' - operation: AND \n' + geom+=' vertices: %s\n' % pts + geom+=' slices: %d\n' % slices + geom+='\n' - vir_z/=float(num_aperture) + # CREATE a receiver entity from "target_g" geometry (primary = 0) + entt='' + entt+='\n- entity:\n' + entt+=' name: target_e\n' + entt+=' primary: 0\n' + if hemisphere=='North': + entt+=' transform: { translation: %s, rotation: %s }\n' % ([x, y, z], [-90.-tilt, 0, 0]) + else: + entt+=' transform: { translation: %s, rotation: %s }\n' % ([x, y, z], [90.+tilt, 0, 0]) + entt+=' geometry: *%s\n' % 'target_g' # CREATE a virtual target entity from "target_g" geometry (primary = 0) - slices = 16 - radius=vir_z*0.5 + pts = [ [-rec_w*10., -rec_h*10.], [-rec_w*10., rec_h*10.], [rec_w*10., rec_h*10.], [rec_w*10.,-rec_h*10.] ] + slices = 4 entt+='\n- entity:\n' entt+=' name: virtual_target_e\n' entt+=' primary: 0\n' - entt+=' transform: { translation: %s}\n' % ([0., 0., vir_z]) - entt+=' geometry: \n' - entt+=' - material: *%s\n' % 'material_virtual' - entt+=' sphere: \n' - entt+=' radius: %s\n' % radius - entt+=' slices: %d\n' % slices + if hemisphere=='North': + entt+=' transform: { translation: %s, rotation: %s }\n' % ([x, y-5., z], [-90.-tilt, 0, 0]) + else: + entt+=' transform: { translation: %s, rotation: %s }\n' % ([x, y+5., z], [90.+tilt, 0, 0]) + entt+=' geometry: \n' + entt+=' - material: *%s\n' % 'material_virtual' + entt+=' plane: \n' + entt+=' clip: \n' + entt+=' - operation: AND \n' + entt+=' vertices: %s\n' % pts + entt+=' slices: %d\n' % slices rcv='' - for i in range(num_aperture): - rcv+='- name: target_e_%.0f \n'%(i) - rcv+=' side: %s \n' % 'FRONT_AND_BACK' - rcv+=' per_primitive: %s \n' % 'INCOMING_AND_ABSORBED' + rcv+='- name: target_e \n' + rcv+=' side: %s \n' % 'FRONT_AND_BACK' + rcv+=' per_primitive: %s \n' % 'INCOMING_AND_ABSORBED' rcv+='- name: virtual_target_e \n' rcv+=' side: %s \n' % 'FRONT' rcv+=' per_primitive: %s \n' % 'INCOMING' return geom, entt, rcv - - - #------------------------------ - diff --git a/solsticepy/input.py b/solsticepy/input.py index 067fa0c..0c172a6 100644 --- a/solsticepy/input.py +++ b/solsticepy/input.py @@ -24,11 +24,11 @@ def __init__(self): def Sun(self): ''' (1) lat : float, latitude of the plant lcation - (2) dni_des : float, dni at design point (W/m2) - (3) sunshape : str, 'pillbox' or 'Buie' - (4) sunsize : float, + (2) dni_des : float, dni at design point (W/m2) + (3) sunshape : str, 'pillbox' or 'Buie' + (4) sunsize : float, -- for 'pillbox': half angle of the pillbox sunshape (degree) - -- for 'Buie': circumsolar ratio (CSR) + -- for 'Buie': circumsolar ratio (CSR) (5) extinction: float, around 1e-6, the extinction coefficient if the atmosphere is surrounded by participant medium ''' @@ -45,35 +45,34 @@ def Sun(self): def Heliostat(self): ''' (1) field_type : str, - -- 'polar', 'surround', 'multi-aperture', 'polar-half' or 'surround-half for desiging a new field - the 'half' option is for simulation of a symmetric field + -- 'polar', 'polar-half', 'surround' or 'surround-half' for desiging a new field -- the directory of the layout file the layout file is a 'csv' file, (n+2, 7) - - n is the total number of heliostats + - n is the total number of heliostats - the 1st row is the number of each column - - the 2nd row is the unit + - the 2nd row is the unit - the first three columns are X, Y, Z coordinates of each heliostat - the fourth column is the focal length - the last three columns are the aiming points (2) Q_in_rcv : float, required heat of the receiver (W) - (3) W_helio : float, width of a heliostat (m) + (3) W_helio : float, width of a heliostat (m) (4) H_helio : float, height of a heliostat (m) (5) slope_error: float, slope error(radians) - (6) rho_helio : float, reflectivity of heliostat + (6) rho_helio : float, reflectivity of heliostat (7) H_tower : float, tower height (m) (8) R_tower : float, radius of tower (m) (9) concret_tower: bool, True-a solid concrete tower, or False-a truss tower (10)single_filed : bool, True-one tower one field, or False-multi-tower - (11)R1 : float, layout parameter, the distance of the first row of heliostat + (11)R1 : float, layout parameter, the distance of the first row of heliostat (12)dsep : float, layout parameter, the separation distance of heliostats (m) - (13)fb : float, in (0-1), a factor to expand the field to reduce block - ---(*) n_helios: int, number of heliostats for the designed field + (13)fb : float, in (0-1), a factor to expand the field to reduce block + ---(*) n_helios: int, number of heliostats for the designed field ---(*) Z_helio : float, the installation height of the heliostat ''' self.field_type='polar' if self.method==1: - self.Q_in_rcv=10e6 # required heat of the receiver + self.Q_in_rcv=10e6 # required heat of the receiver else: self.n_helios=1000 self.W_helio=10. @@ -83,43 +82,36 @@ def Heliostat(self): self.H_tower=100. self.R_tower=0.001 # shading effect of tower is neglected at the moment self.concret_tower=False - self.single_field=True + self.single_field=True - self.R1=90. + self.R1=0.5*self.H_tower self.dsep=0. self.fb=0.7 def Receiver(self): ''' - (1) rcv_type : str, 'flat', 'cylinder', 'particle', 'multi-aperture' or 'stl', type of the receiver - (2) num_aperture: int, number of apertures if it is a multi-aperture receiver - (3) alpha : float, the angular space between two adjacent apertures (except the aperture faces to the South) (deg) - (4) H_rcv : float, height of a flat receiver or radius of a cylindrical receiver (m) - (5) W_rcv : float, width of the receiver (m) - (6) tilt_rcv : float, tilt angle of the receiver (deg), 0 is where the receiver face to the horizontal - (7) alpha_rcv : float, receiver surface absorptivity (0-1), set as 1 if rcv_type='particle' - (8) n_H_rcv : int, number of discretisation of the receiver in the height direction - (9) n_W_rcv : int, number of discretisation of the receiver in the width direction (for a flat receiver) - or in the circular direction (for a cylindrical receiver) - (10) X_rcv : float, x location of the receiver (m) - (11) Y_rcv : float, y location of the receiver (m) - (12) Z_rcv : float, z location of the receiver (m) - (13) num_aperture: int, number of apertures for a multi-aperture configuration - (14) gamma : float, the angular range of the multi-aperture configuration (deg) + (1) rcv_type : str, 'flat', 'cylinder', 'particle' or 'stl', type of the receiver + (2) H_rcv : float, height of a flat receiver or radius of a cylindrical receiver (m) + (3) W_rcv : float, width of the receiver (m) + (4) tilt_rcv : float, tilt angle of the receiver (deg), 0 is where the receiver face to the horizontal + (5) alpha_rcv : float, receiver surface absorptivity (0-1), set as 1 if rcv_type='particle' + (6) n_H_rcv : int, number of discretisation of the receiver in the height direction + (7) n_W_rcv : int, number of discretisation of the receiver in the width direction (for a flat receiver) + or in the circular direction (for a cylindrical receiver) + (8) X_rcv : float, x location of the receiver (m) + (9) Y_rcv : float, y location of the receiver (m) + ---(*) z_rcv : float, z location of the receiver (m) ''' self.rcv_type='flat' self.H_rcv=10. self.W_rcv=10. self.tilt_rcv=0. # receiver tilt angle - self.alpha_rcv=1. + self.alpha_rcv=1. self.n_H_rcv=10 self.n_W_rcv=10 self.X_rcv=0. # receiver location self.Y_rcv=0. - self.Z_rcv=120. - self.num_aperture=1 - self.gamma=0. def simulation(self): @@ -127,9 +119,9 @@ def simulation(self): n_row_oelt: int, number of rows of the lookup table (i.e. simulated days per year) n_col_oelt: int, number of columns of the lookup table (i.e. simulated number of hours per day) n_rays : int, number of rays for the simulation - n_procs : int, number of processors for the mcrt simulation + n_procs : int, number of processors for the mcrt simulation casedir : str, the directory for saving the result files - ''' + ''' self.n_row_oelt=5 self.n_col_oelt=5 self.n_rays=int(5e6) @@ -137,68 +129,61 @@ def simulation(self): self.casedir='.' self.method=1 # 1 - design the field based on the Q_in_rcv # 2 - design the field based on the n_helios - + def dependent_par(self): ''' ''' - self.Z_helio=self.H_helio*0.7 + self.Z_helio=self.H_helio*0.7 + self.Z_rcv=self.H_tower if self.lat>=0: self.hemisphere='North' elif self.lat<0: self.hemisphere='South' # estimate a rough number of large field - eta_field=0.4 # assumed field effieicy at design point + eta_field=0.4 # assumed field efficiency at design point if self.method==1: - if isinstance(self.Q_in_rcv, float): - self.n_helios=self.Q_in_rcv/self.W_helio/self.H_helio/self.dni_des/eta_field - else: - self.n_helios=sum(self.Q_in_rcv)/self.W_helio/self.H_helio/self.dni_des/eta_field - - if self.field_type!='surround': - self.n_helios*=2. + self.n_helios=self.Q_in_rcv/self.W_helio/self.H_helio/self.dni_des/eta_field + if self.field_type=='polar': + self.n_helios*=2. def saveparam(self, savedir): if not os.path.exists(savedir): os.makedirs(savedir) - + param=np.array([ - ['method', self.method, '-'], - ['field', self.field_type, '-'], - ['Q_in_rcv', self.Q_in_rcv, 'W'], - ['n_helios(pre_des if method ==1)', self.n_helios, '-'], - ['W_helio', self.W_helio, 'm'], + ['method', self.method, '-'], + ['field', self.field_type, '-'], + ['Q_in_rcv', self.Q_in_rcv, 'W'], + ['n_helios(pre_des if method ==1)', self.n_helios, '-'], + ['W_helio', self.W_helio, 'm'], ['H_helio', self.H_helio, 'm'], - ['Z_helio', self.Z_helio, 'm'], - ['rho_helio', self.rho_helio, '-'], + ['Z_helio', self.Z_helio, 'm'], + ['rho_helio', self.rho_helio, '-'], ['slope_error', self.slope_error, 'rad'], - ['H_tower', self.H_tower, 'm'], - ['R_tower', self.R_tower, 'm'], - ['concret_tower', self.concret_tower, '-'], - ['single_field', self.single_field, '-'], - ['fb factor', self.fb, '-'], - ['R1', self.R1, 'm'], - ['dsep', self.dsep, 'm'], - ['rcv_type', self.rcv_type, '-'], - ['num_aperture', self.num_aperture, '-'], - ['aperture angular range',self.gamma , 'deg'], + ['H_tower', self.H_tower, 'm'], + ['R_tower', self.R_tower, 'm'], + ['concret_tower', self.concret_tower, '-'], + ['single_field', self.single_field, '-'], + ['fb factor', self.fb, '-'], + ['R1', self.R1, 'm'], + ['dsep', self.dsep, 'm'], + ['rcv_type', self.rcv_type, '-'], ['H_rcv', self.H_rcv, 'm'], - ['W_rcv', self.W_rcv, 'm'], - ['tilt_rcv', self.tilt_rcv, 'deg'], - ['alpha_rcv', self.alpha_rcv, '-'], - ['n_H_rcv', self.n_H_rcv, '-'], - ['n_W_rcv', self.n_W_rcv, '-'], - ['X_rcv', self.X_rcv, 'm'], - ['Y_rcv', self.Y_rcv, 'm'], - ['Z_rcv', self.Z_rcv, 'm'], - ['n_row_oelt', self.n_row_oelt, '-'] , + ['W_rcv', self.W_rcv, 'm'], + ['tilt_rcv', self.tilt_rcv, 'deg'], + ['alpha_rcv', self.alpha_rcv, '-'], + ['n_H_rcv', self.n_H_rcv, '-'], + ['n_W_rcv', self.n_W_rcv, '-'], + ['X_rcv', self.X_rcv, 'm'], + ['Y_rcv', self.Y_rcv, 'm'], + ['Z_rcv', self.Z_rcv, 'm'], + ['n_row_oelt', self.n_row_oelt, '-'] , ['n_col_oelt', self.n_col_oelt, '-'] , ['n_rays', self.n_rays, '-'] , - ['n_procs', self.n_procs, '-'] + ['n_procs', self.n_procs, '-'] ]) np.savetxt(savedir+'/simulated_parameters.csv', param, delimiter=',', fmt='%s') - - diff --git a/solsticepy/master.py b/solsticepy/master.py index e3b2cf4..de95851 100755 --- a/solsticepy/master.py +++ b/solsticepy/master.py @@ -20,7 +20,7 @@ def SPROG(name): def run_prog(name,args,output_file=None,verbose=True): prog = SPROG(name) args1 = [str(a) for a in args] - if verbose: + if verbose: sys.stderr.write("Running '%s' with args: %s\n" % (name," ".join(args1))) if output_file is not None: # any error will cause an exception (and we capture the output to a file) @@ -33,18 +33,13 @@ def run_prog(name,args,output_file=None,verbose=True): class Master: - def __init__(self, casedir='.', nproc=None): + def __init__(self, casedir='.'): """Set up the Solstice simulation, i.e. establishing the case folder, calling the Solstice program and post-processing the results ``Argument`` * casedir (str): the case directory - * nproc (int): number of processors, e.g. nproc=1 will run in serial mode, - nproc=4 will run with 4 processors in parallel - nproc=None will run with any number of processors that are available """ self.casedir=os.path.abspath(casedir) - self.nproc=nproc - if not os.path.exists(self.casedir): os.makedirs(self.casedir) assert os.path.isdir(casedir) @@ -62,24 +57,23 @@ def in_case(self, folder, fn): * a joining directory of the file in the case directory """ - if not os.path.exists(folder): os.makedirs(folder) assert os.path.isdir(folder) return os.path.join(folder,fn) - def run(self, azimuth, elevation, num_rays, rho_mirror, dni, folder, gen_vtk=False, printresult=False, verbose=False, system='crs'): + def run(self, azimuth, elevation, num_rays, rho_mirror, dni, folder, gen_vtk=True, printresult=True, system='crs'): - """Run an optical simulation (one sun position) using Solstice + """Run an optical simulation (one sun position) using Solstice - * `azimuth` (float): the azimuth angle of the ray-tracing simulation in Solstice, counted from East towards to North + * `azimuth` (float): the azimuth angle of the ray-tracing simulation in Solstice, counted from East towards to North * `elevation` (float): the elevation angle of the ray-tracing simulation in Solstice - * `num_rays` (int): number of rays to be cast in the ray-tracing simulation - * `rho_mirror`(float): reflectivity of mirrors, required for results post-processing - * `dni` (float): the direct normal irradiance (W/m2), required to obtain performance of individual heliostat + * `num_rays` (int): number of rays to be cast in the ray-tracing simulation + * `rho_mirror` (float): reflectivity of mirrors, required for results post-processing + * `dni` (float): the direct normal irradiance (W/m2), required to obtain performance of individual heliostat * `gen_vtk` (boolean): if True, generate .vtk files for rendering in Paraview - * `system` (str): 'crs' for a central receiver system, or 'dish' for a parabolic dish system + * `system` (str): 'crs' for a central receiver system, or 'dish' for a parabolic dish system Returns: no return value (results files are created and written) """ @@ -88,13 +82,10 @@ def run(self, azimuth, elevation, num_rays, rho_mirror, dni, folder, gen_vtk=Fal RECV_IN = self.in_case(self.casedir, 'input-rcv.yaml') # main raytrace - if self.nproc==None: - run_prog("solstice",['-D%s,%s'%(azimuth,elevation),'-v','-n',num_rays,'-R',RECV_IN,'-fo',self.in_case(folder, 'simul'),YAML_IN]) - else: - run_prog("solstice",['-D%s,%s'%(azimuth,elevation),'-v', '-t', self.nproc, '-n',num_rays,'-R',RECV_IN,'-fo',self.in_case(folder, 'simul'),YAML_IN]) + run_prog("solstice",['-D%s,%s'%(azimuth,elevation),'-v','-n',num_rays,'-R',RECV_IN,'-fo',self.in_case(folder, 'simul'),YAML_IN]) folder=os.path.abspath(folder) - if gen_vtk and verbose: + if gen_vtk: dirn = os.getcwd() try: os.chdir(folder) @@ -120,26 +111,26 @@ def run(self, azimuth, elevation, num_rays, rho_mirror, dni, folder, gen_vtk=Fal os.chdir(dirn) if system=='dish': - eta=process_raw_results_dish(self.in_case(folder, 'simul'), folder, rho_mirror, dni, verbose=verbose) + eta=process_raw_results_dish(self.in_case(folder, 'simul'), folder, rho_mirror, dni) if printresult: sys.stderr.write('\n' + yellow("Total efficiency: {:f}\n".format(eta))) sys.stderr.write(green("Completed successfully.\n")) return eta else: - if system=='multi-aperture': - eta, performance_hst=process_raw_results_multi_aperture(self.in_case(folder, 'simul'), folder,rho_mirror, dni, verbose=verbose) + if system=='beamdown': + eta, performance_hst=process_raw_results(self.in_case(folder, 'simul'), folder, rho_mirror, dni, verbose=True, num_virt=2) else: - eta, performance_hst=process_raw_results(self.in_case(folder, 'simul'), folder,rho_mirror, dni, verbose=verbose) + eta, performance_hst=process_raw_results(self.in_case(folder, 'simul'), folder, rho_mirror, dni) if printresult: sys.stderr.write('\n' + yellow("Total efficiency: {:f}\n".format(eta))) sys.stderr.write(green("Completed successfully.\n")) return eta, performance_hst - def run_annual(self, nd, nh, latitude, num_rays, num_hst,rho_mirror,dni, gen_vtk=False,verbose=False): + def run_annual(self, nd, nh, latitude, num_rays, num_hst, rho_mirror, dni, gen_vtk=False, system='crs'): - """Run a list of optical simulations to obtain annual performance (lookup table) using Solstice + """Run a list of optical simulations to obtain annual performance (lookup table) using Solstice ``Arguments`` @@ -147,20 +138,15 @@ def run_annual(self, nd, nh, latitude, num_rays, num_hst,rho_mirror,dni, gen_vtk * nh (int): number of columns in the lookup table (discretisation of the solar hour angle) * latitude (float): the latitude angle of the plan location (deg) * num_rays (int): number of rays to be cast in the ray-tracing simulation - * num_hst (int): number of heliostats - * nproc (int): number of processors, e.g. nproc=1 will run in serial mode, - nproc=4 will run with 4 processors in parallel - nproc=None will run with any number of processors that are available - - * rho_mirror (float): reflectivity of mirrors, required for results post-processing + * num_hst (int): number of heliostats + * rho_mirror (float): reflectivity of mirrors, required for results post-processing * dni (float): the direct normal irradiance (W/m2), required to obtain performance of individual heliostat - * gen_vtk (bool): True - perform postprocessing for visualisation of each individual ray-tracing scene (each sun position), False - no postprocessing for visualisation + * gen_vtk (bool): True - perform postprocessing for visualisation of each individual ray-tracing scene (each sun position), False - no postprocessing for visualisation ``Return`` - * table (numpy array), the annual optical efficiency lookup table - * ANNUAL (numpy array), the annual output of each heliostat + * No return value (results files are created and written) """ YAML_IN = self.in_case(self.casedir, 'input.yaml') @@ -172,12 +158,12 @@ def run_annual(self, nd, nh, latitude, num_rays, num_hst,rho_mirror,dni, gen_vtk SOLSTICE_AZI, SOLSTICE_ELE=sun.convert_convention('solstice', AZI, ZENITH) # performance of individual heliostat is recorded - # TODO note, DNI is not varied in the simulation, + # TODO note, DNI is not varied in the simulation, # i.e. performance is not dni-weighted - ANNUAL=np.zeros((num_hst, 9)) + ANNUAL=np.zeros((num_hst, 9)) run=np.r_[0] - for i in range(len(case_list)): + for i in range(len(case_list)): c=int(case_list[i,0].astype(float)) if c not in run: azimuth=SOLSTICE_AZI[c-1] @@ -187,7 +173,7 @@ def run_annual(self, nd, nh, latitude, num_rays, num_hst,rho_mirror,dni, gen_vtk #else: # dni=0. - + sys.stderr.write("\n"+green('Sun position: %s \n'%c)) print('azimuth: %.2f'% azimuth, ', elevation: %.2f'%elevation) @@ -196,9 +182,9 @@ def run_annual(self, nd, nh, latitude, num_rays, num_hst,rho_mirror,dni, gen_vtk # run solstice if elevation<1.: # 1 degree efficiency_total=ufloat(0,0) - performance_hst=np.zeros((num_hst, 9)) + performance_hst=np.zeros((num_hst, 9)) else: - efficiency_total, performance_hst=self.run(azimuth, elevation, num_rays, rho_mirror, dni, folder=onesunfolder, gen_vtk=gen_vtk, printresult=False, verbose=verbose) + efficiency_total, performance_hst=self.run(azimuth, elevation, num_rays, rho_mirror, dni, folder=onesunfolder, gen_vtk=gen_vtk, printresult=False, system=system) sys.stderr.write(yellow("Total efficiency: {:f}\n".format(efficiency_total))) @@ -216,16 +202,12 @@ def run_annual(self, nd, nh, latitude, num_rays, num_hst,rho_mirror,dni, gen_vtk if c==float(val[0]): table[a+3,b+3]=efficiency_total.nominal_value - run=np.append(run,c) + run=np.append(run,c) - annual_title=np.array(['Q_solar','Q_cosine', 'Q_shade', 'Q_hst_abs', 'Q_block', 'Q_atm', 'Q_spil', 'Q_refl', 'Q_rcv_abs']) + annual_title=np.array(['Q_solar','Q_cosine', 'Q_shade', 'Q_hst_abs', 'Q_block', 'Q_atm', 'Q_spil', 'Q_refl', 'Q_rcv_abs']) ANNUAL=np.vstack((annual_title, ANNUAL)) - if verbose: - np.savetxt(self.casedir+'/lookup_table.csv', table, fmt='%s', delimiter=',') - np.savetxt(self.casedir+'/result-heliostats-annual-performance.csv', ANNUAL, fmt='%s', delimiter=',') - + np.savetxt(self.casedir+'/lookup_table.csv', table, fmt='%s', delimiter=',') + np.savetxt(self.casedir+'/result-heliostats-annual-performance.csv', ANNUAL, fmt='%s', delimiter=',') sys.stderr.write("\n"+green("Lookup table saved.\n")) sys.stderr.write(green("Completed successfully.\n"+"\n")) return table, ANNUAL - - diff --git a/solsticepy/output_motab.py b/solsticepy/output_motab.py index fb427bf..55f9e01 100644 --- a/solsticepy/output_motab.py +++ b/solsticepy/output_motab.py @@ -1,10 +1,9 @@ import numpy as np from datetime import datetime -import re def output_motab(table,savedir=None, title=None): ''' - output the .motab table file + output the .motab table field ''' f=open(savedir, 'w') f.write('#1\n') @@ -48,15 +47,15 @@ def output_motab(table,savedir=None, title=None): else: row_i=np.append(declination[i-1], table[t][2+i, 3:]) f.write(" ".join(map(str, row_i))) - f.write("\n") - + f.write("\n") + f.write("\n") f.write("\n") f.close() -def output_matadata_motab(table, field_type, aiming, n_helios, A_helio, eff_design, eff_annual, H_rcv, W_rcv, H_tower, Q_in_rcv, A_land, savedir=None, details_en=None): +def output_matadata_motab(table, field_type, aiming, n_helios, A_helio, eff_design, H_rcv, W_rcv, H_tower, Q_in_rcv, A_land, savedir=None, details_en=None): """Output the .motab file to work with the SolarTherm program ``Arguments`` @@ -66,7 +65,6 @@ def output_matadata_motab(table, field_type, aiming, n_helios, A_helio, eff_desi * n_helios (int): total number of heliostats * A_helio (float): area of each heliostat (m2) * eff_design (float): the optical efficiency at design point - * eff_annual (float): the annual optical efficiency (dni weighted) * H_rcv (float): height of the heliostat (m) * W_rcv (float): width of the heliostat (m) * H_tower (float), tower height (m) @@ -74,18 +72,18 @@ def output_matadata_motab(table, field_type, aiming, n_helios, A_helio, eff_desi * A_land (float): the total land area of the field (m2) * savedir (str): the directory to save this .motab file * details_en (dict): the key of the dict is the name of the breakdown of energy (loss), the value of the dict is a numpy array that contains the value of the energy (loss) with the corresponding declination and hour angles - + ``Return`` write the table(s) to the .motab file """ f=open(savedir, 'w') f.write('#1\n') f.write('#Comments: Field type: %s, Aiming Strategy: %s, Date:%s\n'%(field_type, aiming, datetime.now())) - f.write('#METALABELS,n_helios,A_helio,Eff_design,Eff_annual,H_rcv,W_rcv,H_tower, Q_in_rcv, A_land\n') - f.write('##METAUNITS,real,m2,real,real,m,m,m,W,m2\n') - f.write('#METADATA,%s,%s,%s,%s,%s,%s,%s,%s,%s\n'%(n_helios,A_helio,eff_design,eff_annual,H_rcv,W_rcv,H_tower,Q_in_rcv,A_land)) + f.write('#METALABELS,n_helios,A_helio,Eff_design,H_rcv,W_rcv,H_tower, Q_in_rcv, A_land\n') + f.write('##METAUNITS,real,m2,real,m,m,m,W,m2\n') + f.write('#METADATA,%s,%s,%s,%s,%s,%s,%s,%s\n'%(n_helios,A_helio,eff_design,H_rcv,W_rcv,H_tower,Q_in_rcv,A_land)) - # size of the lookup table + # size of the lookup table m=np.shape(table)[0]-2 n=np.shape(table)[1]-2 f.write('double optics(%s, %s)\n'%(m,n)) @@ -109,7 +107,7 @@ def output_matadata_motab(table, field_type, aiming, n_helios, A_helio, eff_desi for key in details_en: breakdown=key table=details_en[key] - # size of the lookup table + # size of the lookup table m=np.shape(table)[0]-2 n=np.shape(table)[1]-2 f.write('double %s(%s, %s)\n'%(key,m,n)) @@ -119,160 +117,36 @@ def output_matadata_motab(table, field_type, aiming, n_helios, A_helio, eff_desi else: row_i=np.append(declination[i-1], table[2+i, 3:]) f.write(" ".join(map(str, row_i))) - f.write("\n") - f.close() - - - -def output_matadata_motab_multi_aperture(TABLE, eff_design, eff_annual, A_land, H_tower, A_helio, n_helios_total, Q_in_rcv_total, num_aperture, Q_in_rcv, n_helios, H_rcv, W_rcv, savedir=None): - """Output the .motab file to work with the SolarTherm program - - ``Arguments`` - * TABLE (numpy array): the oelt that returned by design_crs.CRS.field_design_annual, listed by declination and hour angles, for each aperture - * eff_design (float): the optical efficiency at design point - * eff_annual (float): the annual optical efficiency (dni weighted) - * A_land (float): the total land area of the field (m2) - * H_tower (float): tower height (m) - * A_helio (float): area of each heliostat (m2) - * n_helios_total (int): total number of heliostats - * Q_in_rcv_total(float): the total required incident power on the receiver at design point (W) - * num_aperture (int): number of apertures in the multi-aperture configuration - * Q_in_rcv (list): a list of the required incident power on all the apertures at design point (W) - * n_helios (list): a list of number of heliostats that associated with each aperture - * H_rcv (list): a list of heights of all the apertures (m) - * W_rcv (list): a list of widthes of all the apertures (m) - * savedir (str): the directory to save this .motab file - * details_en (dict): the key of the dict is the name of the breakdown of energy (loss), the value of the dict is a numpy array that contains the value of the energy (loss) with the corresponding declination and hour angles - - ``Return`` - write the table(s) to the .motab file - """ - - labels='#METALABELS,eff_design,eff_annual,A_land,H_tower,A_helio,n_helios_total,Q_in_rcv_total,num_aperture' - units='##METAUNITS,real,real,m2,m,m2,real,W,real' - data='#METADATA,%s,%s,%s,%s,%s,%s,%s,%s'%(eff_design, eff_annual, A_land, H_tower, A_helio, n_helios_total, Q_in_rcv_total, num_aperture) - - for i in range(num_aperture): - labels+=',Q_in_rcv_%s,n_helios_%s,H_rcv_%s,W_rcv_%s'%(i, i, i, i) - units+=',W,real,m,m' - data+=',%s,%s,%s,%s'%(Q_in_rcv[i], n_helios[i], H_rcv[i], W_rcv[i]) - - labels+='\n' - units+='\n' - data+='\n' - - f=open(savedir, 'w') - f.write('#1\n') - f.write('#Comments: Field type: multi-aperture, Aiming Strategy: to-each-aperture-centre, Date:%s\n'%(datetime.now())) - f.write(labels) - f.write(units) - f.write(data) - f.write('\n') - - for i in range(num_aperture+1): - table=TABLE[i] - - # size of the lookup table - m=np.shape(table)[0]-2 - n=np.shape(table)[1]-2 - - if i==num_aperture: - f.write('double optical_efficiency_total(%s, %s)\n'%(m,n)) - else: - f.write('double optical_efficiency_aperture_%s(%s, %s)\n'%(i, m,n)) - - hour_angle=table[2, 3:] - declination=table[3:,2] - - for i in range(m): - if i ==0: - row_i=np.append(0, hour_angle) - else: - row_i=np.append(declination[i-1], table[2+i, 3:]) - - f.write(" ".join(map(str, row_i))) f.write("\n") - f.write("\n") - f.close() -def read_motab(filename, multi_aperture=False): +def read_motab(filename): with open(filename) as f: content=f.read().splitlines() f.close() res=content[4].split(',') - if multi_aperture: - - eff_des=float(res[1]) - eff_annual=float(res[2]) - A_land=float(res[3]) - n_helios=float(res[6]) - A_helio=float(res[5]) - Q_in_rcv=float(res[7]) - num_aperture=int(res[8]) - - OELT={} - Q_in_rcv_i=[] - n_helios_i=[] - H_rcv_i=[] - W_rcv_i=[] - for i in range(num_aperture+1): - if i Date: Sat, 26 Jun 2021 09:36:07 +1000 Subject: [PATCH 2/9] Implementation of beam down system Implementation of a Beam-Down system with 2D-crossed CPC (Compound Parabolic Concentrator) and secondary reflector (hyperboloid) with vertical axis and equivalent eccentricity in xOz and yOz planes. --- solsticepy/design_bd.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/solsticepy/design_bd.py b/solsticepy/design_bd.py index da88431..e964cb7 100644 --- a/solsticepy/design_bd.py +++ b/solsticepy/design_bd.py @@ -277,11 +277,11 @@ def heliostatfield(self, field, hst_rho, slope, hst_w, hst_h, tower_h, tower_r=0 self.hst_row=layout[:,-1].astype(float) Xmax = max(self.hst_pos[:,0]) - if (x_max>R1) and (Xmax>x_max*1.0): + if (x_max>R1) and (Xmax>x_max*1.05): select_hst=np.array([]) for i in range(len(self.hst_pos)): Xhelio = self.hst_pos[i,0] - if np.abs(Xhelio)R1) and (Ymax>y_max*1.0): + if (y_max>R1) and (Ymax>y_max*1.05): select_hst=np.array([]) for i in range(len(self.hst_pos)): Yhelio = self.hst_pos[i,1] - if np.abs(Yhelio) Date: Sat, 26 Jun 2021 13:04:49 +1000 Subject: [PATCH 3/9] Update gen_bd.py --- example/gen_bd.py | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/example/gen_bd.py b/example/gen_bd.py index 6a787ed..bcc1d53 100644 --- a/example/gen_bd.py +++ b/example/gen_bd.py @@ -47,17 +47,18 @@ num_hst = 300 pm.nd=5 pm.nh=5 -pm.H_rcv=10. -pm.W_rcv=10. +pm.H_rcv=rec_l +pm.W_rcv=rec_w +pm.W_helio=6.1 # ASTRI helio size +pm.H_helio=6.1 pm.dependent_par() pm.saveparam(casefolder) print(pm.fb) print(pm.H_tower) pm.field_type = 'surround' pm.Zhelio = 0. -pm.slope_error = 0. pm.R1 = 30. - +pm.lat = -27.16 #degree # enter the parameters for the beam-down components receiver='beam_down' @@ -74,7 +75,7 @@ slope_error = 2.e-3 # radian # parameters recalculated (pre-optimized before optimization) -secref_vert = None # np.array([[-15,25],[-15,-25],[15,-25],[15,25]]) +secref_vert=None # np.array([[-15,25],[-15,-25],[15,-25],[15,25]]) cpc_h=None # create the environment and scene From c9c055cd45abe41663370b593e04f309ee7e2a20 Mon Sep 17 00:00:00 2001 From: ClothildeCorsi Date: Sat, 26 Jun 2021 13:13:40 +1000 Subject: [PATCH 4/9] Update gen_bd.py correction --- example/gen_bd.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/example/gen_bd.py b/example/gen_bd.py index bcc1d53..93377e6 100644 --- a/example/gen_bd.py +++ b/example/gen_bd.py @@ -47,8 +47,8 @@ num_hst = 300 pm.nd=5 pm.nh=5 -pm.H_rcv=rec_l -pm.W_rcv=rec_w +pm.H_rcv=10. +pm.W_rcv=1.2 pm.W_helio=6.1 # ASTRI helio size pm.H_helio=6.1 pm.dependent_par() @@ -63,8 +63,8 @@ # enter the parameters for the beam-down components receiver='beam_down' # Polygon receiver -rec_w=1.2 -rec_l=10. +rec_w=pm.W_rcv +rec_l=pm.H_rcv rec_z=0. # 2D-crossed cpc with n faces n_CPC_faces=4 From 937735b8bad3bbe1cef9d044c29e667958bdb720 Mon Sep 17 00:00:00 2001 From: ClothildeCorsi Date: Sat, 26 Jun 2021 16:14:13 +1000 Subject: [PATCH 5/9] Update design_bd.py --- solsticepy/design_bd.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/solsticepy/design_bd.py b/solsticepy/design_bd.py index e964cb7..1802fb6 100644 --- a/solsticepy/design_bd.py +++ b/solsticepy/design_bd.py @@ -163,9 +163,9 @@ def receiversystem(self, receiver, rec_abs=1., rec_w=1.2, rec_l=10., rec_z=0., r self.receiver=receiver self.rec_abs=rec_abs - assert cpc_nfaces > 2, f'The number of faces for the CPC should be minimum 3, and it is {cpc_nfaces=}' + assert cpc_nfaces > 2, 'The number of faces for the CPC should be minimum 3, and it is {cpc_nfaces=}' if secref_fratio != None: - assert 0.5<=secref_fratio<1, f'The ratio of the hyperbole distances to foci and apex should be between 0.5 and 1, and it is {secref_fratio=}' + assert 0.5<=secref_fratio<1, 'The ratio of the hyperbole distances to foci and apex should be between 0.5 and 1, and it is {secref_fratio=}' rec_rad_ref, rec_rad = self.cpcradius(rec_w, rec_l, cpc_nfaces) # Calculate the maximum height of the CPC if not defined From b1fab1736549f3b128b58f64d54e2d1b6131a79b Mon Sep 17 00:00:00 2001 From: ClothildeCorsi Date: Sat, 26 Jun 2021 20:35:11 +1000 Subject: [PATCH 6/9] parameters change --- example/gen_bd.py | 4 ++-- solsticepy/design_bd.py | 4 ++-- solsticepy/input.py | 2 +- solsticepy/master.py | 2 +- 4 files changed, 6 insertions(+), 6 deletions(-) diff --git a/example/gen_bd.py b/example/gen_bd.py index 93377e6..7273289 100644 --- a/example/gen_bd.py +++ b/example/gen_bd.py @@ -68,7 +68,7 @@ rec_z=0. # 2D-crossed cpc with n faces n_CPC_faces=4 -rec_grid=200 +rec_grid=50 n_Z=30 # Secondary refector 'hyperboloid' refl_sec = 0.95 @@ -93,7 +93,7 @@ bd.yaml(sunshape=pm.sunshape,csr=pm.crs,half_angle_deg=pm.half_angle_deg,std_dev=pm.std_dev) -oelt, A_land=bd.field_design_annual(dni_des=900., num_rays=int(1e5), nd=pm.nd, nh=pm.nh, weafile=weafile, method=1, Q_in_des=pm.Q_in_rcv, n_helios=None, zipfiles=False, gen_vtk=True, plot=False) +oelt, A_land=bd.field_design_annual(dni_des=900., num_rays=int(5e7), nd=pm.nd, nh=pm.nh, weafile=weafile, method=1, Q_in_des=pm.Q_in_rcv, n_helios=None, zipfiles=False, gen_vtk=True, plot=False) if (A_land==0): diff --git a/solsticepy/design_bd.py b/solsticepy/design_bd.py index 1802fb6..15f8a86 100644 --- a/solsticepy/design_bd.py +++ b/solsticepy/design_bd.py @@ -374,7 +374,7 @@ def field_design_annual(self, dni_des, num_rays, nd, nh, weafile, method, Q_in_ performance_hst=np.zeros((nhst, 9)) efficiency_hst=np.zeros(nhst) else: - efficiency_total, performance_hst=self.master.run(azimuth, elevation, num_rays, self.hst_rho, dni, folder=onesunfolder, gen_vtk=False, printresult=False, system='beamdown') + efficiency_total, performance_hst=self.master.run(azimuth, elevation, num_rays, self.hst_rho, dni, folder=onesunfolder, gen_vtk=gen_vtk, printresult=False, system='beamdown') efficiency_hst=performance_hst[:,-1]/performance_hst[:,0] @@ -406,7 +406,7 @@ def field_design_annual(self, dni_des, num_rays, nd, nh, weafile, method, Q_in_ azi_des, ele_des=self.sun.convert_convention('solstice', azi, zen) sys.stderr.write("\n"+green('Design Point: \n')) - efficiency_total, performance_hst_des=self.master.run(azi_des, ele_des, num_rays, self.hst_rho, dni_des, folder=designfolder, gen_vtk=True, printresult=False, system='beamdown') + efficiency_total, performance_hst_des=self.master.run(azi_des, ele_des, num_rays, self.hst_rho, dni_des, folder=designfolder, gen_vtk=gen_vtk, printresult=False, system='beamdown') Qin=performance_hst_des[:,-1] Qsolar=performance_hst_des[0,0] diff --git a/solsticepy/input.py b/solsticepy/input.py index 0c172a6..31b5231 100644 --- a/solsticepy/input.py +++ b/solsticepy/input.py @@ -124,7 +124,7 @@ def simulation(self): ''' self.n_row_oelt=5 self.n_col_oelt=5 - self.n_rays=int(5e6) + self.n_rays=int(5e7) self.n_procs=1 self.casedir='.' self.method=1 # 1 - design the field based on the Q_in_rcv diff --git a/solsticepy/master.py b/solsticepy/master.py index de95851..a95c634 100755 --- a/solsticepy/master.py +++ b/solsticepy/master.py @@ -119,7 +119,7 @@ def run(self, azimuth, elevation, num_rays, rho_mirror, dni, folder, gen_vtk=Tru else: if system=='beamdown': - eta, performance_hst=process_raw_results(self.in_case(folder, 'simul'), folder, rho_mirror, dni, verbose=True, num_virt=2) + eta, performance_hst=process_raw_results(self.in_case(folder, 'simul'), folder, rho_mirror, dni, verbose=False, num_virt=2) else: eta, performance_hst=process_raw_results(self.in_case(folder, 'simul'), folder, rho_mirror, dni) From 8754008b9fe9da7f0e88a015c7745b5ec81fcd73 Mon Sep 17 00:00:00 2001 From: ClothildeCorsi Date: Sat, 26 Jun 2021 21:39:17 +1000 Subject: [PATCH 7/9] add beam-down parameters in input module + change cpc_h to ratio_cpc_h --- example/gen_bd.py | 11 +++++------ solsticepy/design_bd.py | 13 +++++-------- solsticepy/input.py | 14 ++++++++++++++ 3 files changed, 24 insertions(+), 14 deletions(-) diff --git a/example/gen_bd.py b/example/gen_bd.py index 7273289..e4c177d 100644 --- a/example/gen_bd.py +++ b/example/gen_bd.py @@ -39,7 +39,8 @@ theta_deg=20. # acceptance half angle of the CPC in degree rim_angle = 45. # rim angle of the heliostat field in the xOz plan in degree foci_ratio = 0.6 # ratio of the foci distance and apex distance to the origin [1/2,1] - +ratio_cpc_h=1. # must be inferior to 1 +rec_z=0. # fixed parameters # ========= @@ -56,7 +57,7 @@ print(pm.fb) print(pm.H_tower) pm.field_type = 'surround' -pm.Zhelio = 0. +pm.Z_helio = 0. pm.R1 = 30. pm.lat = -27.16 #degree @@ -65,7 +66,6 @@ # Polygon receiver rec_w=pm.W_rcv rec_l=pm.H_rcv -rec_z=0. # 2D-crossed cpc with n faces n_CPC_faces=4 rec_grid=50 @@ -76,7 +76,6 @@ # parameters recalculated (pre-optimized before optimization) secref_vert=None # np.array([[-15,25],[-15,-25],[15,-25],[15,25]]) -cpc_h=None # create the environment and scene # ========= @@ -85,7 +84,7 @@ bd=BD(latitude=pm.lat, casedir=casefolder) -bd.receiversystem(receiver=receiver, rec_abs=float(pm.alpha_rcv), rec_w=float(rec_w), rec_l=float(rec_l), rec_z=float(rec_z), rec_grid=int(rec_grid), cpc_nfaces=int(n_CPC_faces), cpc_theta_deg=float(theta_deg), cpc_h=cpc_h, cpc_nZ=float(n_Z), field_rim_angle=float(rim_angle), aim_z=float(pm.H_tower), secref_fratio=foci_ratio, refl_sec=float(refl_sec), slope_error=float(slope_error), secref_vert = secref_vert) +bd.receiversystem(receiver=receiver, rec_abs=float(pm.alpha_rcv), rec_w=float(rec_w), rec_l=float(rec_l), rec_z=float(rec_z), rec_grid=int(rec_grid), cpc_nfaces=int(n_CPC_faces), cpc_theta_deg=float(theta_deg), ratio_cpc_h=ratio_cpc_h, cpc_nZ=float(n_Z), field_rim_angle=float(rim_angle), aim_z=float(pm.H_tower), secref_fratio=foci_ratio, refl_sec=float(refl_sec), slope_error=float(slope_error), secref_vert = secref_vert) bd.heliostatfield(field=pm.field_type, hst_rho=pm.rho_helio, slope=pm.slope_error, hst_w=pm.W_helio, hst_h=pm.H_helio, tower_h=pm.H_tower, tower_r=pm.R_tower, hst_z=pm.Z_helio, num_hst=num_hst, R1=pm.R1, fb=pm.fb, dsep=pm.dsep, x_max=150., y_max=150.) @@ -93,7 +92,7 @@ bd.yaml(sunshape=pm.sunshape,csr=pm.crs,half_angle_deg=pm.half_angle_deg,std_dev=pm.std_dev) -oelt, A_land=bd.field_design_annual(dni_des=900., num_rays=int(5e7), nd=pm.nd, nh=pm.nh, weafile=weafile, method=1, Q_in_des=pm.Q_in_rcv, n_helios=None, zipfiles=False, gen_vtk=True, plot=False) +oelt, A_land=bd.field_design_annual(dni_des=900., num_rays=int(5e7), nd=pm.nd, nh=pm.nh, weafile=weafile, method=1, Q_in_des=pm.Q_in_rcv, n_helios=None, zipfiles=False, gen_vtk=False, plot=False) if (A_land==0): diff --git a/solsticepy/design_bd.py b/solsticepy/design_bd.py index 15f8a86..b26196a 100644 --- a/solsticepy/design_bd.py +++ b/solsticepy/design_bd.py @@ -128,7 +128,7 @@ def zhyperboloid(self, x_hyper, y_hyper, vertex_dist, foci_dist): z_hyper = b_para*np.sqrt((x_hyper**2+y_hyper**2)/a_para_sqrt+1) - z0_hyper+g_para/2. return z_hyper - def receiversystem(self, receiver, rec_abs=1., rec_w=1.2, rec_l=10., rec_z=0., rec_grid=200, cpc_nfaces=4, cpc_theta_deg=20., cpc_h=None, + def receiversystem(self, receiver, rec_abs=1., rec_w=1.2, rec_l=10., rec_z=0., rec_grid=200, cpc_nfaces=4, cpc_theta_deg=20., ratio_cpc_h=1., cpc_nZ=20., field_rim_angle=30., aim_z=62., secref_fratio=None, refl_sec=0.95, slope_error=0.0, secref_vert=np.array([[-15,25],[-15,-25],[15,-25],[15,25]])): ''' @@ -148,7 +148,7 @@ def receiversystem(self, receiver, rec_abs=1., rec_w=1.2, rec_l=10., rec_z=0., r # Arguments for Compound Parabolic Concentrator (CPC) (7) cpc_nfaces : int, number of faces of the CPC (8) cpc_theta_deg : float, acceptance angle of CPC (deg) - (9) cpc_h : float, vertical distance of the maximum and minimum point of the parabola + (9) ratio_cpc_h : float, ratio of critical CPC height calculated with theta_deg in local coordinate system of the parabola in the xOy plan WARNING: cpc_h is different from the total height of the CPC (10) cpc_nZ : int, number of number of incrementation for the clipping polygon for the construction of each CPC face @@ -169,10 +169,9 @@ def receiversystem(self, receiver, rec_abs=1., rec_w=1.2, rec_l=10., rec_z=0., r rec_rad_ref, rec_rad = self.cpcradius(rec_w, rec_l, cpc_nfaces) # Calculate the maximum height of the CPC if not defined - if cpc_h is None: - cpc_theta = cpc_theta = cpc_theta_deg*np.pi/180. - cpc_h = rec_rad_ref*(1+1/np.sin(cpc_theta))/np.tan(cpc_theta) - print('height of cpc: ', cpc_h) + cpc_theta = cpc_theta_deg*np.pi/180. + cpc_h = rec_rad_ref*(1+1/np.sin(cpc_theta))/np.tan(cpc_theta)*ratio_cpc_h + print('height of cpc: ', cpc_h) assert aim_z > (cpc_h+rec_z), 'The imaginary foci of the hyperbol is lower than its real foci' @@ -204,13 +203,11 @@ def receiversystem(self, receiver, rec_abs=1., rec_w=1.2, rec_l=10., rec_z=0., r else: y_inter = x_inter print('x_inter: ', x_inter) - print('THEY SHOULD BE EQUAL:', y_inter) secref_vert = np.array([[-x_inter,y_inter],[-x_inter,-y_inter],[x_inter,-y_inter],[x_inter,y_inter]]) # calculate the field maxium radius along x and y axis self.x_max = aim_z*np.tan(field_rim_angle*np.pi/180.) y_inter = max(secref_vert[:,1]) - print('THEY SHOULD BE EQUAL:', y_inter) z_hyper = secref_z + self.zhyperboloid( 0.0, y_inter, vertex_dist, foci_dist) phy = np.arctan(y_inter/(aim_z-z_hyper)) self.y_max = aim_z*np.tan(phy) diff --git a/solsticepy/input.py b/solsticepy/input.py index 31b5231..f04f404 100644 --- a/solsticepy/input.py +++ b/solsticepy/input.py @@ -19,6 +19,7 @@ def __init__(self): self.Sun() self.Heliostat() self.Receiver() + self.BeamDown() self.dependent_par() def Sun(self): @@ -113,6 +114,14 @@ def Receiver(self): self.X_rcv=0. # receiver location self.Y_rcv=0. + def BeamDown(self): + ''' + ''' + self.theta_deg=20. # acceptance half angle of CPC + self.ratio_cpc_h=1. # ratio of CPC critical height [0.5,1] + self.field_rim_angle=45. # field rim angle + self.secref_fratio=0.6 # ratio of the foci distance and apex distance to the origin [0.5,1] + self.rec_z=0. def simulation(self): ''' @@ -173,6 +182,11 @@ def saveparam(self, savedir): ['rcv_type', self.rcv_type, '-'], ['H_rcv', self.H_rcv, 'm'], ['W_rcv', self.W_rcv, 'm'], + ['rec_z', self.rec_z, 'm'], + ['theta_deg', self.theta_deg, 'deg'], + ['ratio_cpc_h', self.ratio_cpc_h, '-'], + ['field_rim_angle', self.field_rim_angle, '-'], + ['secref_fratio', self.secref_fratio, '-'], ['tilt_rcv', self.tilt_rcv, 'deg'], ['alpha_rcv', self.alpha_rcv, '-'], ['n_H_rcv', self.n_H_rcv, '-'], From b7f58c96b485058c2e11c80ed60b8b42cca4200e Mon Sep 17 00:00:00 2001 From: ClothildeCorsi Date: Sat, 26 Jun 2021 22:50:22 +1000 Subject: [PATCH 8/9] correction theta name --- example/gen_bd.py | 4 ++-- solsticepy/input.py | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/example/gen_bd.py b/example/gen_bd.py index e4c177d..e93f987 100644 --- a/example/gen_bd.py +++ b/example/gen_bd.py @@ -36,7 +36,7 @@ # Variables pm.H_tower=65. # tower height or vertical distance to aiming point (located at the center of xOy plan) -theta_deg=20. # acceptance half angle of the CPC in degree +cpc_theta_deg=20. # acceptance half angle of the CPC in degree rim_angle = 45. # rim angle of the heliostat field in the xOz plan in degree foci_ratio = 0.6 # ratio of the foci distance and apex distance to the origin [1/2,1] ratio_cpc_h=1. # must be inferior to 1 @@ -84,7 +84,7 @@ bd=BD(latitude=pm.lat, casedir=casefolder) -bd.receiversystem(receiver=receiver, rec_abs=float(pm.alpha_rcv), rec_w=float(rec_w), rec_l=float(rec_l), rec_z=float(rec_z), rec_grid=int(rec_grid), cpc_nfaces=int(n_CPC_faces), cpc_theta_deg=float(theta_deg), ratio_cpc_h=ratio_cpc_h, cpc_nZ=float(n_Z), field_rim_angle=float(rim_angle), aim_z=float(pm.H_tower), secref_fratio=foci_ratio, refl_sec=float(refl_sec), slope_error=float(slope_error), secref_vert = secref_vert) +bd.receiversystem(receiver=receiver, rec_abs=float(pm.alpha_rcv), rec_w=float(rec_w), rec_l=float(rec_l), rec_z=float(rec_z), rec_grid=int(rec_grid), cpc_nfaces=int(n_CPC_faces), cpc_theta_deg=float(cpc_theta_deg), ratio_cpc_h=ratio_cpc_h, cpc_nZ=float(n_Z), field_rim_angle=float(rim_angle), aim_z=float(pm.H_tower), secref_fratio=foci_ratio, refl_sec=float(refl_sec), slope_error=float(slope_error), secref_vert = secref_vert) bd.heliostatfield(field=pm.field_type, hst_rho=pm.rho_helio, slope=pm.slope_error, hst_w=pm.W_helio, hst_h=pm.H_helio, tower_h=pm.H_tower, tower_r=pm.R_tower, hst_z=pm.Z_helio, num_hst=num_hst, R1=pm.R1, fb=pm.fb, dsep=pm.dsep, x_max=150., y_max=150.) diff --git a/solsticepy/input.py b/solsticepy/input.py index f04f404..8393119 100644 --- a/solsticepy/input.py +++ b/solsticepy/input.py @@ -117,7 +117,7 @@ def Receiver(self): def BeamDown(self): ''' ''' - self.theta_deg=20. # acceptance half angle of CPC + self.cpc_theta_deg=20. # acceptance half angle of CPC self.ratio_cpc_h=1. # ratio of CPC critical height [0.5,1] self.field_rim_angle=45. # field rim angle self.secref_fratio=0.6 # ratio of the foci distance and apex distance to the origin [0.5,1] @@ -183,7 +183,7 @@ def saveparam(self, savedir): ['H_rcv', self.H_rcv, 'm'], ['W_rcv', self.W_rcv, 'm'], ['rec_z', self.rec_z, 'm'], - ['theta_deg', self.theta_deg, 'deg'], + ['cpc_theta_deg', self.cpc_theta_deg, 'deg'], ['ratio_cpc_h', self.ratio_cpc_h, '-'], ['field_rim_angle', self.field_rim_angle, '-'], ['secref_fratio', self.secref_fratio, '-'], From e284b3022f310d52b9bf3addbb52885c9a7c36ed Mon Sep 17 00:00:00 2001 From: ClothildeCorsi Date: Sat, 26 Jun 2021 23:31:02 +1000 Subject: [PATCH 9/9] change back --- example/gen_bd.py | 4 ++-- solsticepy/input.py | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/example/gen_bd.py b/example/gen_bd.py index e93f987..e4c177d 100644 --- a/example/gen_bd.py +++ b/example/gen_bd.py @@ -36,7 +36,7 @@ # Variables pm.H_tower=65. # tower height or vertical distance to aiming point (located at the center of xOy plan) -cpc_theta_deg=20. # acceptance half angle of the CPC in degree +theta_deg=20. # acceptance half angle of the CPC in degree rim_angle = 45. # rim angle of the heliostat field in the xOz plan in degree foci_ratio = 0.6 # ratio of the foci distance and apex distance to the origin [1/2,1] ratio_cpc_h=1. # must be inferior to 1 @@ -84,7 +84,7 @@ bd=BD(latitude=pm.lat, casedir=casefolder) -bd.receiversystem(receiver=receiver, rec_abs=float(pm.alpha_rcv), rec_w=float(rec_w), rec_l=float(rec_l), rec_z=float(rec_z), rec_grid=int(rec_grid), cpc_nfaces=int(n_CPC_faces), cpc_theta_deg=float(cpc_theta_deg), ratio_cpc_h=ratio_cpc_h, cpc_nZ=float(n_Z), field_rim_angle=float(rim_angle), aim_z=float(pm.H_tower), secref_fratio=foci_ratio, refl_sec=float(refl_sec), slope_error=float(slope_error), secref_vert = secref_vert) +bd.receiversystem(receiver=receiver, rec_abs=float(pm.alpha_rcv), rec_w=float(rec_w), rec_l=float(rec_l), rec_z=float(rec_z), rec_grid=int(rec_grid), cpc_nfaces=int(n_CPC_faces), cpc_theta_deg=float(theta_deg), ratio_cpc_h=ratio_cpc_h, cpc_nZ=float(n_Z), field_rim_angle=float(rim_angle), aim_z=float(pm.H_tower), secref_fratio=foci_ratio, refl_sec=float(refl_sec), slope_error=float(slope_error), secref_vert = secref_vert) bd.heliostatfield(field=pm.field_type, hst_rho=pm.rho_helio, slope=pm.slope_error, hst_w=pm.W_helio, hst_h=pm.H_helio, tower_h=pm.H_tower, tower_r=pm.R_tower, hst_z=pm.Z_helio, num_hst=num_hst, R1=pm.R1, fb=pm.fb, dsep=pm.dsep, x_max=150., y_max=150.) diff --git a/solsticepy/input.py b/solsticepy/input.py index 8393119..f04f404 100644 --- a/solsticepy/input.py +++ b/solsticepy/input.py @@ -117,7 +117,7 @@ def Receiver(self): def BeamDown(self): ''' ''' - self.cpc_theta_deg=20. # acceptance half angle of CPC + self.theta_deg=20. # acceptance half angle of CPC self.ratio_cpc_h=1. # ratio of CPC critical height [0.5,1] self.field_rim_angle=45. # field rim angle self.secref_fratio=0.6 # ratio of the foci distance and apex distance to the origin [0.5,1] @@ -183,7 +183,7 @@ def saveparam(self, savedir): ['H_rcv', self.H_rcv, 'm'], ['W_rcv', self.W_rcv, 'm'], ['rec_z', self.rec_z, 'm'], - ['cpc_theta_deg', self.cpc_theta_deg, 'deg'], + ['theta_deg', self.theta_deg, 'deg'], ['ratio_cpc_h', self.ratio_cpc_h, '-'], ['field_rim_angle', self.field_rim_angle, '-'], ['secref_fratio', self.secref_fratio, '-'],