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 *)
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	END absorber_loss;
93
94	MODEL absorber_loss_test;
95	type IS_A symbol_constant;
96	type :== 'vacuum';
97	T_i, T_o IS_A temperature;
98	DNI IS_A intensity;
99	abso IS_A absorber_loss(type,T_i,T_o,DNI);
100	Qd_loss ALIASES abso.Qd_loss;
101	METHODS
102	METHOD on_load;
103	FIX DNI, T_i, T_o;
104	DNI := 1000 {W/m^2};
105	T_i := 300 {K} + 273.15 {K};
106	T_o := 400 {K} + 237.15 {K};
107	END on_load;
108	END absorber_loss_test;
109
110
111	(*
112	Heat absorbed by the solar receiver of cylindrical shape
113
114	TODO implement option to set field orientation EW/NS.
115	*)
116	MODEL parabolic_trough;
117	(* field geometry *)
118	W_apert "collector aperture width" IS_A distance;
119	L_spacing "'length of spacing' between troughs" IS_A distance; (* clarification required, assume centre-to-centre? *)
120	L_SCA "length of a single solar collector assembly (SCA)" IS_A distance;
121	N_SCA "total number of solar collector assemblies" IS_A factor;
122	f "focal length of collector" IS_A distance;
123
124	A_field IS_A area;
125	A_field = N_SCA * L_SCA * W_apert;
126
127	(* sun position *)
128	sun IS_A sunpos_wrapper;
129	(* TODO what do we do with this? *)
130	DST "daylight savings adjust (=1h during DST, else 0)" IS_A time;
131	t_std ALIASES sun.t;
132
133	Qdd_abs "solar radiation absorbed by receiver tubes per aperture area" IS_A power_per_area;
134	DNI "direct normal irradiance" IS_A intensity;
135	IAM "incidence angle modifier" IS_A factor;
136	row_shadow IS_A factor;
137	end_loss IS_A factor;
138	eta_field "averaged field efficiency" IS_A fraction;
139	eta_HCE "averaged heat collection element efficiency" IS_A fraction;
140	avail_SF "solar field availability - on-sun portion" IS_A fraction;
141
142	(* OPTICAL PERFORMANCE *)
143
144	(* absorbed solar radiation, Patnode eq 2.1. factored the cos(theta) into IAM. *)
145	Qdd_abs = DNI * IAM * row_shadow * end_loss * eta_field * eta_HCE * avail_SF;
146	(* NOTE that eta_field and eta_HCE are given fixed values here, though their
147	components factors are shown in Patnode sect 2.2.5 *)
148
149	(* TODO where is Patnode eq 2.8 ? looks like an error? *)
150
151	(* incidence angle modifier, Patnode eq 2.9 *)
152	(* TODO check... do these values assume angles in degrees? *)
153	IAM * cos(sun.theta) = cos(sun.theta) + 8.84e-4 * sun.theta - 5.369e-5 * sun.theta^2;
154
155	(* row shading (Stuetzle), Patnode eq 2.12 *)
156	(* TODO is this the right theta_z and theta? *)
157	row_shadow = (L_spacing/W_apert) * (cos(sun.theta_z)/cos(sun.theta));
158
159	(* end loss (Lippke), Patnode eq 2.13 *)
160	end_loss = 1 - f * tan(sun.theta) / L_SCA;
161
162	(* LOSSES FROM HCE *)
163
164	T_i, T_o IS_A temperature;
165	hce_types IS_A set OF symbol_constant;
166	hce_types :== ['air','vacuum','hydrogen'];
167
168	abso[hce_types] IS_A absorber_loss_base;
169	hce_frac[hce_types] IS_A fraction;
170	FOR i IN hce_types CREATE
171	abso[i] IS_REFINED_TO absorber_loss(i, T_i, T_o, DNI);
172	END FOR;
173
174	Qdd_loss_recv "HCE losses per aperture area" IS_A power_per_area;
175	(* receiver heat loss, Patnode eq 2.19 *)
176	Qdd_loss_recv * W_apert = SUM[abso[i].Qd_loss * hce_frac[i] \| i IN hce_types];
177
178	(* LOSSES FROM FIELD PIPING *)
179
180	T_amb "ambient temperature" IS_A temperature;
181	DT "temperature difference ambient to field piping" IS_A delta_temperature;
182	(* average external temp difference, Patnode eq 2.21 *)
183	T_amb + DT = 0.5 * (T_i + T_o);
184
185	Qdd_loss_pipe "losses from pipework, per solar field aperture area" IS_A power_per_area;
186
187	(* solar field piping loses, Patnode eq 2.20 *)
188	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);
189	(* typically 10 W/m2 or less, apparently *)
190
191	(* NET ENERGY GAIN IN HTF *)
192
193	Vdot_i IS_A volume_rate;
194	Qdd_fluid IS_A power_per_area;
195	(* net absorbed energy, Patnode eq 2.22 *)
196	Qdd_fluid = Qdd_abs - Qdd_loss_pipe - Qdd_loss_recv;
197
198	inlet, outlet IS_A therminol;
199	inlet.T, T_i ARE_THE_SAME;
200	outlet.T, T_o ARE_THE_SAME;
201
202	Q_net IS_A energy_rate;
203	Q_net = Qdd_fluid * A_field;
204
205	(* first-law energy balance, Patnode eq 2.23 *)
206	(outlet.h - inlet.h) * inlet.rho * Vdot_i = Q_net;
207
208	METHODS
209	METHOD specify;
210	RUN sun.specify;
211
212	FIX DNI, L_spacing;
213	FIX W_apert, f, L_SCA, eta_field, eta_HCE, avail_SF;
214
215	FIX inlet.T;
216	FIX T_amb;
217	FIX N_SCA;
218	FIX Vdot_i;
219	FIX hce_frac[hce_types];
220	END specify;
221
222	METHOD values;
223	L_spacing := 15 {m};
224	W_apert := 5 {m};
225	f := 5 {m};
226	L_SCA := 50 {m};
227	N_SCA := 256;
228
229	DNI := 1000 {W/m^2};
230	sun.L_st := -105 {deg};
231	sun.L_loc := -110 {deg};
232	DST := 0 {hour};
233	t_std := 15 {hour};
234	sun.phi := 37.21 {deg};
235
236	(* TODO what's this? *)
237	IF (t_std > 12{hour}) THEN
238	sun.gamma := 90{deg};
239	END IF;
240	IF (t_std <= 12{hour}) THEN
241	sun.gamma := -90{deg};
242	END IF;
243
244	(* fixed field efficiency based on values in Patnode Table 2.1 *)
245	eta_field := 0.857;
246	eta_HCE := 0.832;
247	avail_SF := 1;
248
249	inlet.T := 60{K} + 273.15 {K};
250	T_amb := 30{K} + 273.15 {K};
251	Vdot_i := 400{m^3/hour};
252
253	hce_frac['air'] := 1.0;
254	hce_frac['vacuum'] := 0.0;
255	hce_frac['hydrogen'] := 0.0;
256
257	(* initial guess, needed for solving OK *)
258	T_o := 200 {K} + 273.15{K};
259	END values;
260
261	METHOD on_load;
262	RUN specify;
263	RUN values;
264	END on_load;
265	END parabolic_trough;