models/johnpye/compressible_flow.a4c

(*	ASCEND modelling environment
	Copyright (C) 2011 Carnegie Mellon University

	This program is free software; you can redistribute it and/or modify
	it under the terms of the GNU General Public License as published by
	the Free Software Foundation; either version 2, or (at your option)
	any later version.

	This program is distributed in the hope that it will be useful,
	but WITHOUT ANY WARRANTY; without even the implied warranty of
	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
	GNU General Public License for more details.
	
	You should have received a copy of the GNU General Public License
	along with this program; if not, write to the Free Software
	Foundation, Inc., 59 Temple Place - Suite 330,
	Boston, MA 02111-1307, USA.
*)(*
	Calculations of isentropic shocks in air flow, following approach given
	in Potter & Wiggert, 3rd SI edition.

	Intention is to be able to generalise these models for arbitrary fluids as
	supported by FPROPS. But the present models include ideal gas and constant
	specific heat assumptions.

	Author: John Pye
*)

REQUIRE "johnpye/thermo_types.a4c";


(*
	Ideal gas flow node, with ideal gas law, Mach number and speed of sound
	included. No allowance for enthalpy, entropy at this stage.
*)
MODEL ideal_gas_node;
	M "Mach number" IS_A positive_factor;
	p IS_A pressure;
	T IS_A temperature;
	v IS_A specific_volume;
	rho IS_A mass_density;
	Vel IS_A speed;

	MM IS_A molar_weight_constant;
	R IS_A specific_gas_constant;
	R :== 1{GAS_C}/MM;

	rho*v = 1;
	p*v = R * T;

	k "adiabatic index" IS_A real_constant;

	c IS_A speed;
	c = sqrt(k*R*T);

	M = Vel / c;
END ideal_gas_node;


(*
	Add cross-sectional area to the above model; not needed for normal shock
	calculations, but useful in nozzle calcs.
*)
MODEL ideal_gas_duct_node;
	(* we incorporate state as submodel so we can ARE_THE_SAME non-duct nodes *)
	state IS_A ideal_gas_node;
	M ALIASES state.M;
	p ALIASES state.p;
	T ALIASES state.T;
	v ALIASES state.v;
	rho ALIASES state.rho;
	Vel ALIASES state.Vel;
	MM ALIASES state.MM;
	R ALIASES state.R;
	k ALIASES state.k;
	c ALIASES state.c;

	A IS_A area;
	mdot IS_A mass_rate;
	mdot = rho * A * Vel;
END ideal_gas_duct_node;


(* Specific cases for air...*)

MODEL air_node REFINES ideal_gas_node;
	MM :== 28.97 {kg/kmol};
	k :== 1.4;
END air_node;

MODEL air_duct_node REFINES ideal_gas_duct_node;
	MM :== 28.97 {kg/kmol};
	k :== 1.4;
END air_duct_node;

(*----------------------------STAGNATION CONDITIONS---------------------------*)


MODEL isentropic_stagnation;
	flow IS_A air_duct_node;
	stag IS_A air_node;

	(* conservation of mass *)
	stag.Vel = 0;
	mdot ALIASES flow.mdot;
	k ALIASES flow.k;

	(* conservation of energy *)
	1 = T_rel_0 * (1 + (k-1)/2 * flow.M^2);

	(* isentropic *)
	isen:  (1 + (k-1)/2*flow.M^2)^(k/(k-1)) * p_rel_0 = 1;
	(*(1 + (k-1)/2 * flow.M^2 ) * flow.p^((k-1)/k) = stag.p^((k-1)/k);*)
	(*(stag.rho / flow.rho)^(k-1) = 1 + (k-1)/2*flow.M^2;*)

	p_rel_0 IS_A positive_factor;
	T_rel_0 IS_A positive_factor;
	A_rel_star IS_A positive_factor;
	A_star IS_A area;
	p_0 ALIASES stag.p;
	p ALIASES flow.p;
	T ALIASES flow.T;
	T_0 ALIASES stag.T;
	A ALIASES flow.A;
	M ALIASES flow.M;

	p_rel_0 * stag.p = flow.p;
	T_rel_0 * stag.T = flow.T;
	A_rel_star * A_star = flow.A;

	arat:  A_rel_star * flow.M = ((2 + (k-1)*M^2)/(k+1))^((k+1)/(2*(k-1)));

END isentropic_stagnation;

MODEL isentropic_stagnation_test REFINES isentropic_stagnation;
METHODS
METHOD on_load;
	FIX stag.p, stag.T, flow.A;
	stag.p := 500 {kPa};
	stag.T := 20 {K} + 273.15 {K};
	flow.A := 10 {cm^2};

	FIX M;
	M := 0.80;
END on_load;
END isentropic_stagnation_test;

(*------------------------------EXAMPLE MODELS -------------------------------*)

(*
	Flow from reservoir at 20°C, 500 kPa to receiver at 200 and 300 kPa.
	Potter & Wiggert 3rd SI Edition, Example 9.2.

	We assume the flow is through a converging-only nozzle from the reservoir 
	into the receiver. As such, either the entire flow is subsonic, or else
	there will be choked flow with M = 1 in the exit of the nozzle.

	'crit' uses upstream conditions to determine the reservoir pressure that
	will give M = 1 at the exit. Then, in part (a), the receiver pressure is
	higher than that critical pressure, so subsonic flow occurs. For (b), the
	pressure is lower, so choked flow occurs, with flow rate limited by the
	M_exit = 1 choked-flow limit.

	Tested, works OK -- JP
*)
MODEL example_potter_wiggert_ex_9_2;
	crit, part_a, part_b IS_A isentropic_stagnation;
	crit.A, part_a.A, part_b.A ARE_THE_SAME;
	crit.T_0, part_a.T_0, part_b.T_0 ARE_THE_SAME;
	A ALIASES crit.A;
	T_0 ALIASES crit.T_0;
METHODS
METHOD on_load;
	FIX A; A := 10 {cm^2};
	FIX T_0; T_0 := 20 {K} + 273.15 {K};

	FIX crit.M; crit.M := 1;
	FIX crit.p_0; crit.p_0 := 500 {kPa};

	FIX part_a.p_0; part_a.p_0 := 500 {kPa};
	FIX part_a.p; part_a.p := 300 {kPa};

	(* actually, part (b) ends up being identical to 'crit'! *)
	FIX part_b.p_0; part_b.p_0 := 500 {kPa};
	FIX part_b.M; part_b.M := 1;
END on_load;
END example_potter_wiggert_ex_9_2;


(*
	Potter & Wiggert 3rd SI ed, Example 9.3.
	Converging-diverging nozzle with exit 40 cm² and throat 10 cm², with flow
	coming from reservoir at 20 °C, 500 kPa. Two different exit pressures can
	result in M = 1 at the throat -- we calculate what those pressure are.

	For this problem, saturation states are the upstream reservoir, and we
	assert that A_star is the throat area and A is the exit area. These areas
	then yield the Mach number from eq 9.3.19 (iteratively). We find the two
	solutions by having two instances of the 'isentropic_saturation' model, and
	impose upper and lower bounds on M in each of the two models.

	Tested, works OK -- JP
*)
MODEL example_potter_wiggert_ex_9_3;
	sub, sup IS_A isentropic_stagnation;
	sub.A_star, sup.A_star ARE_THE_SAME;
	A_star ALIASES sub.A_star;
	sub.A, sup.A ARE_THE_SAME;
	A ALIASES sub.A;
	sub.p_0, sup.p_0 ARE_THE_SAME;
	p_0 ALIASES sub.p_0;
	sub.T_0, sup.T_0 ARE_THE_SAME;
	T_0 ALIASES sub.T_0;
METHODS
METHOD on_load;
	FIX A, A_star;
	A_star := 10 {cm^2};
	A := 40 {cm^2};

	(* ensure two different solution regions *)
	sub.M.upper_bound := 1;
	sup.M.lower_bound := 1;

	FIX T_0, p_0;
	T_0 := 20 {K} +  273.15 {K};
	p_0 := 500 {kPa};
END on_load;
END example_potter_wiggert_ex_9_3;

(*-----------------------NORMAL SHOCK IN ISENTROPIC FLOW----------------------*)

(*
	Model of a stationary normal shock in air. Using a moving frame of
	reference, it can be used to model a moving shock wave as well.

	Equations from Potter & Wiggert, 3rd SI edition, sects 9.1, 9.2 and 9.4.
*)
MODEL air_normal_shock;
	S1, S2 IS_A air_node;
	k ALIASES S1.k;
	
	S2.M^2 = ( S1.M^2 + 2/(k-1) ) / (2*k/(k-1)*S1.M^2 - 1);

	S2.p / S1.p = 2*k/(k+1)*S1.M^2 - (k-1)/(k+1);

	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);
END air_normal_shock;

MODEL air_duct_normal_shock;
	S1, S2 IS_A air_duct_node;
	shock IS_A air_normal_shock;
	S1.state, shock.S1 ARE_THE_SAME;
	S2.state, shock.S2 ARE_THE_SAME;
	S1.A, S2.A ARE_THE_SAME;
END air_duct_normal_shock;

(*------------------------------EXAMPLE MODELS -------------------------------*)

(*
	This model reproduces the results of Example 9.5 from Potter & Wiggert, 
	3rd SI edition. The question asks to determine the pressure and temperature 
	conditions downstream of a shock wave passing through ambient air of given
	state.

	Tested, works OK -- JP
*)
MODEL example_potter_wiggert_ex_9_5 REFINES air_normal_shock;
METHODS
METHOD on_load;
	FIX S1.Vel, S1.p, S1.T;
	S1.Vel := 450 {m/s};
	S1.p := 80 {kPa};
	S1.T := 15 {K} + 273.15 {K};
END on_load;
END example_potter_wiggert_ex_9_5;

(*
	This model reproduces the results of Example 9.6 from Potter & Wiggert, 
	3rd SI edition. This problem shows the wind speeds implicit behind a strong
	shock wave such as that arising from a high-powered bomb explosions.

	Although the problem as given in P&W can be solved even without doing so,
	we have added an assumption that the ambient air pressure is 101.3 kPa. This
	allows the model to be 'square' and the pressure and temperature behind the
	shock wave (4.83 bar, 500 K) to also be calculated.

	Tested, works OK -- JP
*)
MODEL example_potter_wiggert_ex_9_6 REFINES air_normal_shock;
	Vel_induced IS_A speed;
	Vel_induced = S2.Vel - S1.Vel;
METHODS
METHOD on_load;
	FIX S1.Vel, S1.T, S1.p;
	S1.Vel := 700 {m/s};
	S1.T := 15 {K} + 273.15 {K};
	S1.p := 101.3 {kPa};
END on_load;
END example_potter_wiggert_ex_9_6;


(*
	This model reproduces the results of Example 9.7 from Potter & Wiggert,
	3rd SI edition. The problem relates to the calculation of outlet conditions 
	and flow rate for a converging-diverging nozzle, given specified upstream
	reservoir pressure and temperature, given the throat and exit areas, and
	subject to the fact that there is a normal shock at the nozzle's exit plane.

	For a shock to exist at the exit, we know that we must have M = 1 at the
	throat, so A* = A_t.

	Tested, works OK -- JP
*)
MODEL example_potter_wiggert_ex_9_7; 
	noz IS_A isentropic_stagnation;
	A_star ALIASES noz.A_star;
	A ALIASES noz.A;

	D_star, D IS_A distance;
	0.25{PI}*D_star^2 = A_star;
	0.25{PI}*D^2 = A;

	p_0 ALIASES noz.p_0;
	T_0 ALIASES noz.T_0;

	shock IS_A air_duct_normal_shock;
	shock.S1, noz.flow ARE_THE_SAME;

	rec IS_A isentropic_stagnation;
	rec.flow, shock.S2 ARE_THE_SAME;

	p_2 ALIASES rec.p;
	mdot ALIASES noz.mdot;

METHODS
METHOD on_load;
	(* nozzle geometry *)
	FIX D_star, D;
	D_star := 5 {cm};
	D := 10 {cm};

	(* note that isentropic flow equations have to be bounded over a reasonable
	range for the Brent	algorithm to converge in this case. Note that P&W Table
	D.1 gives M up to 10 only, so limiting to 100 is still more generous.
	FIXME Perhaps we can make these models more user-friendly by having a
	supersonic and subsonic model that sets these bounds automatically...? *)
	noz.flow.M.lower_bound := 1;
	noz.flow.M.upper_bound := 100;

	(* reservoir conditions *)
	FIX p_0, T_0;
	p_0 := 90 {kPa};
	T_0 := 20 {K} + 273.15 {K};
END on_load;
METHOD self_test;
	(* values from P&W *)
	ASSERT abs(p_2 - 26.6 {kPa}) < 0.05 {kPa};
	ASSERT abs(mdot - 0.417 {kg/s}) < 0.0005 {kg/s};
END self_test;
END example_potter_wiggert_ex_9_7;


MODEL example_potter_wiggert_ex_9_8;
	noz IS_A isentropic_stagnation;
	A_star ALIASES noz.A_star;
	A ALIASES noz.A;

	D_star, D IS_A distance;
	0.25{PI}*D_star^2 = A_star;
	0.25{PI}*D^2 = A;

	p_0 ALIASES noz.p_0;
	T_0 ALIASES noz.T_0;

	shock IS_A air_duct_normal_shock;
	shock.S1, noz.flow ARE_THE_SAME;

	rec IS_A isentropic_stagnation;
	rec.flow, shock.S2 ARE_THE_SAME;
METHODS
METHOD on_load;
	FIX D_star, D;
	D_star := 5 {cm};
	D := 7.5 {cm};

	(* reservoir conditions *)
	FIX p_0, T_0;
	p_0 := 200 {kPa};
	T_0 := 20 {K} + 273.15 {K};

	noz.flow.M.lower_bound := 1;
	noz.flow.M.upper_bound := 100;
END on_load;
END example_potter_wiggert_ex_9_8;
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	pv = R T;
50
51	k "adiabatic index" IS_A real_constant;
52
53	c IS_A speed;
54	c = sqrt(kRT);
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	isen: (1 + (k-1)/2flow.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)/2flow.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	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_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) ) / (2k/(k-1)S1.M^2 - 1);
243
244	S2.p / S1.p = 2k/(k+1)S1.M^2 - (k-1)/(k+1);
245
246	S2.T / S1.T = ( 1 + (k-1)/2S1.M^2) (2k/(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	MODEL example_potter_wiggert_ex_9_8;
364	noz IS_A isentropic_stagnation;
365	A_star ALIASES noz.A_star;
366	A ALIASES noz.A;
367
368	D_star, D IS_A distance;
369	0.25{PI}*D_star^2 = A_star;
370	0.25{PI}*D^2 = A;
371
372	p_0 ALIASES noz.p_0;
373	T_0 ALIASES noz.T_0;
374
375	shock IS_A air_duct_normal_shock;
376	shock.S1, noz.flow ARE_THE_SAME;
377
378	rec IS_A isentropic_stagnation;
379	rec.flow, shock.S2 ARE_THE_SAME;
380	METHODS
381	METHOD on_load;
382	FIX D_star, D;
383	D_star := 5 {cm};
384	D := 7.5 {cm};
385
386	(* reservoir conditions *)
387	FIX p_0, T_0;
388	p_0 := 200 {kPa};
389	T_0 := 20 {K} + 273.15 {K};
390
391	noz.flow.M.lower_bound := 1;
392	noz.flow.M.upper_bound := 100;
393	END on_load;
394	END example_potter_wiggert_ex_9_8;