/[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 2461 - (show annotations) (download) (as text)
Mon Apr 25 09:09:57 2011 UTC (11 years, 2 months ago) by jpye
File MIME type: text/x-ascend
File size: 9486 byte(s)
Working on P&W Ex 9.7. Some minor cleanup of other files.
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 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 stag.Vel = 0;
105 mdot ALIASES flow.mdot;
106 k ALIASES flow.k;
107
108 (* conservation of energy *)
109 1 = T_rel_0 * (1 + (k-1)/2 * flow.M^2);
110
111 (* isentropic *)
112 (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 p_rel_0 * stag.p = flow.p;
128 T_rel_0 * stag.T = flow.T;
129 A_rel_star * A_star = flow.A;
130
131 A_rel_star = 1 / 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_saturation' 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.5 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 This model reproduces the results of Example 9.7 from Potter & Wiggert,
303 3rd SI edition. The problem relates to the calculation of outlet conditions
304 and flow rate for a converging-diverging nozzle, given specified upstream
305 reservoir pressure and temperature, given the throat and exit areas, and
306 subject to the fact that there is a normal shock at the nozzle's exit plane.
307
308 For a shock to exist at the exit, we know that we must have M = 1 at the
309 throat, so A* = A_t.
310 *)
311 MODEL example_potter_wiggert_ex_9_7;
312 noz IS_A isentropic_stagnation;
313 A_star ALIASES noz.A_star;
314 A ALIASES noz.A;
315 p_0 ALIASES noz.p_0;
316 T_0 ALIASES noz.T_0;
317
318 shock IS_A air_duct_normal_shock;
319 shock.S1, noz.flow ARE_THE_SAME;
320
321 rec IS_A isentropic_stagnation;
322 rec.flow, shock.S2 ARE_THE_SAME;
323
324 p_2 ALIASES rec.p;
325 mdot ALIASES noz.mdot;
326 METHODS
327 METHOD on_load;
328 FIX A_star, A;
329 A_star := 5 {cm^2};
330 A := 10 {cm^2};
331 FIX p_0, T_0;
332 p_0 := 90 {kPa};
333 T_0 := 20 {K} + 273.15 {K};
334 END on_load;
335 END example_potter_wiggert_ex_9_7;

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