Chapter 7

[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
N = ureg.N
Q_ = ureg.Quantity

P 7.1

[2]:

from pygasflow.atd.avf.thickness_fp import (
    deltas_lam_ic,
    deltas_tur_ic,
)
Re_u = 1e06 / m
x = 1 * m

[3]:

Re = Re_u * x
Re

[3]:

1000000.0

[4]:

res1_lam = deltas_lam_ic(x, Re)
res1_lam.show()

key      quantity
----------------------
delta    δ [m]             0.00500000
delta_1  δ_1 [m]           0.00172080
delta_2  δ_2 [m]           0.00066410
H12      H12               2.59117603

[5]:

res1_tur = deltas_tur_ic(x, Re)
res1_tur.show()

key       quantity
-----------------------
delta     δ [m]             0.02334542
delta_vs  δ_vs [m]          0.00011569
delta_sc  δ_sc [m]          0.00053538
delta_1   δ_1 [m]           0.00292133
delta_2   δ_2 [m]           0.00227145
H12       H12               1.28611111

P 7.2

[6]:

x = 10 * m
Re = Re_u * x

[7]:

res2_lam = deltas_lam_ic(x, Re)
res2_lam.show()

key      quantity
----------------------
delta    δ [m]             0.01581139
delta_1  δ_1 [m]           0.00544165
delta_2  δ_2 [m]           0.00210007
H12      H12               2.59117603

[8]:

res2_tur = deltas_tur_ic(x, Re)
res2_tur.show()

key       quantity
-----------------------
delta     δ [m]             0.14729965
delta_vs  δ_vs [m]          0.00014565
delta_sc  δ_sc [m]          0.00084852
delta_1   δ_1 [m]           0.01843236
delta_2   δ_2 [m]           0.01433186
H12       H12               1.28611111

P 7.3

(a)

[9]:

res1_lam["delta"] / res2_lam["delta"]

[9]:

0.31622776601683794

[10]:

res1_lam["delta_1"] / res2_lam["delta_1"]

[10]:

0.31622776601683794

[11]:

res1_lam["delta_2"] / res2_lam["delta_2"]

[11]:

0.316227766016838

[12]:

res1_tur["delta"] / res2_tur["delta"]

[12]:

0.15848931924611137

[13]:

res1_tur["delta_1"] / res2_tur["delta_1"]

[13]:

0.15848931924611134

[14]:

res1_tur["delta_2"] / res2_tur["delta_2"]

[14]:

0.15848931924611134

[15]:

res1_tur["delta_sc"] / res2_tur["delta_sc"]

[15]:

0.6309573444801934

[16]:

res1_tur["delta_vs"] / res2_tur["delta_vs"]

[16]:

0.7943282347242816

The thickness of the boundary layer grows as the flow progresses over the flat plate.

(b)

[17]:

import matplotlib.pyplot as plt

x = np.linspace(1e-03, 20, 100) * m
Re = Re_u * x
r1 = deltas_lam_ic(x, Re)
# discard units
r1 = {k: v.magnitude if isinstance(v, Q_) else v for k, v in r1.items()}
x = x.magnitude

fig, ax = plt.subplots()
ax.plot(x, r1["delta"], label=r"$\delta$")
ax.plot(x, r1["delta_1"], label=r"$\delta_{1}$")
ax.plot(x, r1["delta_2"], label=r"$\delta_{2}$")
ax.legend()
ax.set_xlabel("x [m]")
ax.set_ylabel("Boundary-layer Thickness [m]")
ax.set_title("Incompressible laminar flat plate")
plt.show()

[18]:

x = np.linspace(1e-03, 20, 100) * m
Re = Re_u * x
r2 = deltas_tur_ic(x, Re)
# discard units
r2 = {k: v.magnitude if isinstance(v, Q_) else v for k, v in r2.items()}
x = x.magnitude

fig, (ax1, ax2) = plt.subplots(2, 1, sharex=True)
ax1.plot(x, r2["delta"], label=r"$\delta$")
ax1.plot(x, r2["delta_1"], label=r"$\delta_{1}$")
ax1.plot(x, r2["delta_2"], label=r"$\delta_{2}$")
ax1.legend()
ax2.plot(x, r2["delta_vs"], label=r"$\delta_{vs}$")
ax2.plot(x, r2["delta_sc"], label=r"$\delta_{sc}$")
ax2.legend()
ax2.set_xlabel("x [m]")
ax1.set_ylabel("Boundary-layer Thickness")
ax2.set_ylabel("Boundary-layer Thickness")
ax1.set_title("Incompressible turbulent flat plate")
plt.show()

[19]:

import matplotlib.cm as cm
from itertools import cycle

colors = cycle(cm.tab10.colors)
fig, ax = plt.subplots()
c = next(colors)
ax.plot(x, r1["delta"], color=c, label=r"$\delta_{lam}$")
ax.plot(x, r2["delta"], '--', color=c, label=r"$\delta_{tur}$")
c = next(colors)
ax.plot(x, r1["delta_1"], color=c, label=r"$\delta_{1, lam}$")
ax.plot(x, r2["delta_1"], '--', color=c, label=r"$\delta_{1, tur}$")
c = next(colors)
ax.plot(x, r1["delta_2"], color=c, label=r"$\delta_{1, lam}$")
ax.plot(x, r2["delta_2"], '--', color=c, label=r"$\delta_{1, tur}$")
ax.legend()
ax.set_xlabel("x [m]")
ax.set_ylabel("Boundary-layer Thickness Incompressible")
ax.set_title("Laminar vs Turbulent")
ax.set_yscale("log")
plt.show()

P 7.4

The wall shear stress for a compressible flat plate is (eq. 7.145):

\[\begin{split}\begin{aligned} \tau_{w} &= C \frac{\rho_{\infty} u_{\infty}^{2}}{\left(Re_{\infty, x}\right)^{n}} \left(\frac{T^{*}}{T_{\infty}}\right)^{n (1 + \omega) - 1} \\ &= C \frac{\rho_{\infty} u_{\infty}^{2}}{\left(Re_{\infty}^{u}\right)^{n}} \frac{1}{x^{n}} \left(\frac{T^{*}}{T_{\infty}}\right)^{n (1 + \omega) - 1} \end{aligned}\end{split}\]

The friction drag is given by:

\[\begin{split}\begin{aligned} D_{f} &= \int_{A} \tau_{w} dA = \int_{0}^{L_{ref}} \tau_{w} b \, dx = \frac{A_{ref}}{L_{ref}} \int_{0}^{L_{ref}} \tau_{w} \, dx \\ &= \frac{A_{ref}}{L_{ref}} C \frac{\rho_{\infty} u_{\infty}^{2}}{\left(Re_{\infty}^{u}\right)^{n}} \left(\frac{T^{*}}{T_{\infty}}\right)^{n (1 + \omega) - 1} \int_{0}^{L_{ref}} \frac{1}{x^{n}} \, dx \\ &= \frac{A_{ref}}{L_{ref}} C \frac{\rho_{\infty} u_{\infty}^{2}}{\left(Re_{\infty}^{u}\right)^{n}} \left(\frac{T^{*}}{T_{\infty}}\right)^{n (1 + \omega) - 1} \frac{L_{ref}^{1 - n}}{1 - n} \\ &= \frac{A_{ref}}{L_{ref}} C \frac{\rho_{\infty} u_{\infty}^{2}}{\left(Re_{\infty, L}\right)^{n}} \left(\frac{T^{*}}{T_{\infty}}\right)^{n (1 + \omega) - 1} \frac{L_{ref}}{1 - n} \\ &= A_{ref} \frac{C}{1 - n} \frac{\rho_{\infty} u_{\infty}^{2}}{\left(Re_{\infty, L}\right)^{n}} \left(\frac{T^{*}}{T_{\infty}}\right)^{n (1 + \omega) - 1} \\ &= \frac{1}{2} \rho_{\infty} u_{\infty}^{2} C_{D,f} A_{ref} \end{aligned}\end{split}\]

From which the average friction coefficient can be determined:

\[C_{D, f} = \frac{2C}{1 - n} \frac{1}{\left(Re_{\infty, L}\right)^{n}} \left(\frac{T^{*}}{T_{\infty}}\right)^{n (1 + \omega) - 1}\]

P 7.5

[20]:

from pygasflow.atd.avf.wall_shear_stress_fp import wss_c

[21]:

M_inf = 6
H = 30e03 * m
A_ref = 1860 * m**2
L_ref = 80 * m
gamma = 1.4
R = 287.05 * J / (kg * K)
Pr_star = 0.74
Tw1 = 1000 * K
Tw2 = 2000 * K

[22]:

atmosphere = Atmosphere(H.magnitude)
T_inf = atmosphere.temperature[0] * K
rho_inf = atmosphere.density[0] * kg / m**3
T_inf, rho_inf

[22]:

(<Quantity(226.509084, 'kelvin')>,
 <Quantity(0.0184101009, 'kilogram / meter ** 3')>)

[23]:

a_inf = sound_speed(gamma, R, T_inf).to("m/s")
v_inf = M_inf * a_inf
v_inf

[23]:

1810.2429106370967 m/s

[24]:

mu = viscosity_air_power_law(T_inf)
mu

[24]:

1.5765141243692437×10^-5 kg/(m s)

[25]:

Re_infL = Reynolds(rho_inf, v_inf, mu, L_ref)
Re_infL

[25]:

169116173.74147075

[26]:

q_inf = 0.5 * rho_inf * v_inf**2
q_inf

[26]:

30164.76059774928 kg/(m s²)

[27]:

rs = recovery_factor(Pr_star, laminar=True)
rs

[27]:

$\displaystyle 0.860232526704263$

[28]:

T_star1 = reference_temperature(T_inf, Tw1, M_inf, rs)
T_star1

[28]:

921.8977042109084 K

[29]:

# Compute the wall shear stress and friction coefficients for a
# compressible flat plate
res1 = wss_c(rho_inf, v_inf, Re_infL, T_star1 / T_inf, laminar=True)
res1.show()

key     quantity
------------------------------
tau_w   τw [kg / m / s ** 2]       1.20474712
cf      cf []                      0.00003994
CDf     CDf []                     0.00007988

[30]:

CDf1 = res1["CDf"]
D1 = 2 * q_inf * CDf1 * A_ref
D1 = D1.to("N")
D1

[30]:

8963.318554033402 N

[31]:

T_star2 = reference_temperature(T_inf, Tw2, M_inf, rs)
res2 = wss_c(rho_inf, v_inf, Re_infL, T_star2 / T_inf, laminar=True)
CDf2 = res2["CDf"]
D2 = 2 * q_inf * CDf2 * A_ref
D2 = D2.to("N")
D2

[31]:

8308.762554326322 N

P 7.6

[32]:

res1 = wss_c(rho_inf, v_inf, Re_infL, T_star1 / T_inf, laminar=False)
CDf1 = res1["CDf"]
D1 = 2 * q_inf * CDf1 * A_ref
D1 = D1.to("N")
D1

[32]:

73317.46654858341 N

[33]:

T_star2 = reference_temperature(T_inf, Tw2, M_inf, rs)
res2 = wss_c(rho_inf, v_inf, Re_infL, T_star2 / T_inf, laminar=False)
CDf2 = res2["CDf"]
D2 = 2 * q_inf * CDf2 * A_ref
D2 = D2.to("N")
D2

[33]:

54843.23756404523 N

P 7.8

[34]:

M_inf = 6
H = 30e03 * m
alpha = 6 * deg
A_ref = 1860 * m**2
L_ref = 80 * m
gamma = 1.4
Pr_star = 0.74

# windward side
M_w = 5.182
Re_wu = 3.123e06 / m
u_w = 1773.63 * m / s
rho_w = 0.0327 * kg / m**3
T_w = 291.52 * K
T_ww = 1000 * K

# leeward side
M_l = 6.997
Re_lu = 1.294e06 / m
u_l = 1840.01 * m / s
rho_l = 0.00927 * kg / m**3
T_l = 172.12 * K
T_wl = 800 * K

[35]:

atmosphere = Atmosphere(H.magnitude)
T_inf = atmosphere.temperature[0] * K
rho_inf = atmosphere.density[0] * kg / m**3
T_inf, rho_inf

[35]:

(<Quantity(226.509084, 'kelvin')>,
 <Quantity(0.0184101009, 'kilogram / meter ** 3')>)

[36]:

rs = recovery_factor(Pr_star, laminar=False)
Tr_w = recovery_temperature(T_w, M_w, rs, gamma)
Tr_l = recovery_temperature(T_l, M_l, rs, gamma)
Tr_w, Tr_l

[36]:

(<Quantity(1707.65202, 'kelvin')>, <Quantity(1696.50847, 'kelvin')>)

[37]:

Ts_w = reference_temperature(T_w, T_ww, Tr=Tr_w)
Ts_l = reference_temperature(T_l, T_wl, Tr=Tr_l)
Ts_w, Ts_l

[37]:

(<Quantity(957.309044, 'kelvin')>, <Quantity(821.425463, 'kelvin')>)

[38]:

res_w = wss_c(rho_w, u_w, Re_wu * L_ref, Ts_w / T_w, laminar=False)
res_w.show()

key     quantity
------------------------------
tau_w   τw [kg / m / s ** 2]      28.71150588
cf      cf []                      0.00055823
CDf     CDf []                     0.00069779

[39]:

res_l = wss_c(rho_l, u_l, Re_lu * L_ref, Ts_l / T_l, laminar=False)
res_l.show()

key     quantity
------------------------------
tau_w   τw [kg / m / s ** 2]       8.13304710
cf      cf []                      0.00051828
CDf     CDf []                     0.00064785

[40]:

D_w = 0.5 * rho_w * u_w**2 * res_w["CDf"] * A_ref
D_w = D_w.to("N")
D_w

[40]:

66754.25117310249 N

[41]:

D_l = 0.5 * rho_l * u_l**2 * res_l["CDf"] * A_ref
D_l = D_l.to("N")
D_l

[41]:

18909.33450540256 N

[42]:

D = D_w + D_l
Df = D * np.cos(alpha.to("radians"))
Df

[42]:

85194.3115930293 N

[43]:

Lf = -D * np.sin(alpha.to("radians"))
Lf

[43]:

-8954.2829689711 N

P 7.9

[44]:

L_inviscid = 1219344.1612489037 * N
D_inviscid = 128158.2355937337 * N

[45]:

L = Lf + L_inviscid
L

[45]:

1210389.8782799325 N

[46]:

D = Df + D_inviscid
D

[46]:

213352.547186763 N

[47]:

L / D

[47]:

5.673191598787851

[ ]: