/[ascend]/trunk/models/johnpye/compressible_flow.a4c
ViewVC logotype

Contents of /trunk/models/johnpye/compressible_flow.a4c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2466 - (show annotations) (download) (as text)
Thu May 5 04:55:30 2011 UTC (11 years, 2 months ago) by jpye
File MIME type: text/x-ascend
File size: 14044 byte(s)
testing comment on bug 511.
1 (* ASCEND modelling environment
2 Copyright (C) 2011 Carnegie Mellon University
3
4 This program is free software; you can redistribute it and/or modify
5 it under the terms of the GNU General Public License as published by
6 the Free Software Foundation; either version 2, or (at your option)
7 any later version.
8
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
13
14 You should have received a copy of the GNU General Public License
15 along with this program; if not, write to the Free Software
16 Foundation, Inc., 59 Temple Place - Suite 330,
17 Boston, MA 02111-1307, USA.
18 *)(*
19 Calculations of isentropic shocks in air flow, following approach given
20 in Potter & Wiggert, 3rd SI edition.
21
22 Intention is to be able to generalise these models for arbitrary fluids as
23 supported by FPROPS. But the present models include ideal gas and constant
24 specific heat assumptions.
25
26 Author: John Pye
27 *)
28
29 REQUIRE "johnpye/thermo_types.a4c";
30
31 MODEL fluid;
32 component IS_A symbol_constant;
33 (* TODO, generalise for different fluid property correlations, selectable *)
34 END fluid;
35
36 (* XXX how to get hold of the fluid property constants... do we need them?? *)
37
38 (*
39 Ideal gas flow node, with ideal gas law, Mach number and speed of sound
40 included. No allowance for enthalpy, entropy at this stage.
41 *)
42 MODEL ideal_gas_node;
43 cd IS_A fluid;
44 M "Mach number" IS_A positive_factor;
45 p IS_A pressure;
46 T IS_A temperature;
47 v IS_A specific_volume;
48 rho IS_A mass_density;
49 Vel IS_A speed;
50
51 MM IS_A molar_weight_constant;
52 R IS_A specific_gas_constant;
53 R :== 1{GAS_C}/MM;
54
55 rho*v = 1;
56 p*v = R * T;
57
58 k "adiabatic index" IS_A real_constant;
59
60 c IS_A speed;
61 c = sqrt(k*R*T);
62
63 M = Vel / c;
64 END ideal_gas_node;
65
66
67 (*
68 Add cross-sectional area to the above model; not needed for normal shock
69 calculations, but useful in nozzle calcs.
70 *)
71 MODEL ideal_gas_duct_node;
72 (* we incorporate state as submodel so we can ARE_THE_SAME non-duct nodes *)
73 state IS_A ideal_gas_node;
74 M ALIASES state.M;
75 p ALIASES state.p;
76 T ALIASES state.T;
77 v ALIASES state.v;
78 rho ALIASES state.rho;
79 Vel ALIASES state.Vel;
80 MM ALIASES state.MM;
81 R ALIASES state.R;
82 k ALIASES state.k;
83 c ALIASES state.c;
84
85 A IS_A area;
86 mdot IS_A mass_rate;
87 mdoteq: mdot = rho * A * Vel;
88 END ideal_gas_duct_node;
89
90
91 (* Specific cases for air...*)
92
93 MODEL air_node REFINES ideal_gas_node;
94 MM :== 28.97 {kg/kmol};
95 k :== 1.4;
96 END air_node;
97
98 MODEL air_duct_node REFINES ideal_gas_duct_node;
99 MM :== 28.97 {kg/kmol};
100 k :== 1.4;
101 END air_duct_node;
102
103 (*----------------------------STAGNATION CONDITIONS---------------------------*)
104
105
106 MODEL isentropic_stagnation;
107 flow IS_A air_duct_node;
108 stag IS_A air_node;
109
110 (* conservation of mass *)
111 novel: stag.Vel = 0;
112 mdot ALIASES flow.mdot;
113 k ALIASES flow.k;
114
115 (* conservation of energy *)
116 conen: 1 = T_rel_0 * (1 + (k-1)/2 * flow.M^2);
117
118 (* isentropic *)
119 isen: (1 + (k-1)/2*flow.M^2)^(k/(k-1)) * p_rel_0 = 1;
120 (*(1 + (k-1)/2 * flow.M^2 ) * flow.p^((k-1)/k) = stag.p^((k-1)/k);*)
121 (*(stag.rho / flow.rho)^(k-1) = 1 + (k-1)/2*flow.M^2;*)
122
123 p_rel_0 IS_A positive_factor;
124 T_rel_0 IS_A positive_factor;
125 A_rel_star IS_A positive_factor;
126 A_star IS_A area;
127 p_0 ALIASES stag.p;
128 p ALIASES flow.p;
129 T ALIASES flow.T;
130 T_0 ALIASES stag.T;
131 A ALIASES flow.A;
132 M ALIASES flow.M;
133
134 preleq: p_rel_0 * stag.p = flow.p;
135 Treleq: T_rel_0 * stag.T = flow.T;
136 Areleq: A_rel_star * A_star = flow.A;
137
138 arat: A_rel_star * flow.M = ((2 + (k-1)*M^2)/(k+1))^((k+1)/(2*(k-1)));
139
140 END isentropic_stagnation;
141
142 MODEL isentropic_stagnation_test REFINES isentropic_stagnation;
143 METHODS
144 METHOD on_load;
145 FIX stag.p, stag.T, flow.A;
146 stag.p := 500 {kPa};
147 stag.T := 20 {K} + 273.15 {K};
148 flow.A := 10 {cm^2};
149
150 FIX M;
151 M := 0.80;
152 END on_load;
153 END isentropic_stagnation_test;
154
155 (*------------------------------EXAMPLE MODELS -------------------------------*)
156
157 (*
158 Flow from reservoir at 20°C, 500 kPa to receiver at 200 and 300 kPa.
159 Potter & Wiggert 3rd SI Edition, Example 9.2.
160
161 We assume the flow is through a converging-only nozzle from the reservoir
162 into the receiver. As such, either the entire flow is subsonic, or else
163 there will be choked flow with M = 1 in the exit of the nozzle.
164
165 'crit' uses upstream conditions to determine the reservoir pressure that
166 will give M = 1 at the exit. Then, in part (a), the receiver pressure is
167 higher than that critical pressure, so subsonic flow occurs. For (b), the
168 pressure is lower, so choked flow occurs, with flow rate limited by the
169 M_exit = 1 choked-flow limit.
170
171 Tested, works OK -- JP
172 *)
173 MODEL example_potter_wiggert_ex_9_2;
174 crit, part_a, part_b IS_A isentropic_stagnation;
175 crit.A, part_a.A, part_b.A ARE_THE_SAME;
176 crit.T_0, part_a.T_0, part_b.T_0 ARE_THE_SAME;
177 A ALIASES crit.A;
178 T_0 ALIASES crit.T_0;
179 METHODS
180 METHOD on_load;
181 FIX A; A := 10 {cm^2};
182 FIX T_0; T_0 := 20 {K} + 273.15 {K};
183
184 FIX crit.M; crit.M := 1;
185 FIX crit.p_0; crit.p_0 := 500 {kPa};
186
187 FIX part_a.p_0; part_a.p_0 := 500 {kPa};
188 FIX part_a.p; part_a.p := 300 {kPa};
189
190 (* actually, part (b) ends up being identical to 'crit'! *)
191 FIX part_b.p_0; part_b.p_0 := 500 {kPa};
192 FIX part_b.M; part_b.M := 1;
193 END on_load;
194 END example_potter_wiggert_ex_9_2;
195
196
197 (*
198 Potter & Wiggert 3rd SI ed, Example 9.3.
199 Converging-diverging nozzle with exit 40 cm² and throat 10 cm², with flow
200 coming from reservoir at 20 °C, 500 kPa. Two different exit pressures can
201 result in M = 1 at the throat -- we calculate what those pressure are.
202
203 For this problem, saturation states are the upstream reservoir, and we
204 assert that A_star is the throat area and A is the exit area. These areas
205 then yield the Mach number from eq 9.3.19 (iteratively). We find the two
206 solutions by having two instances of the 'isentropic_stagnation' model, and
207 impose upper and lower bounds on M in each of the two models.
208
209 Tested, works OK -- JP
210 *)
211 MODEL example_potter_wiggert_ex_9_3;
212 sub, sup IS_A isentropic_stagnation;
213 sub.A_star, sup.A_star ARE_THE_SAME;
214 A_star ALIASES sub.A_star;
215 sub.A, sup.A ARE_THE_SAME;
216 A ALIASES sub.A;
217 sub.p_0, sup.p_0 ARE_THE_SAME;
218 p_0 ALIASES sub.p_0;
219 sub.T_0, sup.T_0 ARE_THE_SAME;
220 T_0 ALIASES sub.T_0;
221 METHODS
222 METHOD on_load;
223 FIX A, A_star;
224 A_star := 10 {cm^2};
225 A := 40 {cm^2};
226
227 (* ensure two different solution regions *)
228 sub.M.upper_bound := 1;
229 sup.M.lower_bound := 1;
230
231 FIX T_0, p_0;
232 T_0 := 20 {K} + 273.15 {K};
233 p_0 := 500 {kPa};
234 END on_load;
235 END example_potter_wiggert_ex_9_3;
236
237 (*-----------------------NORMAL SHOCK IN ISENTROPIC FLOW----------------------*)
238
239 (*
240 Model of a stationary normal shock in air. Using a moving frame of
241 reference, it can be used to model a moving shock wave as well.
242
243 Equations from Potter & Wiggert, 3rd SI edition, sects 9.1, 9.2 and 9.4.
244 *)
245 MODEL air_normal_shock;
246 S1, S2 IS_A air_node;
247 k ALIASES S1.k;
248
249 S2.M^2 = ( S1.M^2 + 2/(k-1) ) / (2*k/(k-1)*S1.M^2 - 1);
250
251 S2.p / S1.p = 2*k/(k+1)*S1.M^2 - (k-1)/(k+1);
252
253 S2.T / S1.T = ( 1 + (k-1)/2*S1.M^2) * (2*k/(k-1)*S1.M^2 - 1)/ ((k+1)^2/2/(k-1)*S1.M^2);
254 END air_normal_shock;
255
256 MODEL air_duct_normal_shock;
257 S1, S2 IS_A air_duct_node;
258 shock IS_A air_normal_shock;
259 S1.state, shock.S1 ARE_THE_SAME;
260 S2.state, shock.S2 ARE_THE_SAME;
261 S1.A, S2.A ARE_THE_SAME;
262 END air_duct_normal_shock;
263
264 (*------------------------------EXAMPLE MODELS -------------------------------*)
265
266 (*
267 This model reproduces the results of Example 9.5 from Potter & Wiggert,
268 3rd SI edition. The question asks to determine the pressure and temperature
269 conditions downstream of a shock wave passing through ambient air of given
270 state.
271
272 Tested, works OK -- JP
273 *)
274 MODEL example_potter_wiggert_ex_9_5 REFINES air_normal_shock;
275 METHODS
276 METHOD on_load;
277 FIX S1.Vel, S1.p, S1.T;
278 S1.Vel := 450 {m/s};
279 S1.p := 80 {kPa};
280 S1.T := 15 {K} + 273.15 {K};
281 END on_load;
282 END example_potter_wiggert_ex_9_5;
283
284 (*
285 This model reproduces the results of Example 9.6 from Potter & Wiggert,
286 3rd SI edition. This problem shows the wind speeds implicit behind a strong
287 shock wave such as that arising from a high-powered bomb explosions.
288
289 Although the problem as given in P&W can be solved even without doing so,
290 we have added an assumption that the ambient air pressure is 101.3 kPa. This
291 allows the model to be 'square' and the pressure and temperature behind the
292 shock wave (4.83 bar, 500 K) to also be calculated.
293
294 Tested, works OK -- JP
295 *)
296 MODEL example_potter_wiggert_ex_9_6 REFINES air_normal_shock;
297 Vel_induced IS_A speed;
298 Vel_induced = S2.Vel - S1.Vel;
299 METHODS
300 METHOD on_load;
301 FIX S1.Vel, S1.T, S1.p;
302 S1.Vel := 700 {m/s};
303 S1.T := 15 {K} + 273.15 {K};
304 S1.p := 101.3 {kPa};
305 END on_load;
306 END example_potter_wiggert_ex_9_6;
307
308
309 (*
310 This model reproduces the results of Example 9.7 from Potter & Wiggert,
311 3rd SI edition. The problem relates to the calculation of outlet conditions
312 and flow rate for a converging-diverging nozzle, given specified upstream
313 reservoir pressure and temperature, given the throat and exit areas, and
314 subject to the fact that there is a normal shock at the nozzle's exit plane.
315
316 For a shock to exist at the exit, we know that we must have M = 1 at the
317 throat, so A* = A_t.
318
319 Tested, works OK -- JP
320 *)
321 MODEL example_potter_wiggert_ex_9_7;
322 noz IS_A isentropic_stagnation;
323 A_star ALIASES noz.A_star;
324 A ALIASES noz.A;
325
326 D_star, D IS_A distance;
327 0.25{PI}*D_star^2 = A_star;
328 0.25{PI}*D^2 = A;
329
330 p_0 ALIASES noz.p_0;
331 T_0 ALIASES noz.T_0;
332
333 shock IS_A air_duct_normal_shock;
334 shock.S1, noz.flow ARE_THE_SAME;
335
336 rec IS_A isentropic_stagnation;
337 rec.flow, shock.S2 ARE_THE_SAME;
338
339 p_2 ALIASES rec.p;
340 mdot ALIASES noz.mdot;
341
342 METHODS
343 METHOD on_load;
344 (* nozzle geometry *)
345 FIX D_star, D;
346 D_star := 5 {cm};
347 D := 10 {cm};
348
349 (* note that isentropic flow equations have to be bounded over a reasonable
350 range for the Brent algorithm to converge in this case. Note that P&W Table
351 D.1 gives M up to 10 only, so limiting to 100 is still more generous.
352 FIXME Perhaps we can make these models more user-friendly by having a
353 supersonic and subsonic model that sets these bounds automatically...? *)
354 noz.flow.M.lower_bound := 1;
355 noz.flow.M.upper_bound := 100;
356
357 (* reservoir conditions *)
358 FIX p_0, T_0;
359 p_0 := 90 {kPa};
360 T_0 := 20 {K} + 273.15 {K};
361 END on_load;
362 METHOD self_test;
363 (* values from P&W *)
364 ASSERT abs(p_2 - 26.6 {kPa}) < 0.05 {kPa};
365 ASSERT abs(mdot - 0.417 {kg/s}) < 0.0005 {kg/s};
366 END self_test;
367 END example_potter_wiggert_ex_9_7;
368
369 (*
370 Example 9.8 from Potter & Wiggert 3rd SI ed.
371 We have a converging-diverging nozzle with throat 5 cm dia and outlet 10 cm
372 dia. A normal shock is established where the nozzle is 7.5 cm dia. The
373 reservoir is at 200 kPa and 20 °C. We have to calculate the pressure of
374 the receiver that will cause this shock location, and also the pressure at
375 the nozzle exit plane.
376
377 Tested, works OK -- JP
378
379 FIXME this model shows a limitation with repeated use of the
380 isentropic_stagnation model. That model contains a stag.Vel = 0 equation
381 which causes structural singularity if multiple stagnation models are
382 coupled to the same stagnation state.
383
384 FIXME another interesting issue with this model was that the isentropic
385 flow equations between the normal shock and the receiver were difficult
386 to calculate using the isentropic_stagnation model in current form. While we
387 would naturally have written
388 exit.mdot, rec.mdot ARE_THE_SAME;
389 we found that that did not solve well, and instead the following worked much
390 better:
391 exit.A_star, rec.A_star ARE_THE_SAME;
392 Maybe we need a simpler model for an isentropic flow segment that does not
393 incorporate the stagnation state...?
394 *)
395 MODEL example_potter_wiggert_ex_9_8;
396 noz IS_A isentropic_stagnation;
397 A_star ALIASES noz.A_star;
398 A ALIASES noz.A;
399
400 D_star, D, D_exit IS_A distance;
401 0.25{PI}*D_star^2 = A_star;
402 0.25{PI}*D^2 = A;
403
404 p_0 ALIASES noz.p_0;
405 T_0 ALIASES noz.T_0;
406
407 shock IS_A air_duct_normal_shock;
408 shock.S1, noz.flow ARE_THE_SAME;
409
410 rec IS_A isentropic_stagnation;
411 rec.flow, shock.S2 ARE_THE_SAME;
412
413 exit IS_A isentropic_stagnation;
414 exit.stag, rec.stag ARE_THE_SAME;
415 A_exit ALIASES exit.A;
416 0.25{PI}*D_exit^2 = A_exit;
417 exit.A_star, rec.A_star ARE_THE_SAME;
418
419 p_exit ALIASES exit.flow.p;
420 p_rec ALIASES rec.stag.p;
421 METHODS
422 METHOD on_load;
423 FIX D_star, D, D_exit;
424 D_star := 5 {cm};
425 D := 7.5 {cm};
426 D_exit := 10 {cm};
427
428 (* reservoir conditions *)
429 FIX p_0, T_0;
430 p_0 := 200 {kPa};
431 T_0 := 20 {K} + 273.15 {K};
432
433 noz.flow.M.lower_bound := 1;
434 noz.flow.M.upper_bound := 100;
435
436 exit.novel.included := FALSE;
437
438 exit.flow.M.lower_bound := 0;
439 exit.flow.M.upper_bound := 1;
440 END on_load;
441 METHOD self_test;
442 ASSERT abs(p_exit - 109 {kPa}) < 0.5 {kPa};
443 ASSERT abs(p_rec - 114 {kPa}) < 0.5 {kPa};
444 ASSERT abs(exit.M - 0.265) < 0.001;
445 ASSERT abs(exit.A_rel_star - 2.284) < 0.0005;
446 ASSERT abs(shock.S2.M - 0.531) < 0.0005;
447 END self_test;
448 END example_potter_wiggert_ex_9_8;
449
450
451 (*
452 This model is solves Example 9.9 from Potter & Wiggert 3rd SI ed. The
453 problem relates to the performance of a Pitot tube placed in a
454 supersonic flow stream. Given the flow pressure as well as the pitot-tube
455 stagnation pressure, we must determine the free-stream flow velocity. We
456 are also given the temperature at the stagnation point.
457
458 Tested, works OK -- JP
459 *)
460 MODEL example_potter_wiggert_ex_9_9;
461 free ALIASES shock.S1;
462 shock IS_A air_normal_shock;
463 pitot IS_A isentropic_stagnation;
464 shock.S2, pitot.flow.state ARE_THE_SAME;
465 stag ALIASES pitot.stag;
466 Vel ALIASES free.Vel;
467 p ALIASES free.p;
468 p_0 ALIASES stag.p;
469 T_0 ALIASES stag.T;
470 T ALIASES free.T;
471 METHODS
472 METHOD on_load;
473 FIX free.p, stag.p, stag.T;
474 free.p := 75 {kPa};
475 stag.p := 300 {kPa};
476 stag.T := 150 {K} + 273.15 {K};
477
478 (* arbitrarily fix the pitot area to make problem square *)
479 FIX pitot.A;
480 pitot.A := 1 {m^2};
481 END on_load;
482 METHOD self_test;
483 ASSERT abs(free.M - 1.65) < 0.005;
484 ASSERT abs(free.T - 274 {K}) < 0.5 {K};
485 ASSERT abs(free.Vel - 547 {m/s}) < 0.5 {m/s};
486 END self_test;
487 END example_potter_wiggert_ex_9_9;

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