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 ASCEND version by Vikram Kaadam (2011), re-written by 
John Pye (2012).
*)
REQUIRE "atoms.a4l";
REQUIRE "solar/solar_types.a4l";
REQUIRE "johnpye/thermo_types.a4c";
REQUIRE "johnpye/sunpos.a4c";
REQUIRE "solar/therminol.a4c";

(*
	Sun position calculator, using Duffie & Beckman equations.
	
	TODO change to use grena/PSA/NREL algorithm instead of D&B
*)
MODEL sunpos_wrapper REFINES sunpos;
		
	beta = 1.0 * theta_z;

METHODS
	METHOD specify;
		FIX L_st, L_loc, phi; (* time and location *)
		FIX gamma; (* surface orientation *)
		FIX t;			
	END specify;

	METHOD values;
		(* FIXME there should be something here?? *)
	END values;

	METHOD self_test;
		ASSERT abs(theta-35.0{deg}) < 0.15{deg};
		ASSERT abs(delta-(-13.80{deg})) < 0.02{deg};
	END self_test;		
END sunpos_wrapper;

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;
	
	(* sun position *)
	sun        IS_A sunpos_wrapper;
	(* TODO what do we do with this? *)
	DST        "daylight savings adjust (=1h during DST, else 0)" IS_A time;
	t_std 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;
	row_shadow IS_A factor;
	end_loss   IS_A factor;
	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;

	(* OPTICAL PERFORMANCE *)
	
	(* absorbed solar radiation, Patnode eq 2.1. factored the cos(theta) into IAM. *)	
	Qdd_abs = DNI * IAM * row_shadow * end_loss * 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 *)
	(* TODO is this the right theta_z and theta? *)
	row_shadow = (L_spacing/W_apert) * (cos(sun.theta_z)/cos(sun.theta));	
	
	(* end loss (Lippke), Patnode eq 2.13 *)
	end_loss = 1 - f * tan(sun.theta) / L_SCA;
	
	(* 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 := 15 {m};
		W_apert := 5 {m};
		f := 5 {m};
		L_SCA := 50 {m};
		N_SCA := 256;

		DNI := 1000 {W/m^2};
		sun.L_st := -105 {deg};
		sun.L_loc := -110 {deg};
		DST := 0 {hour};
		t_std := 15 {hour};
		sun.phi := 37.21 {deg};

		(* TODO what's this? *)
		IF (t_std > 12{hour}) THEN
			sun.gamma := 90{deg};
		END IF;
		IF (t_std <= 12{hour}) THEN
			sun.gamma := -90{deg};
		END IF;

		(* 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 ASCEND version by Vikram Kaadam (2011), re-written by
7	John Pye (2012).
8	*)
9	REQUIRE "atoms.a4l";
10	REQUIRE "solar/solar_types.a4l";
11	REQUIRE "johnpye/thermo_types.a4c";
12	REQUIRE "johnpye/sunpos.a4c";
13	REQUIRE "solar/therminol.a4c";
14
15	(*
16	Sun position calculator, using Duffie & Beckman equations.
17
18	TODO change to use grena/PSA/NREL algorithm instead of D&B
19	*)
20	MODEL sunpos_wrapper REFINES sunpos;
21
22	beta = 1.0 * theta_z;
23
24	METHODS
25	METHOD specify;
26	FIX L_st, L_loc, phi; (* time and location *)
27	FIX gamma; (* surface orientation *)
28	FIX t;
29	END specify;
30
31	METHOD values;
32	(* FIXME there should be something here?? *)
33	END values;
34
35	METHOD self_test;
36	ASSERT abs(theta-35.0{deg}) < 0.15{deg};
37	ASSERT abs(delta-(-13.80{deg})) < 0.02{deg};
38	END self_test;
39	END sunpos_wrapper;
40
41	MODEL absorber_loss_base;
42	END absorber_loss_base;
43
44	(*
45	Heat loss per metre length of absorber tube, including convective, conductive
46	and radiative losses.
47	*)
48	MODEL absorber_loss(
49	type IS_A symbol_constant;
50	T_i WILL_BE temperature;
51	T_o WILL_BE temperature;
52	DNI WILL_BE intensity;
53	) REFINES absorber_loss_base;
54
55	a[0..3] IS_A real_constant;
56	b[0,1] IS_A real_constant;
57
58	SELECT(type)
59	CASE 'vacuum':
60	(* an intact HCE, evacuated air at 0.0001 Torr *)
61	a[0] :== -2.247372E+01 {W/m^2};
62	a[1] :== 8.374490E-01 {W/m^2/K};
63	a[2] :== 0.00 {W/m^2/K^2};
64	a[3] :== 4.620143E-06 {W/m^2/K^3};
65	b[0] :== 6.983190E-02 {m};
66	b[1] :== 9.312703E-08 {m/K^2};
67	CASE 'hydrogen':
68	(* intact HCE into which hydrogen has diffused to a pressure off 1 Torr *)
69	a[0] :== -2.247372E+01 {W/m^2};
70	a[1] :== 8.374490E-01 {W/m^2/K};
71	a[2] :== 0.00 {W/m^2/K^2};
72	a[3] :== 4.620143E-06 {W/m^2/K^3};
73	b[0] :== 6.983190E-02 {m};
74	b[1] :== 9.312703E-08 {m/K^2};
75	CASE 'air':
76	(* a broken glass envelope; contents will be air at ambient pressure *)
77	a[0] :== -2.247372E+01 {W/m^2};
78	a[1] :== 8.374490E-01 {W/m^2/K};
79	a[2] :== 0.00 {W/m^2/K^2};
80	a[3] :== 4.620143E-06 {W/m^2/K^3};
81	b[0] :== 6.983190E-02 {m};
82	b[1] :== 9.312703E-08 {m/K^2};
83	END SELECT;
84
85	Qd_loss IS_A power_per_length;
86
87	(* heat loss calculated as integral wrt temperature, Patnode eq 2.18. *)
88	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));
89	(* NOTE above eq assumes linear temperature rise with position? is that suff valid? *)
90	(* NOTE above eq assumes ambient temperature of 25 C (Patnode sect 2.3.2) -- should adjust for changes in T_amb?? *)
91	(* NOTE also the equation in Lippke 1995 for bare tubes, which has a different form and requires wind speed *)
92	(* NOTE that Lippke says that end loss must be included in the above equation (eq 7, p. 11) *)
93	END absorber_loss;
94
95	MODEL absorber_loss_test;
96	type IS_A symbol_constant;
97	type :== 'vacuum';
98	T_i, T_o IS_A temperature;
99	DNI IS_A intensity;
100	abso IS_A absorber_loss(type,T_i,T_o,DNI);
101	Qd_loss ALIASES abso.Qd_loss;
102	METHODS
103	METHOD on_load;
104	FIX DNI, T_i, T_o;
105	DNI := 1000 {W/m^2};
106	T_i := 300 {K} + 273.15 {K};
107	T_o := 400 {K} + 237.15 {K};
108	END on_load;
109	END absorber_loss_test;
110
111
112	(*
113	Heat absorbed by the solar receiver of cylindrical shape
114
115	TODO implement option to set field orientation EW/NS.
116	*)
117	MODEL parabolic_trough;
118	(* field geometry *)
119	W_apert "collector aperture width" IS_A distance;
120	L_spacing "'length of spacing' between troughs" IS_A distance; (* clarification required, assume centre-to-centre? *)
121	L_SCA "length of a single solar collector assembly (SCA)" IS_A distance;
122	N_SCA "total number of solar collector assemblies" IS_A factor;
123	f "focal length of collector" IS_A distance;
124
125	A_field IS_A area;
126	A_field = N_SCA * L_SCA * W_apert;
127
128	(* sun position *)
129	sun IS_A sunpos_wrapper;
130	(* TODO what do we do with this? *)
131	DST "daylight savings adjust (=1h during DST, else 0)" IS_A time;
132	t_std ALIASES sun.t;
133
134	Qdd_abs "solar radiation absorbed by receiver tubes per aperture area" IS_A power_per_area;
135	DNI "direct normal irradiance" IS_A intensity;
136	IAM "incidence angle modifier" IS_A factor;
137	row_shadow IS_A factor;
138	end_loss IS_A factor;
139	eta_field "averaged field efficiency" IS_A fraction;
140	eta_HCE "averaged heat collection element efficiency" IS_A fraction;
141	avail_SF "solar field availability - on-sun portion" IS_A fraction;
142
143	(* OPTICAL PERFORMANCE *)
144
145	(* absorbed solar radiation, Patnode eq 2.1. factored the cos(theta) into IAM. *)
146	Qdd_abs = DNI * IAM * row_shadow * end_loss * eta_field * eta_HCE * avail_SF;
147	(* NOTE that eta_field and eta_HCE are given fixed values here, though their
148	components factors are shown in Patnode sect 2.2.5 *)
149
150	(* TODO where is Patnode eq 2.8 ? looks like an error? *)
151
152	(* incidence angle modifier, Patnode eq 2.9 *)
153	(* TODO check... do these values assume angles in degrees? *)
154	IAM * cos(sun.theta) = cos(sun.theta) + 8.84e-4 * sun.theta - 5.369e-5 * sun.theta^2;
155
156	(* row shading (Stuetzle), Patnode eq 2.12 *)
157	(* TODO is this the right theta_z and theta? *)
158	row_shadow = (L_spacing/W_apert) * (cos(sun.theta_z)/cos(sun.theta));
159
160	(* end loss (Lippke), Patnode eq 2.13 *)
161	end_loss = 1 - f * tan(sun.theta) / L_SCA;
162
163	(* LOSSES FROM HCE *)
164
165	T_i, T_o IS_A temperature;
166	hce_types IS_A set OF symbol_constant;
167	hce_types :== ['air','vacuum','hydrogen'];
168
169	abso[hce_types] IS_A absorber_loss_base;
170	hce_frac[hce_types] IS_A fraction;
171	FOR i IN hce_types CREATE
172	abso[i] IS_REFINED_TO absorber_loss(i, T_i, T_o, DNI);
173	END FOR;
174
175	Qdd_loss_recv "HCE losses per aperture area" IS_A power_per_area;
176	(* receiver heat loss, Patnode eq 2.19 *)
177	Qdd_loss_recv * W_apert = SUM[abso[i].Qd_loss * hce_frac[i] \| i IN hce_types];
178
179	(* LOSSES FROM FIELD PIPING *)
180
181	T_amb "ambient temperature" IS_A temperature;
182	DT "temperature difference ambient to field piping" IS_A delta_temperature;
183	(* average external temp difference, Patnode eq 2.21 *)
184	T_amb + DT = 0.5 * (T_i + T_o);
185
186	Qdd_loss_pipe "losses from pipework, per solar field aperture area" IS_A power_per_area;
187
188	(* solar field piping loses, Patnode eq 2.20 *)
189	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);
190	(* typically 10 W/m2 or less, apparently *)
191
192	(* NET ENERGY GAIN IN HTF *)
193
194	Vdot_i IS_A volume_rate;
195	Qdd_fluid IS_A power_per_area;
196	(* net absorbed energy, Patnode eq 2.22 *)
197	Qdd_fluid = Qdd_abs - Qdd_loss_pipe - Qdd_loss_recv;
198
199	inlet, outlet IS_A therminol;
200	inlet.T, T_i ARE_THE_SAME;
201	outlet.T, T_o ARE_THE_SAME;
202
203	Q_net IS_A energy_rate;
204	Q_net = Qdd_fluid * A_field;
205
206	(* first-law energy balance, Patnode eq 2.23 *)
207	(outlet.h - inlet.h) * inlet.rho * Vdot_i = Q_net;
208
209	METHODS
210	METHOD specify;
211	RUN sun.specify;
212
213	FIX DNI, L_spacing;
214	FIX W_apert, f, L_SCA, eta_field, eta_HCE, avail_SF;
215
216	FIX inlet.T;
217	FIX T_amb;
218	FIX N_SCA;
219	FIX Vdot_i;
220	FIX hce_frac[hce_types];
221	END specify;
222
223	METHOD values;
224	L_spacing := 15 {m};
225	W_apert := 5 {m};
226	f := 5 {m};
227	L_SCA := 50 {m};
228	N_SCA := 256;
229
230	DNI := 1000 {W/m^2};
231	sun.L_st := -105 {deg};
232	sun.L_loc := -110 {deg};
233	DST := 0 {hour};
234	t_std := 15 {hour};
235	sun.phi := 37.21 {deg};
236
237	(* TODO what's this? *)
238	IF (t_std > 12{hour}) THEN
239	sun.gamma := 90{deg};
240	END IF;
241	IF (t_std <= 12{hour}) THEN
242	sun.gamma := -90{deg};
243	END IF;
244
245	(* fixed field efficiency based on values in Patnode Table 2.1 *)
246	eta_field := 0.857;
247	eta_HCE := 0.832;
248	avail_SF := 1;
249
250	inlet.T := 60{K} + 273.15 {K};
251	T_amb := 30{K} + 273.15 {K};
252	Vdot_i := 400{m^3/hour};
253
254	hce_frac['air'] := 1.0;
255	hce_frac['vacuum'] := 0.0;
256	hce_frac['hydrogen'] := 0.0;
257
258	(* initial guess, needed for solving OK *)
259	T_o := 200 {K} + 273.15{K};
260	END values;
261
262	METHOD on_load;
263	RUN specify;
264	RUN values;
265	END on_load;
266	END parabolic_trough;