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 REFINES ideal_gas_node;
	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;

(*-----------------------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;

(*------------------------------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.5 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;


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

(*
	Model of isentropic flow in a converging or diverging nozzle over an
	interval in which no shock is present.
*)
MODEL isentropic_flow_segment;
	inlet, outlet IS_A air_duct_node;
	k ALIASES inlet.k;
	inlet.mdot, outlet.mdot ARE_THE_SAME;

	(* conservation of energy *)
	0 = (outlet.Vel^2 - inlet.Vel^2) / 2 + k/(k-1)*(outlet.p * outlet.v - inlet.p * inlet.rho);

	(* isentropic <=> adiabatic <=> no heat transfer *)
	outlet.T / inlet.T = (outlet.p / inlet.p)^((k-1)/k);

	(* conservation of mass *)
	inlet.rho * inlet.A * inlet.Vel = outlet.rho * outlet.A * outlet.Vel;

END isentropic_flow_segment;


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 *)
	stag.T = flow.T * (1 + (k-1)/2 * flow.M^2);

	(* isentropic *)
	(*(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;

END isentropic_stagnation;

(*
	This model calculates the pressure at the throat (critical area) for an
	assumed upstream reservoir pressure, and isentropic flow inbetween.

	Reproduces a part of Potter & Wiggert Example 9.2.
*)
MODEL isentropic_critical_area;
	stag IS_A air_node;
	crit IS_A air_duct_node;
	k ALIASES crit.k;

	crit.T / stag.T = 2 / (k+1);
	(*crit.p / stag.p = (2/(k+1))^(k/(k-1));*)
	crit.rho / stag.rho = (2/(k+1))^(1/(k-1));

	(* stagnation: no velocity *)
	stag.Vel = 0;

	(* supersonic at the critical area *)
	crit.M = 1;
METHODS
METHOD on_load;
	FIX stag.p, stag.T;
	stag.p := 500 {kPa};
	stag.T := 20 {K} + 273.15 {K};

	FIX crit.A;
	crit.A := 10 {cm^2};
END on_load;
END isentropic_critical_area;

MODEL isentropic_stagnation_test_1 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 mdot;
	mdot := 1.167 {kg/s};
END on_load;
END isentropic_stagnation_test_1;

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

	FIX flow.M;
	flow.M := 1;
END on_load;
END isentropic_stagnation_test_2;


(*
	Example 9.7 from Potter & Wiggert, 3rd SI edition. A converging-diverging
	nozzle with normal shock at the exit plane, exiting to ambient conditions.
	Exit diameter and throat diameter given. We calculate receiver pressure and
	flow rate.

	Note:
	* throat must be at M = 1
	* isentropic flow from throat to exit occurs with area ratio given
	* upstream (reservoir) conditions are 90 kPa, 20 °C - ambient.
*)
MODEL example_potter_wiggert_ex_9_7 REFINES air_normal_shock;

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

	FIXME this model isn't working yet.
*)
MODEL example_potter_wiggert_ex_9_2;
	(* wire up a reservoir, through a flowing node, to a receiver *)
	res, rec IS_A isentropic_stagnation;
	inlet ALIASES res.stag;
	outlet ALIASES rec.stag;
	res.flow, rec.flow ARE_THE_SAME;
	exit ALIASES res.flow;
METHODS
METHOD on_load;

	FIX inlet.p, inlet.T;
	inlet.p := 500 {kPa};
	inlet.T := 20 {K} + 273.15 {K};

	FIX exit.A;
	exit.A := 10 {cm^2};

	FIX outlet.p;
	outlet.p := 300 {kPa};
END on_load;
END example_potter_wiggert_ex_9_2;
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 REFINES ideal_gas_node;
65	A IS_A area;
66	mdot IS_A mass_rate;
67	mdot = rho * A * Vel;
68	END ideal_gas_duct_node;
69
70
71	(* Specific cases for air...*)
72
73	MODEL air_node REFINES ideal_gas_node;
74	MM :== 28.97 {kg/kmol};
75	k :== 1.4;
76	END air_node;
77
78	MODEL air_duct_node REFINES ideal_gas_duct_node;
79	MM :== 28.97 {kg/kmol};
80	k :== 1.4;
81	END air_duct_node;
82
83	(-----------------------NORMAL SHOCK IN ISENTROPIC FLOW----------------------)
84
85	(*
86	Model of a stationary normal shock in air. Using a moving frame of
87	reference, it can be used to model a moving shock wave as well.
88
89	Equations from Potter & Wiggert, 3rd SI edition, sects 9.1, 9.2 and 9.4.
90	*)
91	MODEL air_normal_shock;
92	S1, S2 IS_A air_node;
93	k ALIASES S1.k;
94
95	S2.M^2 = ( S1.M^2 + 2/(k-1) ) / (2k/(k-1)S1.M^2 - 1);
96
97	S2.p / S1.p = 2k/(k+1)S1.M^2 - (k-1)/(k+1);
98
99	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);
100	END air_normal_shock;
101
102	(------------------------------EXAMPLE MODELS -------------------------------)
103
104	(*
105	This model reproduces the results of Example 9.5 from Potter & Wiggert,
106	3rd SI edition. The question asks to determine the pressure and temperature
107	conditions downstream of a shock wave passing through ambient air of given
108	state.
109
110	Tested, works OK -- JP
111	*)
112	MODEL example_potter_wiggert_ex_9_5 REFINES air_normal_shock;
113	METHODS
114	METHOD on_load;
115	FIX S1.Vel, S1.p, S1.T;
116	S1.Vel := 450 {m/s};
117	S1.p := 80 {kPa};
118	S1.T := 15 {K} + 273.15 {K};
119	END on_load;
120	END example_potter_wiggert_ex_9_5;
121
122	(*
123	This model reproduces the results of Example 9.5 from Potter & Wiggert,
124	3rd SI edition. This problem shows the wind speeds implicit behind a strong
125	shock wave such as that arising from a high-powered bomb explosions.
126
127	Although the problem as given in P&W can be solved even without doing so,
128	we have added an assumption that the ambient air pressure is 101.3 kPa. This
129	allows the model to be 'square' and the pressure and temperature behind the
130	shock wave (4.83 bar, 500 K) to also be calculated.
131
132	Tested, works OK -- JP
133	*)
134	MODEL example_potter_wiggert_ex_9_6 REFINES air_normal_shock;
135	Vel_induced IS_A speed;
136	Vel_induced = S2.Vel - S1.Vel;
137	METHODS
138	METHOD on_load;
139	FIX S1.Vel, S1.T, S1.p;
140	S1.Vel := 700 {m/s};
141	S1.T := 15 {K} + 273.15 {K};
142	S1.p := 101.3 {kPa};
143	END on_load;
144	END example_potter_wiggert_ex_9_6;
145
146
147	(----------------------------STAGNATION CONDITIONS---------------------------)
148
149	(*
150	Model of isentropic flow in a converging or diverging nozzle over an
151	interval in which no shock is present.
152	*)
153	MODEL isentropic_flow_segment;
154	inlet, outlet IS_A air_duct_node;
155	k ALIASES inlet.k;
156	inlet.mdot, outlet.mdot ARE_THE_SAME;
157
158	(* conservation of energy *)
159	0 = (outlet.Vel^2 - inlet.Vel^2) / 2 + k/(k-1)(outlet.p outlet.v - inlet.p * inlet.rho);
160
161	(* isentropic <=> adiabatic <=> no heat transfer *)
162	outlet.T / inlet.T = (outlet.p / inlet.p)^((k-1)/k);
163
164	(* conservation of mass *)
165	inlet.rho * inlet.A * inlet.Vel = outlet.rho * outlet.A * outlet.Vel;
166
167	END isentropic_flow_segment;
168
169
170	MODEL isentropic_stagnation;
171	flow IS_A air_duct_node;
172	stag IS_A air_node;
173
174	(* conservation of mass *)
175	stag.Vel = 0;
176	mdot ALIASES flow.mdot;
177	k ALIASES flow.k;
178
179	(* conservation of energy *)
180	stag.T = flow.T * (1 + (k-1)/2 * flow.M^2);
181
182	(* isentropic *)
183	((1 + (k-1)/2 flow.M^2 ) * flow.p^((k-1)/k) = stag.p^((k-1)/k);*)
184	(stag.rho / flow.rho)^(k-1) = 1 + (k-1)/2*flow.M^2;
185
186	END isentropic_stagnation;
187
188	(*
189	This model calculates the pressure at the throat (critical area) for an
190	assumed upstream reservoir pressure, and isentropic flow inbetween.
191
192	Reproduces a part of Potter & Wiggert Example 9.2.
193	*)
194	MODEL isentropic_critical_area;
195	stag IS_A air_node;
196	crit IS_A air_duct_node;
197	k ALIASES crit.k;
198
199	crit.T / stag.T = 2 / (k+1);
200	(crit.p / stag.p = (2/(k+1))^(k/(k-1));)
201	crit.rho / stag.rho = (2/(k+1))^(1/(k-1));
202
203	(* stagnation: no velocity *)
204	stag.Vel = 0;
205
206	(* supersonic at the critical area *)
207	crit.M = 1;
208	METHODS
209	METHOD on_load;
210	FIX stag.p, stag.T;
211	stag.p := 500 {kPa};
212	stag.T := 20 {K} + 273.15 {K};
213
214	FIX crit.A;
215	crit.A := 10 {cm^2};
216	END on_load;
217	END isentropic_critical_area;
218
219	MODEL isentropic_stagnation_test_1 REFINES isentropic_stagnation;
220	METHODS
221	METHOD on_load;
222	FIX stag.p, stag.T, flow.A;
223	stag.p := 500 {kPa};
224	stag.T := 20 {K} + 273.15 {K};
225	flow.A := 10 {cm^2};
226
227	FIX mdot;
228	mdot := 1.167 {kg/s};
229	END on_load;
230	END isentropic_stagnation_test_1;
231
232	MODEL isentropic_stagnation_test_2 REFINES isentropic_stagnation;
233	METHODS
234	METHOD on_load;
235	FIX stag.p, flow.A;
236	stag.p := 500 {kPa};
237	stag.T := 20 {K} + 273.15 {K};
238	flow.A := 10 {cm^2};
239
240	FIX flow.M;
241	flow.M := 1;
242	END on_load;
243	END isentropic_stagnation_test_2;
244
245
246
247
248	(*
249	Example 9.7 from Potter & Wiggert, 3rd SI edition. A converging-diverging
250	nozzle with normal shock at the exit plane, exiting to ambient conditions.
251	Exit diameter and throat diameter given. We calculate receiver pressure and
252	flow rate.
253
254	Note:
255	* throat must be at M = 1
256	* isentropic flow from throat to exit occurs with area ratio given
257	* upstream (reservoir) conditions are 90 kPa, 20 °C - ambient.
258	*)
259	MODEL example_potter_wiggert_ex_9_7 REFINES air_normal_shock;
260
261
262
263
264	(*
265	Flow from reservoir at 20°C, 500 kPa to receiver at 300 kPa.
266	Potter & Wiggert 3rd SI Edition, Example 9.2.
267
268	FIXME this model isn't working yet.
269	*)
270	MODEL example_potter_wiggert_ex_9_2;
271	(* wire up a reservoir, through a flowing node, to a receiver *)
272	res, rec IS_A isentropic_stagnation;
273	inlet ALIASES res.stag;
274	outlet ALIASES rec.stag;
275	res.flow, rec.flow ARE_THE_SAME;
276	exit ALIASES res.flow;
277	METHODS
278	METHOD on_load;
279
280	FIX inlet.p, inlet.T;
281	inlet.p := 500 {kPa};
282	inlet.T := 20 {K} + 273.15 {K};
283
284	FIX exit.A;
285	exit.A := 10 {cm^2};
286
287	FIX outlet.p;
288	outlet.p := 300 {kPa};
289	END on_load;
290	END example_potter_wiggert_ex_9_2;