Chapter 6
[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 6.1
[2]:
Tt = 1500 * K
gamma = 1.4
Cp = 1004 * J / kg / K
\[h_{t} = C_{p} T_{t} = C_{p} T + \frac{u^{2}}{2}\]
Maximum velocity happens when \(T = 0\):
[3]:
ht = Cp * Tt
ht
[3]:
1506000.0 J/kg
[4]:
u_max = np.sqrt(2 * ht).to("m/s")
u_max
[4]:
1735.511451993331 m/s
P 6.2
[5]:
# eq 6.19
a_star = np.sqrt(2 * (gamma - 1) / (gamma + 1) * Cp * Tt).to("m/s")
a_star
[5]:
708.519583356734 m/s
P 6.3
[6]:
gamma = 1.4
M_inf = [1e-08, 0.01, 0.1, 0.5, 1]
[7]:
cp = pressure_coefficient(M_inf, gamma=gamma, stagnation=True)
cp
[7]:
array([1. , 1.000025 , 1.0025025 , 1.06407222, 1.27561308])
P 6.4
[8]:
gamma = 1.4
M_inf = [1, 5, 10, 1e05]
cp = pressure_coefficient(M_inf, gamma=gamma, stagnation=True)
cp
[8]:
array([1.27561308, 1.80876996, 1.83167098, 1.83937105])
P 6.5
\[\Delta_{0} = \varepsilon \frac{R_{s}}{1 - \varepsilon}\]
[9]:
from pygasflow.shockwave import density_ratio
Rb = 1 * m
gamma = 1.4
M_inf = [5, 10, 1e10]
[10]:
eps = 1 / density_ratio(M_inf, gamma)
Rs = Rb
delta_0 = eps * Rs / (1 - eps)
delta_0
[10]:
| Magnitude | [0.24999999999999994 0.21212121212121204 0.19999999999999996] |
|---|---|
| Units | m |
P 6.6
[11]:
gamma = 1.2
eps = 1 / density_ratio(M_inf, gamma)
Rs = Rb
delta_0 = eps * Rs / (1 - eps)
delta_0
[11]:
| Magnitude | [0.1458333333333333 0.11111111111111106 0.09999999999999996] |
|---|---|
| Units | m |
P 6.7
[12]:
gamma = 1+1e-05
eps = 1 / density_ratio(M_inf, gamma)
Rs = Rb
delta_0 = eps * Rs / (1 - eps)
delta_0
[12]:
| Magnitude | [0.04167187500000004 0.01010606060606064 5.000000000032766e-06] |
|---|---|
| Units | m |
P 6.8
[13]:
gamma = 1.4
R = 287.05 * J / kg / K
M_inf = 6
H = 30 * km
alpha = 6 * deg
A_ref = 1860 * m**2
L_ref = 40 * m
[14]:
atmosphere = Atmosphere(H.to("m").magnitude)
T_inf = atmosphere.temperature[0] * K
rho_inf = atmosphere.density[0] * kg / m**3
T_inf, rho_inf
[14]:
(<Quantity(226.509084, 'kelvin')>,
<Quantity(0.0184101009, 'kilogram / meter ** 3')>)
[15]:
a_inf = sound_speed(gamma, R, T_inf).to("m/s")
a_inf
[15]:
301.70715177284944 m/s
[16]:
v_inf = a_inf * M_inf
v_inf
[16]:
1810.2429106370967 m/s
[17]:
alpha = alpha.to("radians")
cL = 2 * np.sin(alpha)**2 * np.cos(alpha)
cD = 2 * np.sin(alpha)**3
[18]:
L = 0.5 * rho_inf * v_inf**2 * A_ref * cL
L = L.to("N")
L
[18]:
1219344.1612489037 N
[19]:
D = 0.5 * rho_inf * v_inf**2 * A_ref * cD
D = D.to("N")
D
[19]:
128158.2355937337 N
[20]:
L / D
[20]:
9.514364454222585
[21]:
M = L_ref * (L * np.cos(alpha) + D * np.sin(alpha))
M
[21]:
49042425.99092816 m N
P 6.9
[22]:
gamma = 1.4
R = 287.05 * J / kg / K
M_inf = 6
H = 40 * km
alpha = 6 * deg
A_ref = 1860 * m**2
L_ref = 40 * m
[23]:
atmosphere = Atmosphere(H.to("m").magnitude)
T_inf = atmosphere.temperature[0] * K
rho_inf = atmosphere.density[0] * kg / m**3
T_inf, rho_inf
[23]:
(<Quantity(250.349646, 'kelvin')>,
<Quantity(0.00399565628, 'kilogram / meter ** 3')>)
[24]:
def compute(alpha, A_ref, rho_inf, T_inf):
a_inf = sound_speed(gamma, R, T_inf).to("m/s")
v_inf = a_inf * M_inf
alpha = alpha.to("radians")
cL = 2 * np.sin(alpha)**2 * np.cos(alpha)
cD = 2 * np.sin(alpha)**3
L = 0.5 * rho_inf * v_inf**2 * A_ref * cL
L = L.to("N")
D = 0.5 * rho_inf * v_inf**2 * A_ref * cD
D = D.to("N")
M = L_ref * (L * np.cos(alpha) + D * np.sin(alpha))
return L, D, M
L, D, M = compute(alpha, A_ref, rho_inf, T_inf)
[25]:
L
[25]:
292495.7871300942 N
[26]:
D
[26]:
30742.54602474064 N
[27]:
L / D
[27]:
9.514364454222585
[28]:
M
[28]:
11764277.4278703 m N
P 6.10
(a) - increase angle of attack
[29]:
L, D, M = compute(10 * deg, A_ref, rho_inf, T_inf)
[30]:
L
[30]:
799333.5724989552 N
[31]:
D
[31]:
140944.07541765162 N
[32]:
L / D
[32]:
5.671281819617711
[33]:
M
[33]:
32466583.2515657 m N
The lift increased, but so did the drag. So much so that \(L/D\) decreased rapidly.
(b) - increase surface area
[34]:
L, D, M = compute(alpha, 1.5 * A_ref, rho_inf, T_inf)
[35]:
L
[35]:
438743.6806951413 N
[36]:
D
[36]:
46113.81903711096 N
[37]:
L/D
[37]:
9.514364454222585
[38]:
M
[38]:
17646416.14180545 m N
[ ]: