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

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