Chapter 3 Problems
of “Modern Compressible Flow With Historical Perspective”, Third Edition, John D. Anderson Jr
[1]:
import pint
import pygasflow
from pygasflow.solvers import *
import numpy as np
ureg = pint.UnitRegistry()
# to imperial units
ureg.define("pound_mass = 0.45359237 kg = lbm")
# use "~P" to format units with unicode
ureg.formatter.default_format = "~"
# let pygasflow knows which UnitRegistry to use
pygasflow.defaults.pint_ureg = ureg
# dictionary results offers improved code readability
pygasflow.defaults.solver_to_dict = True
# shorcuts for conveniance
kg = ureg.kg
K = ureg.K
degR = ureg.degR
degC = ureg.degC
J = ureg.J
m = ureg.m
atm = ureg.atm
Pa = ureg.pascal
sec = ureg.s
lb = ureg.lb
ft = ureg.feet
ft_lb = ureg.ft_lb
slug = ureg.slug
Q_ = ureg.Quantity
gamma = 1.4
R = 287.05 * J / (kg * K)
P 3.1
[2]:
M = 0.7
p = 0.9 * atm
T = 250 * K
[3]:
res = isentropic_solver("m", M, gamma, to_dict=True)
res.show()
key quantity
----------------------------
m M 0.70000000
pr P / P0 0.72092786
dr rho / rho0 0.79157879
tr T / T0 0.91074681
prs P / P* 1.36466537
drs rho / rho* 1.24866881
trs T / T* 1.09289617
urs U / U* 0.73179172
ars A / A* 1.09437268
ma Mach Angle nan
pm Prandtl-Meyer nan
[4]:
P0 = (1 / res["pr"]) * p
T0 = (1 / res["tr"]) * T
ps = (1 / res["prs"]) * p
Ts = (1 / res["trs"]) * T
a_star = np.sqrt(gamma * R * Ts)
[5]:
quantities = [P0, T0, ps, Ts, a_star]
labels = ["P0", "T0", "p*", "T*", "a*"]
for q, l in zip(quantities, labels):
print(l, "=", q)
P0 = 1.248391203471733 atm
T0 = 274.49999999999994 K
p* = 0.6595023367404417 atm
T* = 228.74999999999997 K
a* = 303.1959143854019 J ** 0.5 / kg ** 0.5
P 3.2
[6]:
p = 5e04 * Pa
T = 200 * K
p0 = 1.5e06 * Pa
[7]:
pr = p / p0
res = isentropic_solver("pressure", pr, gamma, to_dict=True)
res.show()
key quantity
----------------------------
m M 2.86585027
pr P / P0 0.03333333
dr rho / rho0 0.08808732
tr T / T0 0.37841240
prs P / P* 0.06309764
drs rho / rho* 0.13895254
trs T / T* 0.45409488
urs U / U* 1.93119797
ars A / A* 3.72654782
ma Mach Angle [deg] 20.42228556
pm Prandtl-Meyer [deg] 47.10098339
[8]:
print("Mach =", res["m"])
Mach = 2.8658502699286807
[9]:
T0 = (1 / res["tr"]) * T
T0
[9]:
528.5239107860116 K
P 3.3
[10]:
U = 3000 * ft / sec
T = 500 * degR
[11]:
a = np.sqrt(gamma * R * T.to("kelvin")).to("m/s")
a
[11]:
334.1115914714058 m/s
[12]:
M = U.to("m/s") / a
M
[12]:
2.736810165648673
[13]:
from pygasflow.generic import characteristic_mach_number
M_star = characteristic_mach_number(M)
M_star
[13]:
1.8968667428039343
P 3.4
[14]:
M1 = 3
p1 = 1 * atm
rho1 = 1.23 * kg / m**3
[15]:
T1 = (273.15 + 20) * K
T1
[15]:
293.15 K
[16]:
res_ise = isentropic_solver("m", M1, gamma)
res_ise.show()
key quantity
----------------------------
m M 3.00000000
pr P / P0 0.02722368
dr rho / rho0 0.07622631
tr T / T0 0.35714286
prs P / P* 0.05153250
drs rho / rho* 0.12024251
trs T / T* 0.42857143
urs U / U* 1.96396101
ars A / A* 4.23456790
ma Mach Angle 19.47122063
pm Prandtl-Meyer 49.75734674
[17]:
T01 = (1 / res_ise["tr"]) * T1
T01
[17]:
820.8199999999999 K
[18]:
p01 = (1 / res_ise["pr"]) * p1
p01
[18]:
36.732721804952064 atm
[19]:
res_shock = normal_shockwave_solver("mu", M1, gamma=gamma)
res_shock.show()
key quantity
---------------------
mu Mu 3.00000000
md Md 0.47519096
pr pd/pu 10.33333333
dr rhod/rhou 3.85714286
tr Td/Tu 2.67901235
tpr p0d/p0u 0.32834389
[20]:
p2 = res_shock["pr"] * p1
p2
[20]:
10.333333333333334 atm
[21]:
T2 = res_shock["tr"] * T1
T2
[21]:
785.3524691358024 K
[22]:
rho2 = res_shock["dr"] * rho1
rho2
[22]:
4.744285714285715 kg/m3
[23]:
p02 = res_shock["tpr"] * p01
p02
[23]:
12.060964701266622 atm
[24]:
T02 = T01
T02
[24]:
820.8199999999999 K
[25]:
u2 = (np.sqrt(gamma * R * T1) * res_shock["md"]).to("m / s")
u2
[25]:
163.1007341115972 m/s
P 3.5
[26]:
from pygasflow.shockwave import rayleigh_pitot_formula, m1_from_rayleigh_pitot_pressure_ratio
pr_sonic = rayleigh_pitot_formula(1, gamma)
pr_sonic
[26]:
np.float64(1.892929158737854)
(a)
[27]:
p02 = 1.22e05 * Pa
p1 = 1.01e05 * Pa
[28]:
pr = p02 / p1
pr
[28]:
1.2079207920792079
Because pr < pr_sonic, the Mach number must be subsonic. Hence, isentropic relations:
[29]:
res = isentropic_solver("pressure", p1 / p02)
res.show()
key quantity
----------------------------
m M 0.52656721
pr P / P0 0.82786885
dr rho / rho0 0.87377799
tr T / T0 0.94745903
prs P / P* 1.56709709
drs rho / rho* 1.37833320
trs T / T* 1.13695084
urs U / U* 0.56146754
ars A / A* 1.29217439
ma Mach Angle [deg] nan
pm Prandtl-Meyer [deg] nan
[30]:
M1 = res["m"]
M1
[30]:
np.float64(0.5265672087837046)
(b)
[31]:
p02 = 7222 * lb / ft**2
p1 = 2116 * lb / ft**2
[32]:
pr = p02 / p1
pr
[32]:
3.4130434782608696
Because pr > pr_sonic, there is a normal shock wave:
[33]:
M1 = m1_from_rayleigh_pitot_pressure_ratio(pr, gamma)
M1
[33]:
np.float64(1.4999388067617616)
(c)
[34]:
p02 = 13107 * lb / ft**2
p1 = 1020 * lb / ft**2
[35]:
pr = p02 / p1
pr
[35]:
12.85
Because pr > pr_sonic, there is a normal shock wave:
[36]:
M1 = m1_from_rayleigh_pitot_pressure_ratio(pr, gamma)
M1
[36]:
np.float64(3.1005608737676624)
P 3.6
[37]:
import matplotlib.pyplot as plt
import numpy as np
from pygasflow import shock_compression, isentropic_compression
import pint
ureg = pint.UnitRegistry()
p1 = 1 * atm
v1 = 1 * m**3
v2 = 0.3 * m**3
rho1 = 1 / v1
rho2 = 1 / v2
dr = np.linspace(rho1, rho2) / rho1
pr_ise, dr_ise, tr_ise = isentropic_compression(dr=dr, gamma=1.4, to_dict=False)
pr_shock, dr_shock = shock_compression(dr=dr, gamma=1.4, to_dict=False)
fig, ax = plt.subplots()
ax.plot(1 / (dr_ise * rho1), pr_ise * p1, label="isentropic")
ax.plot(1 / (dr_shock * rho1), pr_shock * p1, ":", label="normal shock")
ax.legend()
ax.set_xlabel("v [$m^{3}/kg$]")
ax.set_ylabel("Pressure [atm]")
plt.show()
/home/davide/Documents/Development/envs/pygas/lib/python3.12/site-packages/matplotlib/cbook.py:1355: UnitStrippedWarning: The unit of the quantity is stripped when downcasting to ndarray.
return np.asarray(x, float)
For a specified change in specific volume, the shock compression is stronger than the isentropic one. However, shock compression is less efficient due to the total pressure loss across the shock wave.
P 3.7
[38]:
M = 38
T = 270 * K
gamma = 1.4
[39]:
res = isentropic_solver("m", M, gamma)
res.show()
key quantity
----------------------------
m M 38.00000000
pr P / P0 0.00000000
dr rho / rho0 0.00000070
tr T / T0 0.00345066
prs P / P* 0.00000000
drs rho / rho* 0.00000110
trs T / T* 0.00414079
urs U / U* 2.44525992
ars A / A* 370653.24671053
ma Mach Angle 1.50795775
pm Prandtl-Meyer 122.92471032
[40]:
T0 = 1 / res["tr"] * T
T0
[40]:
78245.99999999999 K
P 3.8
[41]:
p1 = 1 * atm
T1 = 288 * K
gas = gas_solver("gamma", 1.4, "R", R)
Cp = gas["Cp"]
Cp
[41]:
1004.6750000000002 J/(K kg)
(a)
[42]:
M1 = 2
[43]:
res_ise = isentropic_solver("m", M1)
res_ise.show()
key quantity
----------------------------
m M 2.00000000
pr P / P0 0.12780453
dr rho / rho0 0.23004815
tr T / T0 0.55555556
prs P / P* 0.24192491
drs rho / rho* 0.36288737
trs T / T* 0.66666667
urs U / U* 1.63299316
ars A / A* 1.68750000
ma Mach Angle 30.00000000
pm Prandtl-Meyer 26.37976081
[44]:
T0 = (1 / res_ise["tr"]) * T1
T0
[44]:
518.4 K
[45]:
res = rayleigh_solver("m", M1)
res.show()
key quantity
---------------------
m M 2.00000000
prs P / P* 0.36363636
drs rho / rho* 0.68750000
trs T / T* 0.52892562
tprs P0 / P0* 1.50309598
ttrs T0 / T0* 0.79338843
urs U / U* 1.45454545
eps (s*-s) / R 1.21757521
[46]:
T0s = (1 / res["ttrs"]) * T0
T0s
[46]:
653.4 K
[47]:
Cp
[47]:
1004.6750000000002 J/(K kg)
[48]:
T0s
[48]:
653.4 K
[49]:
qres = specific_heat_solver(Cp=Cp, T01=T0, T02=T0s)
qres.show()
key quantity
--------------------------
q q [J / kg] 135631.12500000
Cp Cp [J / K / kg] 1004.67500000
T01 T01 [K] 518.40000000
T02 T02 [K] 653.40000000
DeltaT0 ΔT0 [K] 135.00000000
q_Cp q / Cp [K] 135.00000000
[50]:
qres["q"]
[50]:
135631.12500000003 J/kg
(b)
[51]:
M1 = 0.2
[52]:
res_ise = isentropic_solver("m", M1)
T0 = (1 / res_ise["tr"]) * T1
T0
[52]:
290.304 K
[53]:
res = rayleigh_solver("m", M1)
T0s = (1 / res["ttrs"]) * T0
T0s
[53]:
1672.7039999999995 K
[54]:
q = specific_heat_solver(Cp=Cp, T01=T0, T02=T0s)
q
[54]:
{'q': <Quantity(1388862.72, 'joule / kilogram')>,
'Cp': <Quantity(1004.675, 'joule / kilogram / kelvin')>,
'T01': <Quantity(290.304, 'kelvin')>,
'T02': <Quantity(1672.704, 'kelvin')>,
'DeltaT0': <Quantity(1382.4, 'kelvin')>,
'q_Cp': <Quantity(1382.4, 'kelvin')>}
P 3.9
[55]:
p1 = 10 * atm
T1 = 1000 * degR
M1 = 0.2
q_fuel = 4.5e08 * ft_lb / slug
f_a = 0.06
[56]:
res_ise = isentropic_solver("m", M1, gamma)
res_ise.show()
key quantity
----------------------------
m M 0.20000000
pr P / P0 0.97249670
dr rho / rho0 0.98027668
tr T / T0 0.99206349
prs P / P* 1.84086737
drs rho / rho* 1.54632859
trs T / T* 1.19047619
urs U / U* 0.21821789
ars A / A* 2.96352000
ma Mach Angle nan
pm Prandtl-Meyer nan
[57]:
T01 = (1 / res_ise["tr"]) * T1
T01
[57]:
1008.0 °R
[58]:
res_ray = rayleigh_solver("m", M1, gamma)
res_ray.show()
key quantity
---------------------
m M 0.20000000
prs P / P* 2.27272727
drs rho / rho* 11.00000000
trs T / T* 0.20661157
tprs P0 / P0* 1.23459588
ttrs T0 / T0* 0.17355372
urs U / U* 0.09090909
eps (s*-s) / R 6.34018207
[59]:
T1s = (1 / res_ray["trs"]) * T1
T1s
[59]:
4839.999999999999 °R
[60]:
p1s = (1 / res_ray["prs"]) * p1
p1s
[60]:
4.4 atm
[61]:
T01s = (1 / res_ray["ttrs"]) * T01
T01s
[61]:
5807.999999999999 °R
[62]:
gas = gas_solver("gamma", gamma, "R", R)
Cp = gas["Cp"]
Cp
[62]:
1004.6750000000002 J/(K kg)
[63]:
from sympy import symbols, solve, Eq
a, f = symbols("a, f", positive=True)
e1 = Eq(f / a, f_a)
e2 = Eq(1, f + a)
sol = solve([e1, e2], [f, a])
sol
[63]:
{a: 0.943396226415094, f: 0.0566037735849057}
[64]:
mdot = 100 * (slug / sec)
mdot_air = mdot * float(sol[a])
mdot_fuel = mdot * float(sol[f])
print("Air mass flow rate:", mdot_air)
print("Fuel mass flow rate:", mdot_fuel)
Air mass flow rate: 94.33962264150944 slug / s
Fuel mass flow rate: 5.660377358490567 slug / s
[65]:
qdot = q_fuel * mdot_fuel
qdot
[65]:
2547169811.320755 ft_lb/s
[66]:
q = qdot / mdot
q
[66]:
25471698.11320755 ft_lb/slug
[67]:
T02 = (q / Cp + T01).to("rankine")
T02
[67]:
5247.6961600742925 °R
[68]:
res_ray_2 = rayleigh_solver("total_temperature_super", T02 / T01s, gamma)
res_ray_2.show()
key quantity
---------------------
m M 1.52281760
prs P / P* 0.56516296
drs rho / rho* 0.76301054
trs T / T* 0.74070138
tprs P0 / P0* 1.13297391
ttrs T0 / T0* 0.90352895
urs U / U* 1.31059789
eps (s*-s) / R 0.47991089
[69]:
P2 = res_ray_2["prs"] * p1s
P2
[69]:
2.48671702192768 atm
[70]:
T2 = res_ray_2["trs"] * T1s
T2
[70]:
3584.9946793591835 °R
P 3.10
[71]:
rayleigh_solver("m", 1).show()
key quantity
---------------------
m M 1.00000000
prs P / P* 1.00000000
drs rho / rho* 1.00000000
trs T / T* 1.00000000
tprs P0 / P0* 1.00000000
ttrs T0 / T0* 1.00000000
urs U / U* 1.00000000
eps (s*-s) / R -0.00000000
[72]:
T02s = T01s
T02s
[72]:
5807.999999999999 °R
[73]:
T02 = T02s
T02
[73]:
5807.999999999999 °R
[74]:
qres = specific_heat_solver(Cp=Cp, T01=T01, T02=T02)
q = qres["q"].to("footpound / slug")
q
[74]:
28837951.19442091 ft_lb/slug
[75]:
qdot = mdot * q
qdot
[75]:
2883795119.442091 ft_lb/s
[76]:
mdot_fuel = qdot / q_fuel
mdot_fuel
[76]:
6.408433598760202 slug/s
[77]:
mdot_air = mdot - mdot_fuel
mdot_air
[77]:
93.5915664012398 slug/s
[78]:
new_f_a = mdot_fuel / mdot_air
new_f_a
[78]:
0.06847234045946377
P 3.11
[79]:
f_a = 0.03
T2 = 4800 * degR
gamma = 1.4
q_fuel = 4.5e08 * ft_lb / slug
[80]:
from sympy import symbols, solve, Eq
a, f = symbols("a, f", positive=True)
e1 = Eq(f / a, f_a)
e2 = Eq(1, f + a)
sol = solve([e1, e2], [f, a])
sol
[80]:
{a: 0.970873786407767, f: 0.0291262135922330}
[81]:
mdot = 100 * (slug / sec)
mdot_air = mdot * float(sol[a])
mdot_fuel = mdot * float(sol[f])
print("Air mass flow rate:", mdot_air)
print("Fuel mass flow rate:", mdot_fuel)
Air mass flow rate: 97.0873786407767 slug / s
Fuel mass flow rate: 2.912621359223301 slug / s
[82]:
qdot = q_fuel * mdot_fuel
qdot
[82]:
1310679611.6504855 ft_lb/s
[83]:
gas = gas_solver("gamma", gamma, "R", R)
Cp = gas["Cp"]
Cp = Cp.to("foot * force_pound / (lbm * degR)")
Cp
[83]:
186.7314425258093 ft lbf/(°R lbm)
[84]:
Delta_T0 = (qdot / mdot / Cp).to("degR")
Delta_T0
[84]:
2181.591227999394 °R
[85]:
res_ise_2 = isentropic_solver("m", 1, gamma=gamma)
res_ise_2.show()
key quantity
----------------------------
m M 1.00000000
pr P / P0 0.52828179
dr rho / rho0 0.63393815
tr T / T0 0.83333333
prs P / P* 1.00000000
drs rho / rho* 1.00000000
trs T / T* 1.00000000
urs U / U* 1.00000000
ars A / A* 1.00000000
ma Mach Angle 90.00000000
pm Prandtl-Meyer 0.00000000
[86]:
T02 = (1 / res_ise_2["tr"]) * T2
T02
[86]:
5760.0 °R
\[\begin{split}\begin{aligned}
q &= C_{P} (T_{02} - T_{01}) \\
\Delta T_{0} &= T_{02} - T_{01}
\end{aligned}\end{split}\]
[87]:
T01 = T02 - Delta_T0
T01
[87]:
3578.408772000606 °R
Assuming \(T_{1} = 25°C\):
[88]:
T1 = Q_(25, ureg.degC).to("degR")
T1
[88]:
536.67 °R
[89]:
res_ise_1 = isentropic_solver("temperature", T1 / T01, gamma=gamma)
res_ise_1.show()
key quantity
----------------------------
m M 5.32343922
pr P / P0 0.00130635
dr rho / rho0 0.00871051
tr T / T0 0.14997448
prs P / P* 0.00247284
drs rho / rho* 0.01374031
trs T / T* 0.17996938
urs U / U* 2.25835186
ars A / A* 32.22640460
ma Mach Angle [deg] 10.82725089
pm Prandtl-Meyer [deg] 79.79470958
[90]:
M1 = res_ise_1["m"]
M1
[90]:
np.float64(5.323439216045449)
P 3.12
[91]:
D = 0.02 * m
L = 40 * m
M2 = 0.5
p2 = 1 * atm
T2 = 270 * K
f = 0.005
[92]:
res_ise_2 = isentropic_solver("m", M2)
res_ise_2.show()
key quantity
----------------------------
m M 0.50000000
pr P / P0 0.84301918
dr rho / rho0 0.88517013
tr T / T0 0.95238095
prs P / P* 1.59577558
drs rho / rho* 1.39630363
trs T / T* 1.14285714
urs U / U* 0.53452248
ars A / A* 1.33984375
ma Mach Angle nan
pm Prandtl-Meyer nan
[93]:
T02 = (1 / res_ise_2["tr"]) * T2
p02 = (1 / res_ise_2["pr"]) * p2
print(T02)
print(p02)
283.5 K
1.1862126380443982 atm
[94]:
res_fanno_2 = fanno_solver("m", M2)
res_fanno_2.show()
key quantity
---------------------
m M 0.50000000
prs P / P* 2.13808994
drs rho / rho* 1.87082869
trs T / T* 1.14285714
tprs P0 / P0* 1.33984375
urs U / U* 0.53452248
fps 4fL* / D 1.06906031
eps (s*-s) / R 0.29255300
[95]:
crit_fric_L2s = res_fanno_2["fps"]
crit_fric_L2s
[95]:
np.float64(1.0690603127182559)
[96]:
crit_fric_L = 4 * f * L / D
crit_fric_L
[96]:
40.0
[97]:
crit_fric_L1s = crit_fric_L2s + crit_fric_L
crit_fric_L1s
[97]:
41.06906031271826
[98]:
res_fanno_1 = fanno_solver("friction_sub", crit_fric_L1s)
res_fanno_1.show()
key quantity
---------------------
m M 0.12573061
prs P / P* 8.69889577
drs rho / rho* 7.27199877
trs T / T* 1.19621799
tprs P0 / P0* 4.64652171
urs U / U* 0.13751377
fps 4fL* / D 41.06906031
eps (s*-s) / R 1.53611892
[99]:
M1 = res_fanno_1["m"]
Ts = (1 / res_fanno_2["trs"]) * T2
ps = (1 / res_fanno_2["prs"]) * p2
T1 = res_fanno_1["trs"] * Ts
p1 = res_fanno_1["prs"] * ps
print(M1)
print(T1)
print(p1)
0.12573061452167186
282.60650069373264 K
4.0685359499124925 atm
P 3.13
[100]:
D = 0.2 * ft
L = 3 * ft
M1 = 2.5
p1 = 0.5 * atm
T1 = 520 * degR
f = 0.005
[101]:
res_fanno_1 = fanno_solver("m", M1, gamma=gamma)
res_fanno_1.show()
key quantity
---------------------
m M 2.50000000
prs P / P* 0.29211870
drs rho / rho* 0.54772256
trs T / T* 0.53333333
tprs P0 / P0* 2.63671875
urs U / U* 1.82574186
fps 4fL* / D 0.43197669
eps (s*-s) / R 0.96953525
[102]:
res_ise_1 = isentropic_solver("m", M1, gamma=gamma)
res_ise_1.show()
key quantity
----------------------------
m M 2.50000000
pr P / P0 0.05852766
dr rho / rho0 0.13168724
tr T / T0 0.44444444
prs P / P* 0.11078872
drs rho / rho* 0.20772885
trs T / T* 0.53333333
urs U / U* 1.82574186
ars A / A* 2.63671875
ma Mach Angle 23.57817848
pm Prandtl-Meyer 39.12356383
[103]:
p0 = (1 / res_ise_1["pr"]) * p1
p0
[103]:
8.542968750000004 atm
[104]:
Ts = (1 / res_fanno_1["trs"]) * T1
ps = (1 / res_fanno_1["prs"]) * p1
p0s = (1 / res_fanno_1["tprs"]) * p0
print(Ts)
print(ps)
print(p0s)
975.0 °R
1.7116329922036442 atm
3.24 atm
[105]:
crit_fric_L1s = res_fanno_1["fps"]
crit_fric_L1s
[105]:
np.float64(0.4319766894222311)
[106]:
crit_fric_L = 4 * f * L / D
crit_fric_L
[106]:
0.3
[107]:
crit_fric_L2s = crit_fric_L1s - crit_fric_L
crit_fric_L2s
[107]:
0.13197668942223112
[108]:
res_fanno_2 = fanno_solver("friction_super", crit_fric_L2s, gamma=gamma)
res_fanno_2.show()
key quantity
---------------------
m M 1.48883667
prs P / P* 0.61243670
drs rho / rho* 0.73662199
trs T / T* 0.83141247
tprs P0 / P0* 1.16870081
urs U / U* 1.35754840
fps 4fL* / D 0.13197669
eps (s*-s) / R 0.15589271
[109]:
M2 = res_fanno_2["m"]
T2 = res_fanno_2["trs"] * Ts
p2 = res_fanno_2["prs"] * ps
p02 = res_fanno_2["tprs"] * p0s
print(M2)
print(T2)
print(p2)
print(p02)
1.48883666686136
810.6271562243851 °R
1.0482668680069764 atm
3.786590625873568 atm
P 3.14
[110]:
D = 4 * ureg.inch
p01_p02 = 10
p01 = 100 * atm
f = 0.005
[111]:
p02 = (1 / p01_p02) * p01
p02
[111]:
10.0 atm
For chocking to happen:
\[\frac{P_{02}}{P_{0}^{*}} = 1\]
[112]:
p0s = p02
p0s
[112]:
10.0 atm
[113]:
p01_p0s = p01 / p0s
p01_p0s
[113]:
10.0
[114]:
res_fanno_1 = fanno_solver("total_pressure_sub", p01_p0s)
res_fanno_1.show()
key quantity
---------------------
m M 0.05798720
prs P / P* 18.88480389
drs rho / rho* 15.74791999
trs T / T* 1.19919354
tprs P0 / P0* 10.00000000
urs U / U* 0.06350045
fps 4fL* / D 206.98591486
eps (s*-s) / R 2.30258509
[115]:
crit_fric_L1s = res_fanno_1["fps"]
crit_fric_L2s = 0
crit_fric_L = crit_fric_L1s - crit_fric_L2s
crit_fric_L
[115]:
np.float64(206.98591485503835)
[116]:
L = crit_fric_L * D / (4 * f)
L
[116]:
41397.18297100767 inch
[117]:
L.to("m")
[117]:
1051.4884474635949 meter
P 3.16
[118]:
M1 = 2.5
\[\begin{split}\begin{aligned}
\frac{3}{10} h_{1} &= q = c_{p} * (T_{02} - T_{01}) \\
\frac{3}{10} c_{p} T_{1} &= c_{p} * (T_{02} - T_{01}) \\
\frac{3}{10} &= \frac{T_{02} - T_{01}}{T_{1}} \\
\frac{3}{10} &= \frac{T_{02}}{T_{1}} - \frac{T_{01}}{T_{1}} \\
\end{aligned}\end{split}\]
[119]:
res_ise_1 = isentropic_solver("m", M1)
res_ise_1.show()
key quantity
----------------------------
m M 2.50000000
pr P / P0 0.05852766
dr rho / rho0 0.13168724
tr T / T0 0.44444444
prs P / P* 0.11078872
drs rho / rho* 0.20772885
trs T / T* 0.53333333
urs U / U* 1.82574186
ars A / A* 2.63671875
ma Mach Angle 23.57817848
pm Prandtl-Meyer 39.12356383
[120]:
T01_T1 = (1 / res_ise_1["tr"])
T01_T1
[120]:
np.float64(2.25)
[121]:
T02_T1 = T01_T1 + 0.3
T02_T1
[121]:
np.float64(2.55)
Assuming \(T_{1} = 25°C\):
[122]:
T1 = Q_(25, ureg.degC).to("K")
T01 = T01_T1 * T1
T02 = T02_T1 * T1
print(T01)
print(T02)
670.8375 K
760.2824999999999 K
[123]:
res_ray_1 = rayleigh_solver("m", M1)
res_ray_1.show()
key quantity
---------------------
m M 2.50000000
prs P / P* 0.24615385
drs rho / rho* 0.65000000
trs T / T* 0.37869822
tprs P0 / P0* 2.22183129
ttrs T0 / T0* 0.71005917
urs U / U* 1.53846154
eps (s*-s) / R 1.99675616
[124]:
T0s = (1 / res_ray_1["ttrs"]) * T01
T0s
[124]:
944.7628125 K
[125]:
T02_T0s = T02 / T0s
T02_T0s
[125]:
0.8047337278106508
[126]:
res_ray_2 = rayleigh_solver("total_temperature_super", T02_T0s)
res_ray_2.show()
key quantity
---------------------
m M 1.94447058
prs P / P* 0.38135479
drs rho / rho* 0.69353445
trs T / T* 0.54987145
tprs P0 / P0* 1.44609837
ttrs T0 / T0* 0.80473373
urs U / U* 1.44188943
eps (s*-s) / R 1.12922255
[127]:
M2 = res_ray_2["m"]
M2
[127]:
np.float64(1.94447057632683)
[ ]: