Chapter 3

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

P 3.1

ht, hinf, vinf, Minf = symbols(r"h_t, h_\infty, v_\infty, M_\infty")
h = Eqn(ht, hinf + vinf**2 / 2)
h

$\displaystyle h_{t} = h_{\infty} + \frac{v_{\infty}^{2}}{2}$

cv, cp, gamma, R, Tt, Tinf, ainf = symbols(r"c_v, c_p, gamma, R, T_t, T_\infty, a_\infty")
e1 = Eqn(gamma, cp / cv)
e2 = Eqn(cv, R / (gamma - 1))

h = h.subs({ht: cp * Tt, hinf: cp * Tinf})
h

$\displaystyle T_{t} c_{p} = T_{\infty} c_{p} + \frac{v_{\infty}^{2}}{2}$

h = (h / cp).expand()
h

$\displaystyle T_{t} = T_{\infty} + \frac{v_{\infty}^{2}}{2 c_{p}}$

e3 = solve(e1.subs(e2), cp)[0]
e3

$\displaystyle c_{p} = \frac{R \gamma}{\gamma - 1}$

h = h.subs(e3)
h

$\displaystyle T_{t} = T_{\infty} + \frac{v_{\infty}^{2} \left(\gamma - 1\right)}{2 R \gamma}$

h = h.subs(vinf, Minf * ainf)
h

$\displaystyle T_{t} = \frac{M_{\infty}^{2} a_{\infty}^{2} \left(\gamma - 1\right)}{2 R \gamma} + T_{\infty}$

h = h.subs(ainf, sqrt(gamma * R * Tinf))
h

$\displaystyle T_{t} = \frac{M_{\infty}^{2} T_{\infty} \left(\gamma - 1\right)}{2} + T_{\infty}$

h = h.collect(Tinf)
h

$\displaystyle T_{t} = T_{\infty} \left(\frac{M_{\infty}^{2} \left(\gamma - 1\right)}{2} + 1\right)$

P 3.2

alpha = 0 * deg
vinf = 1 * km / s
H = 30 * km
gamma = 1.4
R = 287.05 * J / (kg * K)

atmosphere = Atmosphere(H.to("m").magnitude)
Tinf = atmosphere.temperature[0] * K
Tinf

226.50908361133006 K

ainf = sound_speed(gamma, R, Tinf).to("km / s")
ainf

0.30170715177284946 km/s

Minf = vinf / ainf
Minf

3.3144723090716925

flow = isentropic_solver("m", Minf)
flow.show()

key     quantity
----------------------------
m       M                        3.31447231
pr      P / P0                   0.01711322
dr      rho / rho0               0.05471344
tr      T / T0                   0.31277902
prs     P / P*                   0.03239411
drs     rho / rho*               0.08630723
trs     T / T*                   0.37533483
urs     U / U*                   2.03059741
ars     A / A*                   5.70596370
ma      Mach Angle [deg]        17.56016750
pm      Prandtl-Meyer [deg]     55.46993616

T_Tt = flow["tr"]
Tt = (1 / T_Tt) * Tinf
Tt

724.1824605740292 K

Compute the recovery factors:

Pr = Prandtl(gamma)
r_lam = recovery_factor(Pr, laminar=True)
r_tur = recovery_factor(Pr, laminar=False)
print("Recovery factor for laminar boundary layer:", r_lam)
print("Recovery factor for turbulent boundary layer:", r_tur)

Recovery factor for laminar boundary layer: 0.8583950752789521
Recovery factor for turbulent boundary layer: 0.9032157004286028

Compute the recovery temperatures:

Tr_lam = recovery_temperature(Tinf, Minf, r_lam)
Tr_lam

653.7094594935564 K

Tr_tur = recovery_temperature(Tinf, Minf, r_tur)
Tr_tur

676.0154913693625 K

P 3.3

Tra, Tr, eps, Re, L, n, C, x_L = symbols(r"T_ra, T_r, \varepsilon, Re, L, n, C, (x/L)", positive=True)
e0 = Eqn(
    Tra**4,
    Tra**(2 * n - 1) / eps * Re**(1 - n) / x_L**n / L * Tr * (1 - Tra / Tr)
)
e0

$\displaystyle T_{ra}^{4} = \frac{(x/L)^{- n} Re^{1 - n} T_{r} T_{ra}^{2 n - 1} \left(1 - \frac{T_{ra}}{T_{r}}\right)}{L \varepsilon}$

Assuming \(T_{ra} << T_{r}\):

e = e0.subs(Tra / Tr, 0)
e

$\displaystyle T_{ra}^{4} = \frac{(x/L)^{- n} Re^{1 - n} T_{r} T_{ra}^{2 n - 1}}{L \varepsilon}$

Bringing \(T_{ra}^{2 n - 1}\) to the LHS:

e = (e / Tra**(2 * n - 1)).apply(sp.powsimp)
e

$\displaystyle T_{ra}^{5 - 2 n} = \frac{(x/L)^{- n} Re^{1 - n} T_{r}}{L \varepsilon}$

e = e.apply(lambda t: sp.root(t, 5-2*n).expand().powsimp())
e

$\displaystyle T_{ra} = (x/L)^{- \frac{n}{5 - 2 n}} Re^{\frac{1 - n}{5 - 2 n}} \left(\frac{T_{r}}{L \varepsilon}\right)^{\frac{1}{5 - 2 n}}$

toe = table_of_expressions(e.rhs)

idx	args
0	\((x/L)^{- \frac{n}{5 - 2 n}}\)
1	\(Re^{\frac{1 - n}{5 - 2 n}}\)
2	\(\left(\frac{T_{r}}{L \varepsilon}\right)^{\frac{1}{5 - 2 n}}\)

C = sp.symbols("C")
c = Eqn(C, sp.Mul(*toe[1:]))
c

$\displaystyle C = Re^{\frac{1 - n}{5 - 2 n}} \left(\frac{T_{r}}{L \varepsilon}\right)^{\frac{1}{5 - 2 n}}$

e = e.subs(c.swap)
e

$\displaystyle T_{ra} = (x/L)^{- \frac{n}{5 - 2 n}} C$

Let’s solve for \(C\):

C_eq = solve(e, C)[0]
C_eq

$\displaystyle C = (x/L)^{\frac{n}{5 - 2 n}} T_{ra}$

And create the numerical functions so that they can be evaluated with pint quantities:

C_func = sp.lambdify([x_L, Tra, n], C_eq.rhs)
T_func = sp.lambdify([x_L, C, n], e.rhs)

For the laminar case, \(n=0.5\):

n_lam = 0.5
x_L1 = 0.1
Tw1 = 1100 * K # from Fig. 3.3
C_val = C_func(x_L1, Tw1, n_lam)
C_val

824.8836302657014 K

x_L2 = 0.75
Tw2 = 880 * K # from Fig. 3.3
Tw2_comp = T_func(x_L2, C_val, n_lam)
Tw2, Tw2_comp

(<Quantity(880, 'kelvin')>, <Quantity(855.086455, 'kelvin')>)

Tw2 - Tw2_comp

24.91354462818333 K

It clearly underestimates the result, but not by much. Maybe it has something to do with the assumption that \(T_{ra} << T_{r}\).

P 3.4

Tw1 = 656 * K
C_val = C_func(x_L1, Tw1, n_lam)
C_val

491.93060132209104 K

Tw2 = 522 * K
Tw2_comp = T_func(x_L2, C_val, n_lam)
Tw2_comp

509.94246793082885 K

Tw2 - Tw2_comp

12.057532069171145 K

P 3.5

n_tur = 0.2

Tw1 = 874 * K
C_val = C_func(x_L1, Tw1, n_tur)
C_val

790.7390229813286 K

Tw2 = 820 * K
Tw2_comp = T_func(x_L2, C_val, n_tur)
Tw2_comp

800.6916338680519 K

Tw2 - Tw2_comp

19.308366131948105 K

P 3.6

n_tur = 0.1

Tw1 = 874 * K
C_val = C_func(x_L1, Tw1, n_tur)
C_val

833.0634859472043 K

Tw2 = 820 * K
Tw2_comp = T_func(x_L2, C_val, n_tur)
Tw2_comp

838.0713410649367 K

Tw2 - Tw2_comp

-18.071341064936746 K

P 3.7

H = 60.65 * km
Minf = 17
x_L = 0.5
gamma = 1.3

atmosphere = Atmosphere(H.to("m").magnitude)
Tinf = atmosphere.temperature[0] * K
Tinf

245.2349423942596 K

e0

$\displaystyle T_{ra}^{4} = \frac{(x/L)^{- n} Re^{1 - n} T_{r} T_{ra}^{2 n - 1} \left(1 - \frac{T_{ra}}{T_{r}}\right)}{L \varepsilon}$

Assuming \(T_{ra} << T_{r}\):

e = e0.subs(Tra / Tr, 0)
e

$\displaystyle T_{ra}^{4} = \frac{(x/L)^{- n} Re^{1 - n} T_{r} T_{ra}^{2 n - 1}}{L \varepsilon}$

e = (e / Tra**(2*n - 1)).apply(sp.powsimp)
e

$\displaystyle T_{ra}^{5 - 2 n} = \frac{(x/L)^{- n} Re^{1 - n} T_{r}}{L \varepsilon}$

e = e.apply(sp.root, e.lhs.exp)
e

$\displaystyle T_{ra} = \left(\frac{(x/L)^{- n} Re^{1 - n} T_{r}}{L \varepsilon}\right)^{\frac{1}{5 - 2 n}}$

e = e.expand()
e

$\displaystyle T_{ra} = (x/L)^{- \frac{n}{5 - 2 n}} L^{- \frac{1}{5 - 2 n}} Re^{- \frac{n}{5 - 2 n} + \frac{1}{5 - 2 n}} T_{r}^{\frac{1}{5 - 2 n}} \varepsilon^{- \frac{1}{5 - 2 n}}$

toe = table_of_expressions(e.rhs)

idx	args
0	\((x/L)^{- \frac{n}{5 - 2 n}}\)
1	\(L^{- \frac{1}{5 - 2 n}}\)
2	\(Re^{- \frac{n}{5 - 2 n} + \frac{1}{5 - 2 n}}\)
3	\(T_{r}^{\frac{1}{5 - 2 n}}\)
4	\(\varepsilon^{- \frac{1}{5 - 2 n}}\)

c_eq = Eqn(C, e.rhs / toe[3])
c_eq

$\displaystyle C = (x/L)^{- \frac{n}{5 - 2 n}} L^{- \frac{1}{5 - 2 n}} Re^{- \frac{n}{5 - 2 n} + \frac{1}{5 - 2 n}} \varepsilon^{- \frac{1}{5 - 2 n}}$

e = e.subs(c_eq.swap)
e

$\displaystyle T_{ra} = C T_{r}^{\frac{1}{5 - 2 n}}$

Assume:

• \(Pr = 1\): typical order of magnitude for gases

• \(\gamma = 1.3\): take into account for chemical reacting gases.

Pr = 1
gamma = 1.3
r_lam = recovery_factor(Pr, laminar=True)
r_lam

$\displaystyle 1.0$

Tr_lam_1 = recovery_temperature(Tinf, Minf, r_lam, gamma)
Tr_lam_2 = recovery_temperature(Tinf, 15.7, r_lam, gamma)
Tr_lam_1, Tr_lam_2

(<Quantity(10876.1697, 'kelvin')>, <Quantity(9312.42909, 'kelvin')>)

Let’s compute the ratio \(T_{ra, 1} / T_{ra, 2}\):

\[\frac{T_{ra, 1}}{T_{ra, 2}} = \frac{C T_{r, 1}^{\frac{1}{5 - 2 n}}}{C T_{r, 2}^{\frac{1}{5 - 2 n}}} = \frac{T_{r, 1}^{\frac{1}{5 - 2 n}}}{T_{r, 2}^{\frac{1}{5 - 2 n}}}\]

Tr_expr = e.rhs / C
Tr_expr

$\displaystyle T_{r}^{\frac{1}{5 - 2 n}}$

Tr_func = sp.lambdify([Tr, n], Tr_expr)
n_lam = 0.5
Tra1_M17 = Tr_func(Tr_lam_1, n_lam)
Tra2_M15_7 = Tr_func(Tr_lam_2, n_lam)
Tra1_M17, Tra2_M15_7

(<Quantity(10.2121925, 'kelvin ** 0.25')>,
 <Quantity(9.82348858, 'kelvin ** 0.25')>)

ratio = Tra1_M17 / Tra2_M15_7
ratio

1.0395688295125671

From Figure 3.3, at \(M_{\infty} = 15.7\) and \(x/L = 0.5\):

Tw_M15_7 = 916 * K

Then, for \(M_{\infty} = 17\) at \(x / L = 0.5\):

Tw_M17 = ratio * Tw_M15_7
Tw_M17

952.2450478335115 K

P 3.8

Tr_lam_1 = recovery_temperature(Tinf, 14, r_lam, gamma)
Tr_lam_2 = recovery_temperature(Tinf, 15.7, r_lam, gamma)
Tr_lam_1, Tr_lam_2

(<Quantity(7455.14225, 'kelvin')>, <Quantity(9312.42909, 'kelvin')>)

Tra1_M14 = Tr_func(Tr_lam_1, n_lam)
Tra2_M15_7 = Tr_func(Tr_lam_2, n_lam)
Tra1_M14, Tra2_M15_7

(<Quantity(9.29210232, 'kelvin ** 0.25')>,
 <Quantity(9.82348858, 'kelvin ** 0.25')>)

ratio = Tra1_M14 / Tra2_M15_7
ratio

0.9459065634888746

Tw_M15_7 = 916 * K
Tw_M14 = ratio * Tw_M15_7
Tw_M14

866.4504121558091 K