/[ascend]/trunk/models/steam/dsgsat3.a4c
ViewVC logotype

Contents of /trunk/models/steam/dsgsat3.a4c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1312 - (show annotations) (download) (as text)
Sun Mar 4 14:04:39 2007 UTC (15 years, 4 months ago) by johnpye
File MIME type: text/x-ascend
File size: 8455 byte(s)
Fixed the model again :-(
1 REQUIRE "ivpsystem.a4l";
2 REQUIRE "atoms.a4l";
3 REQUIRE "johnpye/thermo_types.a4c";
4
5 (*
6 An attempt to model direct steam generation in pipe flow, limited to the
7 saturated regime, and with constant-valued friction factor. External heat
8 loss is also simplified.
9 *)
10 REQUIRE "steam/satsteamstream.a4c";
11
12 MODEL dsgsat3;
13 n IS_A integer_constant;
14 n :== 7;(* with L = 10m: 5,6,7,8,9,10,11 *)
15 (* with L = 5m: 2,3,4,5,7,9,11,12,13,1415,16 *)
16
17 dz IS_A real_constant;
18 L IS_A real_constant;
19 L :== 10 {m};
20 dz :== L / (n-1);
21
22 nodes,butfirst1,upwind4,central IS_A set OF integer_constant;
23 nodes :== [1..n];
24 butfirst1 :== nodes - [1];
25 upwind4 :== nodes - [1,2,n];
26 central :== nodes - [1,n];
27
28 (* temporal derivatives *)
29 drho_dt[butfirst1] IS_A density_rate;
30 dmdot_dt[butfirst1] IS_A mass_rate_rate;
31 du_dt[butfirst1] IS_A specific_energy_rate;
32 dTw_dt[butfirst1] IS_A temperature_rate;
33
34 (* wall properties *)
35 rho_w IS_A mass_density;
36 D, D_2 IS_A distance;
37 c_w IS_A specific_heat_capacity;
38 A, A_w IS_A area;
39 h_int IS_A heat_transfer_coefficient; (* internal *)
40 h_ext IS_A heat_transfer_coefficient; (* external *)
41 z_A: A = 1{PI}*D^2/4;
42 z_Aw: A_w = 1{PI}*(D_2^2 - D^2)/4;
43
44 (* fluid properties *)
45 node[nodes] IS_A satsteamstream;
46
47 (* flow properties *)
48 vel[nodes] IS_A speed;
49 T_w[butfirst1] IS_A temperature;
50 T[nodes] IS_A temperature;
51
52 (* constant, for the moment: *)
53 f IS_A positive_factor;
54 (* mu_f IS_A viscosity; *)
55 T_amb IS_A temperature;
56
57 (* system dynamics *)
58 qdot_t[butfirst1], qdot_l[butfirst1] IS_A power_per_length;
59 qdot_s IS_A power_per_length;
60
61 FOR i IN nodes CREATE
62 z_vel[i]: vel[i] = v[i]*mdot[i]/A;
63 END FOR;
64
65 (* some aliases just for easier review of the state of the model *)
66 x[nodes] IS_A fraction;
67 mdot[nodes] IS_A mass_rate;
68 p[nodes] IS_A pressure;
69 rho[nodes] IS_A mass_density;
70 u[nodes] IS_A specific_energy;
71 v[nodes] IS_A specific_volume;
72 FOR i IN nodes CREATE
73 x[i], node[i].x ARE_THE_SAME;
74 mdot[i], node[i].mdot ARE_THE_SAME;
75 p[i], node[i].p ARE_THE_SAME;
76 T[i], node[i].T ARE_THE_SAME;
77 rho[i], node[i].rho ARE_THE_SAME;
78 u[i], node[i].u ARE_THE_SAME;
79 v[i], node[i].v ARE_THE_SAME;
80 END FOR;
81
82 en_upwind4,en_central,mom_upwind4,mom_central,mass_upwind4,mass_central IS_A set OF integer_constant;
83
84 (* mass conservation *)
85 mass_upwind4 :== upwind4;
86 mass_central :== [];
87 FOR i IN mass_upwind4 CREATE (* 4-pt upwind biased *)
88 z_massbal1[i]: A * drho_dt[i] * dz =
89 - (mdot[i+1] + 6.*mdot[i] - 3.*mdot[i-1] - 2.*mdot[i-2]) / 6.;
90 END FOR;
91 FOR i IN mass_central CREATE
92 z_massbal2[i]: A * drho_dt[i] * dz =
93 - (mdot[i+1] - mdot[i-1]) / 2.;
94 END FOR;
95 FOR i IN butfirst1 - mass_upwind4 - mass_central CREATE
96 z_massbal[i]: A * drho_dt[i] * dz = - (mdot[i] - mdot[i-1]);
97 END FOR;
98
99 (* energy conservation *)
100 en_upwind4 :== [];
101 en_central :== central;
102 FOR i IN en_upwind4 CREATE
103 z_enbal2[i]: dz * (qdot_t[i] - rho[i] * A * du_dt[i]) =
104 + mdot[i] * (node[i+1].u + 6.*u[i] - 3.*u[i-1] - 2.*u[i-1]) / 6.
105 + (p[i+1]*node[i+1].v*mdot[i+1] - p[i-1]*v[i-1]*mdot[i-1]) / 2.;
106 END FOR;
107 FOR i IN en_central CREATE
108 z_enbal1[i]: dz * (qdot_t[i] - rho[i] * A * du_dt[i]) =
109 + mdot[i] * (u[i] - u[i-1]) (* NOTE: not central *)
110 + (p[i+1]*v[i+1]*mdot[i+1] - p[i-1]*v[i-1]*mdot[i-1]) / 2.;
111 END FOR;
112 FOR i IN butfirst1 - en_upwind4 - en_central CREATE
113 z_enbal[i]: dz * (qdot_t[i] - rho[i] * A * du_dt[i]) =
114 + mdot[i] * (u[i] - u[i-1])
115 + (p[i]*v[i]*mdot[i] - p[i-1]*v[i-1]*mdot[i-1]);
116 END FOR;
117
118 (* momentum conservation *)
119 mom_upwind4 :== [];
120 mom_central :== central;
121 FOR i IN mom_upwind4 CREATE
122 z_mombal2[i]: - dz/A * dmdot_dt[i]
123 = (p[i]-p[i-1]) (* backdiff for pressure *)
124 + dz * f/D/2 * rho[i] * vel[i]^2
125 + (rho[i+1]*vel[i+1]^2 + 6.*rho[i]*vel[i]^2 - 3.*rho[i-1]*vel[i-1]^2 - 2.*rho[i-2]*vel[i-2]^2) / 6.;
126 END FOR;
127 FOR i IN mom_central CREATE
128 z_mombal1[i]: - dz/A * dmdot_dt[i]
129 = (p[i+1]-p[i-1]) / 2.
130 + dz * f/D/2 * rho[i] * vel[i]^2
131 + (rho[i+1]*vel[i+1]^2 - rho[i-1]*vel[i-1]^2) / 2.;
132 END FOR;
133 FOR i IN butfirst1 - mom_upwind4 - mom_central CREATE
134 z_mombal[i]: - dz/A * dmdot_dt[i]
135 = (p[i]-p[i-1])
136 + dz * f/D/2 * rho[i] * vel[i]^2
137 + (rho[i]*vel[i]^2 - rho[i-1]*vel[i-1]^2);
138 END FOR;
139
140 (* internal/external convection, and thermal mass of wall -- no spatial derivs here *)
141 FOR i IN butfirst1 CREATE
142 z_wall[i]: rho_w*A_w*c_w*dTw_dt[i] = qdot_s - qdot_l[i] - qdot_t[i];
143 z_loss[i]: qdot_l[i] = h_ext*(1{PI}*D_2)*(T_w[i] - T_amb);
144 z_trans[i]: qdot_t[i] = h_int*(1{PI}*D) *(T_w[i] - T[i]);
145 END FOR;
146
147 t IS_A time;
148 METHODS
149 METHOD bound_self;
150 vel[nodes].upper_bound := 100 {m/s};
151 qdot_l[butfirst1].lower_bound := 0 {W/m};
152 FOR i IN nodes DO
153 RUN node[i].bound_self;
154 END FOR;
155 END bound_self;
156 METHOD default;
157 (* these are initial guesses only; fixed parameters are overwritten by 'values' below *)
158 D := 0.06 {m};
159 D_2 := 0.07 {m};
160 A_w := 0.25{PI}*D_2^2 -0.25{PI}*D^2;
161 A := 1 {m^2};
162 T_amb := 298 {K};
163 f := 0.0;
164 h_ext := 10 {W/m^2/K};
165 h_int := 100 {W/m^2/K};
166 qdot_s := 1000 {W/m};
167 rho_w := 1000 {kg/m^3};
168 t := 0 {s};
169 FOR i IN nodes DO
170 T[i] := 298 {K};
171 vel[i] := 1 {m/s};
172 RUN node[i].default_self;
173 END FOR;
174 FOR i IN butfirst1 DO
175 T_w[i] := 298 {K};
176 drho_dt[i] := 0 {kg/m^3/s};
177 dmdot_dt[i] := 0 {kg/s/s};
178 du_dt[i] := 0 {kJ/kg/s};
179 dTw_dt[i] := 0 {K/s};
180 qdot_t[i] := 0 {W/m};
181 qdot_l[i] := 0 {W/m};
182 x[i] := x[1];
183 END FOR;
184 END default;
185 METHOD specify;
186 (* change to a proper steady-state problem, with fluid properties FREEd *)
187 FOR i IN nodes DO
188 RUN node[i].specify;
189 FIX dTw_dt[i]; FREE T_w[i];
190 END FOR;
191 FIX p[1];
192 FREE T[1];
193 FIX qdot_s;
194 FIX D, D_2;
195 FIX h_int, c_w, rho_w, h_ext;
196 FIX f;
197 (* FIX mu_f; *)
198 FIX T_amb;
199 (* fix derivatives to zero *)
200 FOR i IN butfirst1 DO
201 FREE x[i]; FIX p[i];
202 FIX drho_dt[i]; FREE p[i];
203 FIX du_dt[i]; FREE T[i];
204 FREE mdot[i]; FIX dmdot_dt[i];
205 END FOR;
206 END specify;
207 METHOD values;
208 h_int := 100 {W/m^2/K};
209 h_ext := 20 {W/m^2/K};
210 f := 0.03;
211 mdot[1] := 0.26 {kg/s};
212 p[1] := 10 {bar};
213 x[1] := 0.23;
214 qdot_s := 1000 {W/m^2} * D_2 * 10;
215 FOR i IN butfirst1 DO
216 dmdot_dt[i] := 0.0 {kg/s/s};
217 du_dt[i] := 0 {kJ/kg/s};
218 v[i] := 0.2 {L/kg};
219 rho[i] := 6 {kg/L};
220 node[i].dp_dT := +0.5 {kPa/K};
221 p[i] := 5 {bar};
222 END FOR;
223 END values;
224 METHOD on_load;
225 RUN configure_steady;
226 RUN ode_init;
227 END on_load;
228 (*---------------- a physically sensible steady-state configuration-----------*)
229 METHOD configure_steady;
230 RUN default_self;
231 RUN ClearAll;
232 RUN specify;
233 RUN bound_steady;
234 RUN values;
235 END configure_steady;
236 METHOD bound_steady;
237 RUN bound_self;
238 T_w[butfirst1].upper_bound := 1000 {K};
239 END bound_steady;
240 (*------------------------- the dynamic problem ------------------------------*)
241 METHOD configure_dynamic;
242 FOR i IN butfirst1 DO
243 FREE drho_dt[i]; FIX rho[i];
244 FREE dmdot_dt[i]; FIX mdot[i];
245 FREE du_dt[i]; FIX u[i];
246 FREE dTw_dt[i]; FIX T_w[i];
247 FREE x[i];
248 FREE T[i];
249 END FOR;
250 t := 0 {s};
251 END configure_dynamic;
252 METHOD free_states;
253 FOR i IN butfirst1 DO
254 FREE rho[i];
255 FREE mdot[i];
256 FREE u[i];
257 FREE T_w[i];
258 END FOR;
259 END free_states;
260 METHOD ode_init;
261 (* add the necessary meta data to allow solving with the integrator *)
262 t.ode_type := -1;
263
264 FOR i IN butfirst1 DO
265 drho_dt[i].ode_id := 4*i; rho[i].ode_id := 4*i;
266 drho_dt[i].ode_type := 2; rho[i].ode_type := 1;
267
268 dmdot_dt[i].ode_id := 4*i+1; mdot[i].ode_id := 4*i+1;
269 dmdot_dt[i].ode_type := 2; mdot[i].ode_type := 1;
270
271 du_dt[i].ode_id := 4*i+2; u[i].ode_id := 4*i+2;
272 du_dt[i].ode_type := 2; u[i].ode_type := 1;
273
274 dTw_dt[i].ode_id := 4*i+3; T_w[i].ode_id := 4*i+3;
275 dTw_dt[i].ode_type := 2; T_w[i].ode_type := 1;
276 END FOR;
277
278 FOR i IN nodes DO
279 (* p[i].obs_id := 1 + 10*i; *)
280 (* x[i].obs_id := 2 + 10*i; *)
281 (* mdot[i].obs_id := 4 + 10*i; *)
282 (* node[i].h.obs_id := 3 + 10*i; *)
283 END FOR;
284 FOR i IN butfirst1 DO
285 (* qdot_t[i].obs_id := 3 + 10*i; *)
286 (* T_w[i].obs_id := 5 + 10*i; *)
287 (* T[i].obs_id := 6 + 10*i;*)
288 END FOR;
289
290 mdot[n].obs_id := 1;
291 x[n].obs_id := 1;
292 p[n].obs_id := 1;
293 vel[n].obs_id := 1;
294
295 END ode_init;
296 METHOD fix_outlet_quality;
297 FIX x[n];
298 FREE mdot[1];
299 END fix_outlet_quality;
300
301 END dsgsat3;
302 ADD NOTES IN dsgsat2;
303 'QRSlv' iterationlimit {50}
304 END NOTES;

john.pye@anu.edu.au
ViewVC Help
Powered by ViewVC 1.1.22