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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

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