Pressure distribution on the Apollo Command Module

In the following notebook we are quickly going to check the pressure distribution on the Apollo Command Module computed with the modified Newtonian theory, and compare it with laboratory data.

For more information check:

  • “Hypersonic Aerothermodynamics”, Chapter 6, by John J. Bertin

  • “The effect of protuberances, cavities and angle of attack on the wind-tunnel pressure and heat-transfer distribution for the Apollo Command Module”, John J. Bertin

Note: the results will not be commented here. Refer to the above source to better understand the differences between computed results and lab-data.

Apollo Command Module geometry

[1]:

import matplotlib.pyplot as plt
import numpy as np

R1 = 4.694
R2 = 0.196
R3 = 0.39
Rb = 1.956
# half-angle of the conical section
theta1 = 33
# angle of the intersection between arc of radius R1 with
# arc of radius R2
theta2 = 23.045

A quick drawing to visualize the profile:

[2]:

t = np.deg2rad(np.linspace(0, theta2, 100))
x1 = R1 * np.cos(t)
y1 = R1 * np.sin(t)
t = np.deg2rad(np.linspace(theta2, 90 + theta1, 100))
x2 = R2 * np.cos(t) + x1[-1] - R2 * np.cos(np.deg2rad(theta2))
y2 = R2 * np.sin(t) + y1[-1] - R2 * np.sin(np.deg2rad(theta2))
t = np.deg2rad(np.linspace(90 + theta1, 180, 100))
d = 1.91627
x3 = R3 * np.cos(t) + d
y3 = R3 * np.sin(t)

x = np.concatenate([x1, x2, x3])
x = np.concatenate([x, np.flip(x)])
y = np.concatenate([y1, y2, y3])
y = np.concatenate([y, -np.flip(y)])

plt.figure()
plt.plot(x, y)
plt.arrow(0, 0, R1 - 0.12, 0, width=0.0005, head_width=0.08)
plt.annotate("$R_{1}$", (2 * R1 / 3, 0.1))
plt.arrow(4.13937, 1.760, 0, R2, width=0.0005)
plt.arrow(4.13937, 1.760 + 3 * R2 + 0.1, 0, -2 * R2, width=0.0005, head_width=0.08)
plt.annotate("$R_{2}$", (4.13937 + 0.05, 1.760 + 2 * R2))
plt.arrow(4.13937, 0, 0, -Rb + 0.12, width=0.0005, head_width=0.08)
plt.annotate("$R_{b}$", (4.13937 + 0.05, -Rb / 2))

t = np.deg2rad(np.linspace(0, theta1))
x4 = 1.06934 + (2.2 - 1.06934) * np.cos(t)
y4 = (2.2 - 1.06934) * np.sin(t)
plt.plot(x4, y4, 'k:')
plt.annotate(r"$\theta_{1}$", (2.25, 0.35))
plt.plot([0, R1 * np.cos(np.deg2rad(theta2))], [0, R1 * np.sin(np.deg2rad(theta2))], 'k--', alpha=0.5)

t = np.deg2rad(np.linspace(0, theta2))
x5, y5 = 3.85 * np.cos(t), 3.85 * np.sin(t)
plt.plot(x5, y5, 'k:')
plt.annotate(r"$\theta_{2}$", (3.85, 0.75))

t = np.deg2rad(np.linspace(0, 5))
x6, y6 = 4.8 * np.cos(t), 4.8 * np.sin(t)
plt.plot(x6, y6, 'k')
plt.arrow(x6[-1], y6[-1], -0.005, 0.1, head_width=0.08)
plt.annotate("s", (4.85, 0.25))

plt.xlabel("x [m]")
plt.ylabel("y [m]")
plt.title("Apollo Command Module profile")
plt.axis('scaled')
plt.show()
../_images/examples_tut-10_4_0.png

Let’s now compute the length \(s\) and the associated local body angles \(\theta_{b}\):

[3]:

theta3 = 90 - theta2
# arc length of radius R1
sR1 = np.deg2rad(theta2) * R1
# total length up to the maximum body radius Rb
s = sR1 + np.deg2rad(theta3) * R2
# discretize the length
s_arr = np.linspace(0, s, 200)
# local body angle
theta_b = np.zeros_like(s_arr)
for i, v in enumerate(s_arr):
    if v <= sR1:
        theta_b[i] = np.pi / 2 - s_arr[i] / R1
    else:
        theta_b[i] = np.pi / 2 - ((s_arr[i] - sR1) / R2 + np.deg2rad(theta2))
s_arr = np.concatenate([np.flip(-s_arr), s_arr])
theta_b = np.concatenate([np.pi - np.flip(theta_b), theta_b])

Note:

  • \(s / R_{b} = 0.965\) defines the tangency point between the arc of radius \(R_{1}\) and the arc of radius \(R_{2}\).

  • \(s / R_{b} = 1.082\) corresponds to the maximum body radius \(R_{b}\).

0deg angle of attack

The experimental data is given here. The first column represents \(s / R_{B}\) while the second column is \(p_{s} / p_{t2}\).

[4]:

# free stream Mach number
Minf = 10.18

data_0deg = np.array([
    [-1.027237354085603, 0.1834808259587022],
    [-0.9478599221789883, 0.5775811209439529],
    [-0.7844357976653697, 0.8442477876106195],
    [-0.2334630350194553, 0.983480825958702],
    [0.03268482490272362, 0.9929203539823008],
    [0.2894941634241248, 0.966961651917404],
    [0.5649805447470819, 0.9244837758112094],
    [0.7003891050583655, 0.8418879056047198],
    [0.8357976653696495, 0.8135693215339233],
    [0.891828793774319, 0.7734513274336283],
    [0.943190661478599, 0.7097345132743362],
    [0.9992217898832685, 0.5256637168141592],
    [1.045914396887159, 0.17640117994100302],
    [1.0972762645914398, 0.018289085545722727]
])

Let’s now compute the pressure distribution with a modified newtonian flow model:

[5]:

from pygasflow.newton import modified_newtonian_pressure_ratio
ps_pt2 = modified_newtonian_pressure_ratio(Minf, theta_b)

plt.figure()
plt.plot(data_0deg[:, 0], data_0deg[:, 1], 's', label="lab data")
plt.plot(s_arr / Rb, ps_pt2, label="Modified Newtonian theory")
plt.xlabel(r"$s / R_{b}$")
plt.ylabel(r"$p_{s} / p_{t2}$")
plt.title(r"$M_{\infty} = 10.18 \qquad \alpha = 0°$")
plt.legend()
plt.show()
../_images/examples_tut-10_11_0.png

33deg angle of attack

The experimental data is given here. The first column represents \(s / R_{B}\) while the second column is \(p_{s} / p_{t2}\).

[6]:

data_33deg = np.array([
    [-1.0217821782178214, 0.06461538461538452],
    [-0.9623762376237623, 0.24615384615384595],
    [-0.7960396039603954, 0.36307692307692285],
    [-0.2554455445544548, 0.6030769230769228],
    [0.00594059405940639, 0.723076923076923],
    [0.27326732673267395, 0.8307692307692307],
    [0.5465346534653475, 0.9353846153846151],
    [0.6831683168316838, 0.9969230769230768],
    [0.8198019801980205, 1.006153846153846],
    [0.8613861386138622, 1.003076923076923],
    [0.9267326732673276, 0.9938461538461537],
    [0.9801980198019811, 0.9323076923076921],
    [1.0455445544554465, 0.6492307692307691],
    [1.0990099009901, 0.2338461538461536]
])

Let’s now compute the pressure distribution with a modified newtonian flow model:

[7]:

ps_pt2 = modified_newtonian_pressure_ratio(Minf, theta_b, alpha=np.deg2rad(33), beta=np.deg2rad(0))

plt.figure()
plt.plot(data_33deg[:, 0], data_33deg[:, 1], 's', label="lab data")
plt.plot(s_arr / Rb, ps_pt2, label="Modified Newtonian theory")
plt.xlabel(r"$s / R_{b}$")
plt.ylabel(r"$p_{s} / p_{t2}$")
plt.title(r"$M_{\infty} = 10.18 \qquad \alpha = 33°$")
plt.legend()
plt.show()
../_images/examples_tut-10_15_0.png

Location of Stagnation Point as a function of the angle of attack

First, let’s load a few experimental data:

[8]:

data1 = np.array([
    [0, -0.11709601873536712],
    [0.5274418604651163, 15.05854800936767],
    [0.5302325581395348, 20.257611241217788],
    [0.6697674418604653, 25.17564402810305],
    [0.6697674418604653, 31.21779859484777],
    [0.8037209302325581, 25.31615925058547],
    [0.8009302325581396, 28.266978922716618],
    [0.798139534883721, 30.093676814988285],
    [0.8009302325581396, 31.49882903981264],
    [0.8065116279069766, 33.18501170960187],
    [0.8483720930232559, 33.18501170960187],
    [0.8762790697674419, 33.18501170960187],
    [0.8483720930232559, 30.23419203747072],
    [0.9069767441860466, 39.929742388758775],
    [0.9683720930232558, 44.98829039812646],
    [0.9627906976744187, 50.46838407494145]
])

data2 = np.array([
    [0, -0.11709601873536712],
    [0.041860465116279055, 2.4121779859484747],
    [0.08372093023255811, 4.941451990632316],
    [0.10883720930232554, 6.065573770491795],
    [0.1590697674418604, 8.032786885245898],
    [0.22046511627906962, 9.578454332552681],
    [0.2790697674418603, 10.983606557377044],
    [0.32651162790697663, 12.107728337236523],
    [0.37953488372093025, 13.793911007025756],
    [0.43255813953488387, 15.761124121779858],
    [0.4800000000000002, 17.728337236533946],
    [0.521860465116279, 19.69555035128805],
    [0.5804651162790697, 22.646370023419195],
    [0.6418604651162791, 26.018735362997646],
    [0.6948837209302325, 29.39110070257611],
    [0.7479069767441859, 32.903981264637],
    [0.798139534883721, 36.557377049180324],
    [0.8483720930232559, 40.35128805620609],
])

Now, we are going to use the modified newtonian theory to estimate the location of the stagnation point:

[9]:

from scipy.optimize import minimize
Mfs = 10.18
beta = 0
gamma = 1.4

[10]:

def location_of_stagnation_point(s_arr, theta_b, theta):
    j = 0
    go_on = True
    t1, t2 = None, None
    s1, s2 = None, None
    while (j < len(theta_b) - 1) and go_on:
        if (theta <= theta_b[j]) and (theta >= theta_b[j + 1]):
            t1 = theta_b[j]
            t2 = theta_b[j + 1]
            s1 = s_arr[j]
            s2 = s_arr[j + 1]
            go_on = False
        j += 1
    if s1 == s2:
        return s1
    # linear interpolation: the finer the discretization the more
    # accurate the result
    t_percent = (theta - t1) / (t2 - t1)
    return s1 + t_percent * (s2 - s1)

All we have to do is to find the \(\theta_{b}\) value at which \(P_{s} / P_{t2}\) is maximum.

[11]:

# function to minimize
func = lambda theta_b, M, alpha=0, beta=0, gamma=1.4: 1.0 / modified_newtonian_pressure_ratio(M, theta_b, alpha=alpha, beta=beta, gamma=gamma)

alpha = np.deg2rad(np.linspace(0, 70, 100))
loc = np.zeros_like(alpha)
for i, a in enumerate(alpha):
    res = minimize(func, np.pi/2, bounds=[(0, np.pi)], args=(Mfs, a), tol=1e-10)
    loc[i] = location_of_stagnation_point(s_arr, theta_b, res.x[0])

plt.figure()
plt.plot(data1[:, 0], data1[:, 1], 's', label="lab data 1")
plt.plot(data2[:, 0], data2[:, 1], '--', label="lab data 2")
plt.plot(loc / Rb, np.rad2deg(alpha), label="Modified Newtonian theory")
plt.xlabel("$s / R_{b}$")
plt.ylabel(r"$\alpha$ [deg]")
plt.title("Location of the Stagnation Point")
plt.legend()
plt.show()
../_images/examples_tut-10_22_0.png

Again, the exaplanation of the result can be found in the aforementioned sources.

[ ]: