models/johnpye/cavity.a4c

(*
	ASCEND model library: radiation heat transfer model
	Copyright (C) John Pye 2006
	http://pye.dyndns.org/

	This program is free software; you can redistribute it
	and/or modify it under the terms of the GNU General Public
	License as published by the Free Software Foundation; either
	version 2 of the License, or (at your option) any later
	version.

	This program is distributed in the hope that it will be
	useful, but WITHOUT ANY WARRANTY; without even the implied
	warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
	PURPOSE.  See the GNU General Public License for more
	details.

	You should have received a copy of the GNU General Public
	License along with this program; if not, write to the Free
	Software Foundation, Inc., 59 Temple Place, Suite 330,
	Boston, MA 02111-1307  USA
	-------------------------------------------------------------------

	This is a radiation heat transfer model that is being used
	to simulate (just a rough model) the absorber of a large solar
	thermal energy system. For more information see 
	http://pye.dyndns.org or google for 'CLFR'. This is a linear
	concentrator, so the cavity is simulated with a two-dimensional 
	model. The cavity is trapezoidal in shape, wider at the base
	through which radiation enters.

	The upper inside surface 'W' is the absorber. The cavity has a flat
	base 'B', and side walls 'L' and 'R' (of length N). This model
	calculates view factors between these surfaces then implements a
	'network' model to determine the equilibrium temperatures and
	heat transfer between the surfaces inside the cavity. External heat
	losses are allowed for: radiation and convention at the cavity
	base, plus convection off the side-walls.
*)

REQUIRE "atoms.a4l";

ATOM energy_rate_per_length REFINES solver_var
	DIMENSION M*L/T^3
	DEFAULT 1000{W/m};
	lower_bound := -1e50{W/m};
	upper_bound := 1e50{W/m};
	nominal := 1000{W/m};
END energy_rate_per_length;

ATOM energy_flux REFINES solver_var
	DIMENSION M/T^3
	DEFAULT 1000{W/m^2};
	lower_bound := -1e50{W/m^2};
	upper_bound := 1e50{W/m^2};
	nominal := 1000{W/m^2};
END energy_flux;

ATOM heat_transfer_coefficient REFINES solver_var
		DIMENSION M/T^3/TMP
		DEFAULT 5{W/m^2/K};

	lower_bound := 0{W/m^2/K};
	upper_bound := 1e50{W/m^2/K};
	nominal := 5{W/m^2/K};
END heat_transfer_coefficient;

ATOM specific_heat_capacity REFINES solver_var
		DIMENSION L^2/T^2/TMP
		DEFAULT 4.0{kJ/kg/K};
	
	lower_bound := 0{kJ/kg/K};
	upper_bound := 1e50{kJ/kg/K};
	nominal := 4.2{kJ/kg/K};
END specific_heat_capacity;

MODEL cavity;
	W,B,S,N,C,E,D IS_A distance;
	theta, phi, psi IS_A angle;
	
	F_WN, F_WW, F_WB IS_A factor;
	F_NB, F_NW, F_NN, F_BW, F_BN IS_A factor;
	F_WS IS_A factor;
	
	z_N: N * sin(theta) = D;
	z_S: S * tan(theta) = D;
	z_B: B = W + 2 * S;
	z_E: E^2 = W^2 + D^2;
	z_phi: phi = arctan(D/W);
	z_psi: psi = 1{PI} - phi;
	z_C: C^2 = E^2 + S^2 - 2* E * S * cos(psi);
	
	z_F_WN: F_WN = (W+N-C)/2/W;
	z_F_WW: F_WW = (2*E-2*D)/2/W; (* from top to directly opp part of bottom *)
	z_F_WB: F_WB = 1 - 2 * F_WN;
	z_F_WS: F_WS = (1 - 2 * F_WN - F_WB)/2;
	z_F_NB: F_NB = (N+B-C)/2/N;
	z_F_NW: F_NW = (N+W-C)/2/N;
	
	z_F_BN: F_BN = F_NB*N/B;
	z_F_NN: F_NN = 1 - F_NW - F_NB;
	z_F_BW: F_BW = F_WB*W/B;
		
	n IS_A set OF symbol_constant;
	n :== ['W','B','L','R'];

	(* Put lengths into a vector *)
	
	A[n] IS_A distance;
	A['W'], W ARE_THE_SAME;
	A['B'], B ARE_THE_SAME;
	A['L'], N ARE_THE_SAME;
	A['R'], N ARE_THE_SAME;
	
	(* View factor matrix *)

	F[n][n] IS_A factor;
	
	F['W']['L'], F['W']['R'], F_WN ARE_THE_SAME;
	F['W']['B'], F_WB ARE_THE_SAME;
	z_F_WW1: F['W']['W'] = 0;
	
	F['B']['L'], F['B']['R'], F_BN ARE_THE_SAME;
	F['B']['W'], F_BW ARE_THE_SAME;
	z_F_BB1: F['B']['B'] = 0;
	
	F['L']['R'], F['R']['L'], F_NN ARE_THE_SAME;
	F['L']['B'], F['R']['B'], F_NB ARE_THE_SAME;
	F['L']['W'], F['R']['W'], F_NW ARE_THE_SAME;
	z_F_LL: F['L']['L'] = 0;
	z_F_RR: F['R']['R'] = 0;

	(* Radiation equations *)
	q[n] IS_A energy_rate_per_length;
	
	E_b[n] IS_A energy_flux;
	J[n] IS_A energy_flux;	
	T[n] IS_A temperature;
	eps[n] IS_A factor; (* emissivity *)
	
	FOR i IN n CREATE
		z_q[i]: q[i] =  SUM[ (J[i]-J[j])*(A[i]*F[i][j]) | j IN n];
	END FOR;
	
	FOR i IN n CREATE
		z_E_b[i]: E_b[i] = 1{SIGMA_C} * T[i]^4;
	END FOR;

	FOR i IN n CREATE
		z_J[i]: q[i] * (1-eps[i]) = (E_b[i] - J[i]) * (eps[i]*A[i]);
	END FOR;
		
METHODS
METHOD default_self;
	RUN reset; RUN values; RUN bound_self;
END default_self;

METHOD bound_self;
	phi.lower_bound := 0 {deg};
	phi.upper_bound := 90 {deg};
	psi.lower_bound := 90 {deg};
	psi.upper_bound := 180 {deg};
END bound_self;

METHOD specify;
	FIX T['W', 'B'];
	FIX q['L','R'];
	FIX W,D,theta;
	FIX eps[n];
END specify;

METHOD values;
	T['W'] := 550 {K};
	T['B'] := 373.15 {K};
	q['L','R'] := 0 {W/m};
	
	W := 500 {mm};
	D := 300 {mm};
	theta := 30 {deg};
	
	eps['W'] := 0.49;
	eps['B'] := 0.9;
	eps['L','R'] := 0.1;
	
END values;

END cavity;

(*
	The following model adds external convection coefficients to
	the model, and an ambient temperature.
	
	We also calculate the 'F_rad' correlation parameter (see
	http://pye.dyndns.org for more information).
*)
MODEL cavity_losses REFINES cavity;
	h_B, h_N IS_A heat_transfer_coefficient;
	T_amb IS_A temperature;
	
	(* external heat loss *)
	
	q_ext_rad IS_A energy_rate_per_length;
	q_ext_conv IS_A energy_rate_per_length;
	q_ext_conv = h_B * B * (T['B'] - T_amb);
	q_ext_rad = B * eps['B'] * 1{SIGMA_C} * (T['B']^4 - T_amb^4);
	z_q_B: - q['B'] = q_ext_conv + q_ext_rad;
	z_q_L: - q['L'] = h_N * N * (T['L'] - T_amb);
	z_q_R:- q['R'] = h_N * N * (T['L'] - T_amb);

	(* Determine 'F_rad' from overall heatloss... *)
	
	F_rad IS_A factor;
	F_rad_1 IS_A factor;
	
	z_F_rad: q['W'] = F_rad * W * eps['W'] * 1{SIGMA_C} * (T['W']^4 - T['B']^4);
	
	z_F_rad_1: q['W'] = F_rad_1 * W * 1{SIGMA_C} * (T['W']^4 - T['B']^4) /
		(1/eps['B'] + 1/eps['W'] - 1);

METHODS
METHOD specify;
	FIX T['W'], T_amb;
	FIX h_B, h_N;
	FIX W,D,theta;
	FIX eps[n];
END specify;
METHOD values;
	T['W'] := 550 {K};
	T_amb := 290 {K};
	W := 500 {mm};
	D := 300 {mm};
	theta := 30 {deg};
	eps['W'] := 0.49;
	eps['B'] := 0.9;
	eps['L','R'] := 0.1;
	h_B := 10 {W/m^2/K};
	h_N := 0.5 {W/m^2/K};
	(* free values *)
	T['L','R'] := 500 {K};
	T['B'] := 400 {K};
END values;
END cavity_losses;

(*
	This final model will examine the possible convection
	behaviour inside the stratified cavity
*)
MODEL cavity_convection;
	W, D_strat IS_A distance;
	k IS_A thermal_conductivity;
	Q IS_A energy_rate_per_length;
	T[1,2] IS_A temperature;
	
	Q * D_strat = W * k * (T[1] - T[2]); (* conduction through stratified zone*)

(*	rho IS_A mass_density;
	cp IS_A specific_heat_capacity;
	
	Pr IS_A factor;
	Pr = mu / rho / alpha;
	
	mu IS_A viscosity;
	alpha IS_A diffusivity;
	alpha = k / rho / cp;
*)
	
METHODS
METHOD specify;
	FIX k, mu, cp, rho;
	FIX T[1,2];
	FIX W, D_strat;
END specify;
METHOD default_self;
	RUN reset; RUN values;
END default_self;
METHOD values;
	D_strat := 0.15 {m};
	W := 0.5 {m};
	k := 0.039 {W/m/K};
	T[1] := 550 {K};
	T[2] := 332 {K};
END values;
END cavity_convection;

	
1	(*
2	ASCEND model library: radiation heat transfer model
3	Copyright (C) John Pye 2006
4	http://pye.dyndns.org/
5
6	This program is free software; you can redistribute it
7	and/or modify it under the terms of the GNU General Public
8	License as published by the Free Software Foundation; either
9	version 2 of the License, or (at your option) any later
10	version.
11
12	This program is distributed in the hope that it will be
13	useful, but WITHOUT ANY WARRANTY; without even the implied
14	warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
15	PURPOSE. See the GNU General Public License for more
16	details.
17
18	You should have received a copy of the GNU General Public
19	License along with this program; if not, write to the Free
20	Software Foundation, Inc., 59 Temple Place, Suite 330,
21	Boston, MA 02111-1307 USA
22	-------------------------------------------------------------------
23
24	This is a radiation heat transfer model that is being used
25	to simulate (just a rough model) the absorber of a large solar
26	thermal energy system. For more information see
27	http://pye.dyndns.org or google for 'CLFR'. This is a linear
28	concentrator, so the cavity is simulated with a two-dimensional
29	model. The cavity is trapezoidal in shape, wider at the base
30	through which radiation enters.
31
32	The upper inside surface 'W' is the absorber. The cavity has a flat
33	base 'B', and side walls 'L' and 'R' (of length N). This model
34	calculates view factors between these surfaces then implements a
35	'network' model to determine the equilibrium temperatures and
36	heat transfer between the surfaces inside the cavity. External heat
37	losses are allowed for: radiation and convention at the cavity
38	base, plus convection off the side-walls.
39	*)
40
41	REQUIRE "atoms.a4l";
42
43	ATOM energy_rate_per_length REFINES solver_var
44	DIMENSION M*L/T^3
45	DEFAULT 1000{W/m};
46	lower_bound := -1e50{W/m};
47	upper_bound := 1e50{W/m};
48	nominal := 1000{W/m};
49	END energy_rate_per_length;
50
51	ATOM energy_flux REFINES solver_var
52	DIMENSION M/T^3
53	DEFAULT 1000{W/m^2};
54	lower_bound := -1e50{W/m^2};
55	upper_bound := 1e50{W/m^2};
56	nominal := 1000{W/m^2};
57	END energy_flux;
58
59	ATOM heat_transfer_coefficient REFINES solver_var
60	DIMENSION M/T^3/TMP
61	DEFAULT 5{W/m^2/K};
62
63	lower_bound := 0{W/m^2/K};
64	upper_bound := 1e50{W/m^2/K};
65	nominal := 5{W/m^2/K};
66	END heat_transfer_coefficient;
67
68	ATOM specific_heat_capacity REFINES solver_var
69	DIMENSION L^2/T^2/TMP
70	DEFAULT 4.0{kJ/kg/K};
71
72	lower_bound := 0{kJ/kg/K};
73	upper_bound := 1e50{kJ/kg/K};
74	nominal := 4.2{kJ/kg/K};
75	END specific_heat_capacity;
76
77	MODEL cavity;
78	W,B,S,N,C,E,D IS_A distance;
79	theta, phi, psi IS_A angle;
80
81	F_WN, F_WW, F_WB IS_A factor;
82	F_NB, F_NW, F_NN, F_BW, F_BN IS_A factor;
83	F_WS IS_A factor;
84
85	z_N: N * sin(theta) = D;
86	z_S: S * tan(theta) = D;
87	z_B: B = W + 2 * S;
88	z_E: E^2 = W^2 + D^2;
89	z_phi: phi = arctan(D/W);
90	z_psi: psi = 1{PI} - phi;
91	z_C: C^2 = E^2 + S^2 - 2* E * S * cos(psi);
92
93	z_F_WN: F_WN = (W+N-C)/2/W;
94	z_F_WW: F_WW = (2E-2D)/2/W; (* from top to directly opp part of bottom *)
95	z_F_WB: F_WB = 1 - 2 * F_WN;
96	z_F_WS: F_WS = (1 - 2 * F_WN - F_WB)/2;
97	z_F_NB: F_NB = (N+B-C)/2/N;
98	z_F_NW: F_NW = (N+W-C)/2/N;
99
100	z_F_BN: F_BN = F_NB*N/B;
101	z_F_NN: F_NN = 1 - F_NW - F_NB;
102	z_F_BW: F_BW = F_WB*W/B;
103
104	n IS_A set OF symbol_constant;
105	n :== ['W','B','L','R'];
106
107	(* Put lengths into a vector *)
108
109	A[n] IS_A distance;
110	A['W'], W ARE_THE_SAME;
111	A['B'], B ARE_THE_SAME;
112	A['L'], N ARE_THE_SAME;
113	A['R'], N ARE_THE_SAME;
114
115	(* View factor matrix *)
116
117	F[n][n] IS_A factor;
118
119	F['W']['L'], F['W']['R'], F_WN ARE_THE_SAME;
120	F['W']['B'], F_WB ARE_THE_SAME;
121	z_F_WW1: F['W']['W'] = 0;
122
123	F['B']['L'], F['B']['R'], F_BN ARE_THE_SAME;
124	F['B']['W'], F_BW ARE_THE_SAME;
125	z_F_BB1: F['B']['B'] = 0;
126
127	F['L']['R'], F['R']['L'], F_NN ARE_THE_SAME;
128	F['L']['B'], F['R']['B'], F_NB ARE_THE_SAME;
129	F['L']['W'], F['R']['W'], F_NW ARE_THE_SAME;
130	z_F_LL: F['L']['L'] = 0;
131	z_F_RR: F['R']['R'] = 0;
132
133	(* Radiation equations *)
134	q[n] IS_A energy_rate_per_length;
135
136	E_b[n] IS_A energy_flux;
137	J[n] IS_A energy_flux;
138	T[n] IS_A temperature;
139	eps[n] IS_A factor; (* emissivity *)
140
141	FOR i IN n CREATE
142	z_q[i]: q[i] = SUM[ (J[i]-J[j])(A[i]F[i][j]) \| j IN n];
143	END FOR;
144
145	FOR i IN n CREATE
146	z_E_b[i]: E_b[i] = 1{SIGMA_C} * T[i]^4;
147	END FOR;
148
149	FOR i IN n CREATE
150	z_J[i]: q[i] * (1-eps[i]) = (E_b[i] - J[i]) * (eps[i]*A[i]);
151	END FOR;
152
153	METHODS
154	METHOD default_self;
155	RUN reset; RUN values; RUN bound_self;
156	END default_self;
157
158	METHOD bound_self;
159	phi.lower_bound := 0 {deg};
160	phi.upper_bound := 90 {deg};
161	psi.lower_bound := 90 {deg};
162	psi.upper_bound := 180 {deg};
163	END bound_self;
164
165	METHOD specify;
166	FIX T['W', 'B'];
167	FIX q['L','R'];
168	FIX W,D,theta;
169	FIX eps[n];
170	END specify;
171
172	METHOD values;
173	T['W'] := 550 {K};
174	T['B'] := 373.15 {K};
175	q['L','R'] := 0 {W/m};
176
177	W := 500 {mm};
178	D := 300 {mm};
179	theta := 30 {deg};
180
181	eps['W'] := 0.49;
182	eps['B'] := 0.9;
183	eps['L','R'] := 0.1;
184
185	END values;
186
187	END cavity;
188
189	(*
190	The following model adds external convection coefficients to
191	the model, and an ambient temperature.
192
193	We also calculate the 'F_rad' correlation parameter (see
194	http://pye.dyndns.org for more information).
195	*)
196	MODEL cavity_losses REFINES cavity;
197	h_B, h_N IS_A heat_transfer_coefficient;
198	T_amb IS_A temperature;
199
200	(* external heat loss *)
201
202	q_ext_rad IS_A energy_rate_per_length;
203	q_ext_conv IS_A energy_rate_per_length;
204	q_ext_conv = h_B * B * (T['B'] - T_amb);
205	q_ext_rad = B * eps['B'] * 1{SIGMA_C} * (T['B']^4 - T_amb^4);
206	z_q_B: - q['B'] = q_ext_conv + q_ext_rad;
207	z_q_L: - q['L'] = h_N * N * (T['L'] - T_amb);
208	z_q_R:- q['R'] = h_N * N * (T['L'] - T_amb);
209
210	(* Determine 'F_rad' from overall heatloss... *)
211
212	F_rad IS_A factor;
213	F_rad_1 IS_A factor;
214
215	z_F_rad: q['W'] = F_rad * W * eps['W'] * 1{SIGMA_C} * (T['W']^4 - T['B']^4);
216
217	z_F_rad_1: q['W'] = F_rad_1 * W * 1{SIGMA_C} * (T['W']^4 - T['B']^4) /
218	(1/eps['B'] + 1/eps['W'] - 1);
219
220	METHODS
221	METHOD specify;
222	FIX T['W'], T_amb;
223	FIX h_B, h_N;
224	FIX W,D,theta;
225	FIX eps[n];
226	END specify;
227	METHOD values;
228	T['W'] := 550 {K};
229	T_amb := 290 {K};
230	W := 500 {mm};
231	D := 300 {mm};
232	theta := 30 {deg};
233	eps['W'] := 0.49;
234	eps['B'] := 0.9;
235	eps['L','R'] := 0.1;
236	h_B := 10 {W/m^2/K};
237	h_N := 0.5 {W/m^2/K};
238	(* free values *)
239	T['L','R'] := 500 {K};
240	T['B'] := 400 {K};
241	END values;
242	END cavity_losses;
243
244	(*
245	This final model will examine the possible convection
246	behaviour inside the stratified cavity
247	*)
248	MODEL cavity_convection;
249	W, D_strat IS_A distance;
250	k IS_A thermal_conductivity;
251	Q IS_A energy_rate_per_length;
252	T[1,2] IS_A temperature;
253
254	Q * D_strat = W * k * (T[1] - T[2]); (* conduction through stratified zone*)
255
256	(* rho IS_A mass_density;
257	cp IS_A specific_heat_capacity;
258
259	Pr IS_A factor;
260	Pr = mu / rho / alpha;
261
262	mu IS_A viscosity;
263	alpha IS_A diffusivity;
264	alpha = k / rho / cp;
265	*)
266
267	METHODS
268	METHOD specify;
269	FIX k, mu, cp, rho;
270	FIX T[1,2];
271	FIX W, D_strat;
272	END specify;
273	METHOD default_self;
274	RUN reset; RUN values;
275	END default_self;
276	METHOD values;
277	D_strat := 0.15 {m};
278	W := 0.5 {m};
279	k := 0.039 {W/m/K};
280	T[1] := 550 {K};
281	T[2] := 332 {K};
282	END values;
283	END cavity_convection;
284
285