models/solar/solar_field.a4l

(*
Model of a parabolic trough solar thermal collector field, of the type
of SEGS VI plant, closely following the approach of Patnode (2006).
https://www.nrel.gov/analysis/sam/pdfs/thesis_patnode06.pdf

First version by Vikram Kaadam (2011) as part of GSOC2011.
Second version by John Pye (2012).
*)
REQUIRE "atoms.a4l";
REQUIRE "solar/solar_types.a4l";
REQUIRE "johnpye/thermo_types.a4c";
REQUIRE "solar/tracker.a4l";
REQUIRE "solar/therminol.a4c";

(* refinement base, used in the parabolic_trough model... *)
MODEL absorber_loss_base;
END absorber_loss_base;

(*
	Heat loss per metre length of absorber tube, including convective, conductive 
	and radiative losses.	
*) 
MODEL absorber_loss(
	type IS_A symbol_constant;
	T_i WILL_BE temperature;
	T_o WILL_BE temperature;
	DNI WILL_BE intensity;
) REFINES absorber_loss_base;

	a[0..3] IS_A real_constant;
	b[0,1] IS_A real_constant;
	
	SELECT(type)
	CASE 'vacuum':
		(* an intact HCE, evacuated air at 0.0001 Torr *)
		a[0] :== -2.247372E+01 {W/m^2};
		a[1] :==  8.374490E-01 {W/m^2/K};
		a[2] :==  0.00         {W/m^2/K^2};
		a[3] :==  4.620143E-06 {W/m^2/K^3};
		b[0] :==  6.983190E-02	{m};
		b[1] :==  9.312703E-08 {m/K^2};
	CASE 'hydrogen':
		(* intact HCE into which hydrogen has diffused to a pressure off 1 Torr *)
		a[0] :== -2.247372E+01 {W/m^2};
		a[1] :==  8.374490E-01 {W/m^2/K};
		a[2] :==  0.00         {W/m^2/K^2};
		a[3] :==  4.620143E-06 {W/m^2/K^3};
		b[0] :==  6.983190E-02	{m};
		b[1] :==  9.312703E-08 {m/K^2};
	CASE 'air':
		(* a broken glass envelope; contents will be air at ambient pressure *)
		a[0] :== -2.247372E+01 {W/m^2};
		a[1] :==  8.374490E-01 {W/m^2/K};
		a[2] :==  0.00         {W/m^2/K^2};
		a[3] :==  4.620143E-06 {W/m^2/K^3};
		b[0] :==  6.983190E-02	{m};
		b[1] :==  9.312703E-08 {m/K^2};		
	END SELECT;
	
	Qd_loss IS_A power_per_length;  

	(* heat loss calculated as integral wrt temperature, Patnode eq 2.18. *)
	Qd_loss * (T_o - T_i) = SUM[a[i]*(T_o^(i+1) - T_i^(i+1))/(i+1) | i IN [0..3]] + DNI * (b[0] + b[1]/3*(T_o^3 - T_i^3));
	(* NOTE above eq assumes linear temperature rise with position? is that suff valid? *)
	(* NOTE above eq assumes ambient temperature of 25 C (Patnode sect 2.3.2) -- should adjust for changes in T_amb?? *)
	(* NOTE also the equation in Lippke 1995 for bare tubes, which has a different form and requires wind speed *)
	(* NOTE that Lippke says that end loss must be included in the above equation (eq 7, p. 11) *)
END absorber_loss;

MODEL absorber_loss_test;
	type IS_A symbol_constant;
	type :== 'vacuum';
	T_i, T_o IS_A temperature;
	DNI IS_A intensity;
	abso IS_A absorber_loss(type,T_i,T_o,DNI);
	Qd_loss ALIASES abso.Qd_loss;
METHODS 
METHOD on_load;
	FIX DNI, T_i, T_o;
	DNI := 1000 {W/m^2};
	T_i := 300 {K} +  273.15 {K};
	T_o := 400 {K} + 237.15 {K};
END on_load;
END absorber_loss_test;


(*
	Heat absorbed by the solar receiver of cylindrical shape
	
	TODO implement option to set field orientation EW/NS.
*) 
MODEL parabolic_trough;
	(* field geometry *)
	W_apert    "collector aperture width" IS_A distance;
	L_spacing  "'length of spacing' between troughs" IS_A distance; (* clarification required, assume centre-to-centre? *)
	L_SCA      "length of a single solar collector assembly (SCA)" IS_A distance;
	N_SCA      "total number of solar collector assemblies" IS_A factor;
	f          "focal length of collector" IS_A distance;
	
	A_field    IS_A area;
	A_field = N_SCA * L_SCA * W_apert;
	
	(* tracker and sun position *)
	sun    IS_A solar_tracker_single_horiz;
	t          "local time (standard, not DST correction)" ALIASES sun.t;

	Qdd_abs    "solar radiation absorbed by receiver tubes per aperture area" IS_A power_per_area;
	DNI        "direct normal irradiance" IS_A intensity;
	IAM        "incidence angle modifier" IS_A factor;
	eta_shading IS_A fraction;
	eta_endloss IS_A fraction;
	eta_field  "averaged field efficiency" IS_A fraction;
	eta_HCE    "averaged heat collection element efficiency" IS_A fraction;
	avail_SF   "solar field availability - on-sun portion" IS_A fraction;

	(* constants that define the location -- can we do this better?? *)
	(* Kramer Junction *)
	sun.loc.latitude :== 37.21 {deg};
	sun.loc.longitude :== -117.022 {deg};
	sun.loc.elevation :== 755 {ft};

	(* OPTICAL PERFORMANCE *)
	
	(* absorbed solar radiation, Patnode eq 2.1. factored the cos(theta) into IAM. *)	
	absheat_eq: Qdd_abs = DNI * IAM * eta_shading * eta_endloss * eta_field * eta_HCE * avail_SF;
	(* NOTE that eta_field and eta_HCE are given fixed values here, though their 
	components factors are shown in Patnode sect 2.2.5 *)

	(* TODO where is Patnode eq 2.8 ? looks like an error? *)
	
	(* incidence angle modifier, Patnode eq 2.9 *)
	(* TODO check... do these values assume angles in degrees? *)
	IAM * cos(sun.theta) = cos(sun.theta) + 8.84e-4 * sun.theta - 5.369e-5 * sun.theta^2;
	
	(* row shading (Stuetzle), Patnode eq 2.12 *)
	eta_shading_1, eta_shading_2 IS_A factor; (* before being normalised to [0,1] *)
	shading_1: eta_shading_1 = (L_spacing/W_apert) * (cos(sun.zenith)/cos(sun.theta));
	(* limit first to max of 1, then min of 0 *)
	shading_2: eta_shading_2 = 1. + 0.5*((eta_shading_1 - 1) - abs(eta_shading_1 - 1));
	shading_3: eta_shading = (eta_shading_2 + abs(eta_shading_2))/2;

	(* end loss (Lippke), Patnode eq 2.13 *)
	eta_endloss_1, eta_endloss_2 IS_A factor;
	endloss_1: eta_endloss_1 = 1 - f * tan(sun.theta) / L_SCA;
	(* limit first to max of 1, then min of 0 *)
	endloss_2: eta_endloss_2 = 1. + 0.5*((eta_endloss_1 - 1) - abs(eta_endloss_1 - 1));
	endloss_3: eta_endloss = (eta_endloss_2 + abs(eta_endloss_2))/2;

	(* LOSSES FROM HCE *)
	
	T_i, T_o IS_A temperature;
	hce_types IS_A set OF symbol_constant;
	hce_types :== ['air','vacuum','hydrogen'];

	abso[hce_types] IS_A absorber_loss_base;
	hce_frac[hce_types] IS_A fraction;
	FOR i IN hce_types CREATE
		abso[i] IS_REFINED_TO absorber_loss(i, T_i, T_o, DNI);
	END FOR;
	
	Qdd_loss_recv "HCE losses per aperture area" IS_A power_per_area;
	(* receiver heat loss, Patnode eq 2.19 *)
	Qdd_loss_recv * W_apert = SUM[abso[i].Qd_loss * hce_frac[i] | i IN hce_types];
	
	(* LOSSES FROM FIELD PIPING *)
	
	T_amb "ambient temperature" IS_A temperature;
	DT "temperature difference ambient to field piping" IS_A delta_temperature;
	(* average external temp difference, Patnode eq 2.21 *)
	T_amb + DT = 0.5 * (T_i + T_o);

	Qdd_loss_pipe "losses from pipework, per solar field aperture area" IS_A power_per_area;
	
	(* solar field piping loses, Patnode eq 2.20 *)
	Qdd_loss_pipe = 0.01693{W/m^2/K}*DT - 0.0001683{W/m^2/K^2}*(DT^2) + 0.78e-7{W/m^2/K^3}*(DT^3);	
	(* typically 10 W/m2 or less, apparently *)
	
	(* NET ENERGY GAIN IN HTF *)
	
	Vdot_i IS_A volume_rate;
	Qdd_fluid IS_A power_per_area;
	(* net absorbed energy, Patnode eq 2.22 *)
	Qdd_fluid = Qdd_abs - Qdd_loss_pipe - Qdd_loss_recv;
	
	inlet, outlet IS_A therminol;
	inlet.T, T_i ARE_THE_SAME;
	outlet.T, T_o ARE_THE_SAME;
	
	Q_net IS_A energy_rate;
	Q_net = Qdd_fluid * A_field;
	
	(* first-law energy balance, Patnode eq 2.23 *)
	(outlet.h - inlet.h) * inlet.rho * Vdot_i = Q_net;

METHODS
	METHOD specify;
		RUN sun.specify;
		FIX DNI, L_spacing;
		FIX W_apert, f, L_SCA, eta_field, eta_HCE, avail_SF;
		
		FIX inlet.T;
		FIX T_amb;			
		FIX N_SCA;
		FIX Vdot_i;
		FIX hce_frac[hce_types];		
	END specify;

	METHOD values;		
		L_spacing := 13 {m};
		W_apert := 4.83 {m};
		f := 5 {m};
		L_SCA := 50 {m};
		N_SCA := 256;

		DNI := 1000 {W/m^2};

		(* tracker, including sun position and geographical location *)
		RUN sun.values;
		sun.offset.tz := -8{h};
		sun.offset.h := 0;
		t := 43.5 {d};

		(* fixed field efficiency based on values in Patnode Table 2.1 *)
		eta_field := 0.857;		
		eta_HCE := 0.832;
		avail_SF := 1;

		inlet.T := 60{K} + 273.15 {K};
		T_amb := 30{K} + 273.15 {K};
		Vdot_i := 400{m^3/hour};

		hce_frac['air'] := 1.0;
		hce_frac['vacuum'] := 0.0;
		hce_frac['hydrogen'] := 0.0;

		(* initial guess, needed for solving OK  *) 
		T_o := 200 {K} + 273.15{K};
	END values;

	METHOD on_load;
		RUN specify;
		RUN values;
	END on_load;	
END parabolic_trough;
1	(*
2	Model of a parabolic trough solar thermal collector field, of the type
3	of SEGS VI plant, closely following the approach of Patnode (2006).
4	https://www.nrel.gov/analysis/sam/pdfs/thesis_patnode06.pdf
5
6	First version by Vikram Kaadam (2011) as part of GSOC2011.
7	Second version by John Pye (2012).
8	*)
9	REQUIRE "atoms.a4l";
10	REQUIRE "solar/solar_types.a4l";
11	REQUIRE "johnpye/thermo_types.a4c";
12	REQUIRE "solar/tracker.a4l";
13	REQUIRE "solar/therminol.a4c";
14
15	(* refinement base, used in the parabolic_trough model... *)
16	MODEL absorber_loss_base;
17	END absorber_loss_base;
18
19	(*
20	Heat loss per metre length of absorber tube, including convective, conductive
21	and radiative losses.
22	*)
23	MODEL absorber_loss(
24	type IS_A symbol_constant;
25	T_i WILL_BE temperature;
26	T_o WILL_BE temperature;
27	DNI WILL_BE intensity;
28	) REFINES absorber_loss_base;
29
30	a[0..3] IS_A real_constant;
31	b[0,1] IS_A real_constant;
32
33	SELECT(type)
34	CASE 'vacuum':
35	(* an intact HCE, evacuated air at 0.0001 Torr *)
36	a[0] :== -2.247372E+01 {W/m^2};
37	a[1] :== 8.374490E-01 {W/m^2/K};
38	a[2] :== 0.00 {W/m^2/K^2};
39	a[3] :== 4.620143E-06 {W/m^2/K^3};
40	b[0] :== 6.983190E-02 {m};
41	b[1] :== 9.312703E-08 {m/K^2};
42	CASE 'hydrogen':
43	(* intact HCE into which hydrogen has diffused to a pressure off 1 Torr *)
44	a[0] :== -2.247372E+01 {W/m^2};
45	a[1] :== 8.374490E-01 {W/m^2/K};
46	a[2] :== 0.00 {W/m^2/K^2};
47	a[3] :== 4.620143E-06 {W/m^2/K^3};
48	b[0] :== 6.983190E-02 {m};
49	b[1] :== 9.312703E-08 {m/K^2};
50	CASE 'air':
51	(* a broken glass envelope; contents will be air at ambient pressure *)
52	a[0] :== -2.247372E+01 {W/m^2};
53	a[1] :== 8.374490E-01 {W/m^2/K};
54	a[2] :== 0.00 {W/m^2/K^2};
55	a[3] :== 4.620143E-06 {W/m^2/K^3};
56	b[0] :== 6.983190E-02 {m};
57	b[1] :== 9.312703E-08 {m/K^2};
58	END SELECT;
59
60	Qd_loss IS_A power_per_length;
61
62	(* heat loss calculated as integral wrt temperature, Patnode eq 2.18. *)
63	Qd_loss * (T_o - T_i) = SUM[a[i](T_o^(i+1) - T_i^(i+1))/(i+1) \| i IN [0..3]] + DNI (b[0] + b[1]/3*(T_o^3 - T_i^3));
64	(* NOTE above eq assumes linear temperature rise with position? is that suff valid? *)
65	(* NOTE above eq assumes ambient temperature of 25 C (Patnode sect 2.3.2) -- should adjust for changes in T_amb?? *)
66	(* NOTE also the equation in Lippke 1995 for bare tubes, which has a different form and requires wind speed *)
67	(* NOTE that Lippke says that end loss must be included in the above equation (eq 7, p. 11) *)
68	END absorber_loss;
69
70	MODEL absorber_loss_test;
71	type IS_A symbol_constant;
72	type :== 'vacuum';
73	T_i, T_o IS_A temperature;
74	DNI IS_A intensity;
75	abso IS_A absorber_loss(type,T_i,T_o,DNI);
76	Qd_loss ALIASES abso.Qd_loss;
77	METHODS
78	METHOD on_load;
79	FIX DNI, T_i, T_o;
80	DNI := 1000 {W/m^2};
81	T_i := 300 {K} + 273.15 {K};
82	T_o := 400 {K} + 237.15 {K};
83	END on_load;
84	END absorber_loss_test;
85
86
87	(*
88	Heat absorbed by the solar receiver of cylindrical shape
89
90	TODO implement option to set field orientation EW/NS.
91	*)
92	MODEL parabolic_trough;
93	(* field geometry *)
94	W_apert "collector aperture width" IS_A distance;
95	L_spacing "'length of spacing' between troughs" IS_A distance; (* clarification required, assume centre-to-centre? *)
96	L_SCA "length of a single solar collector assembly (SCA)" IS_A distance;
97	N_SCA "total number of solar collector assemblies" IS_A factor;
98	f "focal length of collector" IS_A distance;
99
100	A_field IS_A area;
101	A_field = N_SCA * L_SCA * W_apert;
102
103	(* tracker and sun position *)
104	sun IS_A solar_tracker_single_horiz;
105	t "local time (standard, not DST correction)" ALIASES sun.t;
106
107	Qdd_abs "solar radiation absorbed by receiver tubes per aperture area" IS_A power_per_area;
108	DNI "direct normal irradiance" IS_A intensity;
109	IAM "incidence angle modifier" IS_A factor;
110	eta_shading IS_A fraction;
111	eta_endloss IS_A fraction;
112	eta_field "averaged field efficiency" IS_A fraction;
113	eta_HCE "averaged heat collection element efficiency" IS_A fraction;
114	avail_SF "solar field availability - on-sun portion" IS_A fraction;
115
116	(* constants that define the location -- can we do this better?? *)
117	(* Kramer Junction *)
118	sun.loc.latitude :== 37.21 {deg};
119	sun.loc.longitude :== -117.022 {deg};
120	sun.loc.elevation :== 755 {ft};
121
122	(* OPTICAL PERFORMANCE *)
123
124	(* absorbed solar radiation, Patnode eq 2.1. factored the cos(theta) into IAM. *)
125	absheat_eq: Qdd_abs = DNI * IAM * eta_shading * eta_endloss * eta_field * eta_HCE * avail_SF;
126	(* NOTE that eta_field and eta_HCE are given fixed values here, though their
127	components factors are shown in Patnode sect 2.2.5 *)
128
129	(* TODO where is Patnode eq 2.8 ? looks like an error? *)
130
131	(* incidence angle modifier, Patnode eq 2.9 *)
132	(* TODO check... do these values assume angles in degrees? *)
133	IAM * cos(sun.theta) = cos(sun.theta) + 8.84e-4 * sun.theta - 5.369e-5 * sun.theta^2;
134
135	(* row shading (Stuetzle), Patnode eq 2.12 *)
136	eta_shading_1, eta_shading_2 IS_A factor; (* before being normalised to [0,1] *)
137	shading_1: eta_shading_1 = (L_spacing/W_apert) * (cos(sun.zenith)/cos(sun.theta));
138	(* limit first to max of 1, then min of 0 *)
139	shading_2: eta_shading_2 = 1. + 0.5*((eta_shading_1 - 1) - abs(eta_shading_1 - 1));
140	shading_3: eta_shading = (eta_shading_2 + abs(eta_shading_2))/2;
141
142	(* end loss (Lippke), Patnode eq 2.13 *)
143	eta_endloss_1, eta_endloss_2 IS_A factor;
144	endloss_1: eta_endloss_1 = 1 - f * tan(sun.theta) / L_SCA;
145	(* limit first to max of 1, then min of 0 *)
146	endloss_2: eta_endloss_2 = 1. + 0.5*((eta_endloss_1 - 1) - abs(eta_endloss_1 - 1));
147	endloss_3: eta_endloss = (eta_endloss_2 + abs(eta_endloss_2))/2;
148
149	(* LOSSES FROM HCE *)
150
151	T_i, T_o IS_A temperature;
152	hce_types IS_A set OF symbol_constant;
153	hce_types :== ['air','vacuum','hydrogen'];
154
155	abso[hce_types] IS_A absorber_loss_base;
156	hce_frac[hce_types] IS_A fraction;
157	FOR i IN hce_types CREATE
158	abso[i] IS_REFINED_TO absorber_loss(i, T_i, T_o, DNI);
159	END FOR;
160
161	Qdd_loss_recv "HCE losses per aperture area" IS_A power_per_area;
162	(* receiver heat loss, Patnode eq 2.19 *)
163	Qdd_loss_recv * W_apert = SUM[abso[i].Qd_loss * hce_frac[i] \| i IN hce_types];
164
165	(* LOSSES FROM FIELD PIPING *)
166
167	T_amb "ambient temperature" IS_A temperature;
168	DT "temperature difference ambient to field piping" IS_A delta_temperature;
169	(* average external temp difference, Patnode eq 2.21 *)
170	T_amb + DT = 0.5 * (T_i + T_o);
171
172	Qdd_loss_pipe "losses from pipework, per solar field aperture area" IS_A power_per_area;
173
174	(* solar field piping loses, Patnode eq 2.20 *)
175	Qdd_loss_pipe = 0.01693{W/m^2/K}DT - 0.0001683{W/m^2/K^2}(DT^2) + 0.78e-7{W/m^2/K^3}*(DT^3);
176	(* typically 10 W/m2 or less, apparently *)
177
178	(* NET ENERGY GAIN IN HTF *)
179
180	Vdot_i IS_A volume_rate;
181	Qdd_fluid IS_A power_per_area;
182	(* net absorbed energy, Patnode eq 2.22 *)
183	Qdd_fluid = Qdd_abs - Qdd_loss_pipe - Qdd_loss_recv;
184
185	inlet, outlet IS_A therminol;
186	inlet.T, T_i ARE_THE_SAME;
187	outlet.T, T_o ARE_THE_SAME;
188
189	Q_net IS_A energy_rate;
190	Q_net = Qdd_fluid * A_field;
191
192	(* first-law energy balance, Patnode eq 2.23 *)
193	(outlet.h - inlet.h) * inlet.rho * Vdot_i = Q_net;
194
195	METHODS
196	METHOD specify;
197	RUN sun.specify;
198	FIX DNI, L_spacing;
199	FIX W_apert, f, L_SCA, eta_field, eta_HCE, avail_SF;
200
201	FIX inlet.T;
202	FIX T_amb;
203	FIX N_SCA;
204	FIX Vdot_i;
205	FIX hce_frac[hce_types];
206	END specify;
207
208	METHOD values;
209	L_spacing := 13 {m};
210	W_apert := 4.83 {m};
211	f := 5 {m};
212	L_SCA := 50 {m};
213	N_SCA := 256;
214
215	DNI := 1000 {W/m^2};
216
217	(* tracker, including sun position and geographical location *)
218	RUN sun.values;
219	sun.offset.tz := -8{h};
220	sun.offset.h := 0;
221	t := 43.5 {d};
222
223	(* fixed field efficiency based on values in Patnode Table 2.1 *)
224	eta_field := 0.857;
225	eta_HCE := 0.832;
226	avail_SF := 1;
227
228	inlet.T := 60{K} + 273.15 {K};
229	T_amb := 30{K} + 273.15 {K};
230	Vdot_i := 400{m^3/hour};
231
232	hce_frac['air'] := 1.0;
233	hce_frac['vacuum'] := 0.0;
234	hce_frac['hydrogen'] := 0.0;
235
236	(* initial guess, needed for solving OK *)
237	T_o := 200 {K} + 273.15{K};
238	END values;
239
240	METHOD on_load;
241	RUN specify;
242	RUN values;
243	END on_load;
244	END parabolic_trough;