Chapter 5
[1]:
import sympy as sp
from sympy import symbols, sqrt, init_printing
from sympy_equation import Eqn, solve, table_of_expressions
from ambiance import Atmosphere
import numpy as np
import pint
import pygasflow
from pygasflow import *
from pygasflow.atd import *
init_printing()
pygasflow.defaults.solver_to_dict = True
ureg = pint.UnitRegistry()
# use "~P" to format units with unicode
ureg.formatter.default_format = "~"
pygasflow.defaults.pint_ureg = ureg
K = ureg.K
m = ureg.m
km = ureg.km
s = ureg.s
J = ureg.J
W = ureg.W
kg = ureg.kg
deg = ureg.deg
atm = ureg.atm
P 5.2
[2]:
ht = 20e06 * J / kg
T_exit = 1000 * K
Cp = 1004 * J/kg/K
R = 287.05 * J / kg / K
(a)
[3]:
hts, Cps, Ts, us = symbols("h_t, C_p, T, u", positive=True, real=True)
e = Eqn(hts, Cps * Ts + us**2 / 2)
e
[3]:
$\displaystyle h_{t} = C_{p} T + \frac{u^{2}}{2}$
Outside the reservoir, in order to achieve the maximum speed it must be \(T_{t} = 0\):
[4]:
u = (e.subs({hts: ht.magnitude, Ts: 0}) * 2).swap.apply(sp.sqrt)
u
[4]:
$\displaystyle u = 6324.55532033676$
[5]:
u_max = float(u.rhs) * m / s
u_max
[5]:
6324.555320336759 m/s
(b)
[6]:
f = 6 # Lighthill model
gamma = (f + 2) / f
Cp = (f + 2) / 2 * R
[7]:
gamma
[7]:
$\displaystyle 1.33333333333333$
[8]:
Cp
[8]:
1148.2 J/(K kg)
[9]:
u = solve(e.subs({hts: ht.magnitude, Ts: T_exit.magnitude, Cps: Cp.magnitude}))[0]
u
[9]:
$\displaystyle u = 6140.32572425926$
[10]:
u = float(u.rhs) * m / s
u
[10]:
6140.32572425926 m/s
(c)
If 20% of the entalphy is frozen, then: \(h_{t, avail} = 0.8 h_{t}\):
[11]:
ht_avail = 0.8 * ht
ht_avail
[11]:
16000000.0 J/kg
[12]:
u = solve(e.subs({hts: ht_avail.magnitude, Ts: T_exit.magnitude, Cps: Cp.magnitude}))[0]
u = float(u.rhs) * m / s
u
[12]:
5450.100916496868 m/s
P 5.3
[13]:
v_inf = 2 * km / s
H = 30 * km
(a)
[14]:
atmosphere = Atmosphere(H.to("m").magnitude)
T = atmosphere.temperature[0] * K
T
[14]:
226.50908361133006 K
[15]:
gamma = 1.4
R = 287.05 * J / kg / K
a = sound_speed(gamma, R, T).to("m/s")
a
[15]:
301.70715177284944 m/s
[16]:
M_inf = v_inf.to("m/s") / a
M_inf
[16]:
6.628944618143386
[17]:
res = isentropic_solver("m", M_inf, gamma=gamma)
res.show()
key quantity
----------------------------
m M 6.62894462
pr P / P0 0.00034079
dr rho / rho0 0.00333580
tr T / T0 0.10215985
prs P / P* 0.00064508
drs rho / rho* 0.00526203
trs T / T* 0.12259182
urs U / U* 2.32099998
ars A / A* 81.87873693
ma Mach Angle [deg] 8.67639576
pm Prandtl-Meyer [deg] 88.92735528
[18]:
Tt = (1 / res["tr"]) * T
Tt
[18]:
2217.2025914621277 K
(b)
[19]:
gamma_eff = 1.2
res = isentropic_solver("m", M_inf, gamma=gamma_eff)
Tt = (1 / res["tr"]) * T
Tt
[19]:
1221.8558375367288 K
(c)
[20]:
gamma_eff = 1 + 1e-06
res = isentropic_solver("m", M_inf, gamma=gamma_eff)
Tt = (1 / res["tr"]) * T
Tt
[20]:
226.51406034509927 K
[ ]: