Compressible Flow & Aerodynamics

Question 1

(a) Area–velocity relation derivation (8 marks)

Continuity (steady, 1-D): $\rho A V = \text{const}$. Log-differentiate: $\frac{d\rho}{\rho} + \frac{dA}{A} + \frac{dV}{V} = 0 \quad (1)$ (2 marks)

Euler (momentum), inviscid, no body force: $dP = -\rho V\,dV$. (1 mark)

Isentropic sound speed: $a^2 = \dfrac{dP}{d\rho}\Big|_s \Rightarrow d\rho = \dfrac{dP}{a^2} = \dfrac{-\rho V\,dV}{a^2}$, so $\frac{d\rho}{\rho} = -\frac{V\,dV}{a^2} = -M^2\frac{dV}{V}\quad (2)$ (2 marks)

Substitute (2) into (1): $-M^2\frac{dV}{V} + \frac{dA}{A} + \frac{dV}{V} = 0 \Rightarrow \boxed{\frac{dA}{A} = (M^2-1)\frac{dV}{V}}$ (2 marks)

Physical interpretation: For $M<1$, $(M^2-1)<0$: area must decrease to accelerate ($dV>0 \Rightarrow dA<0$). For $M>1$, $(M^2-1)>0$: area must increase to accelerate. To pass from subsonic to supersonic, $dA$ must change sign, which requires $dA=0$ at $M=1$ — i.e. a throat where sonic conditions occur. (1 mark)

(b) Supersonic design exit conditions (10 marks)

Area–Mach relation: $\frac{A}{A^*} = \frac{1}{M}\left[\frac{2}{\gamma+1}\left(1+\frac{\gamma-1}{2}M^2\right)\right]^{\frac{\gamma+1}{2(\gamma-1)}} = 2.5$

Solving the supersonic root gives $M_e \approx 2.44$. (4 marks)

Isentropic ratios at $M_e = 2.44$:

• $\dfrac{T_0}{T_e} = 1 + 0.2 M_e^2 = 1 + 0.2(5.9536) = 2.1907 \Rightarrow T_e = 600/2.1907 = 273.9\ \text{K}$ (2 marks)
• $\dfrac{P_0}{P_e} = (T_0/T_e)^{3.5} = 2.1907^{3.5} = 15.53 \Rightarrow P_e = 10^6/15.53 = 64.4\ \text{kPa}$ (2 marks)

Mass flow (choked): $\dot m = \frac{P_0 A^*}{\sqrt{T_0}}\sqrt{\frac{\gamma}{R}}\left(\frac{2}{\gamma+1}\right)^{\frac{\gamma+1}{2(\gamma-1)}}$ $= \frac{10^6 \times 10^{-3}}{\sqrt{600}}\sqrt{\frac{1.4}{287}}(0.8333)^{3}$ $= \frac{1000}{24.49}(0.06986)(0.5787) = 40.83 \times 0.04043 = 1.651\ \text{kg/s}$ $\dot m \approx 1.65\ \text{kg/s}$. (2 marks)

(c) Normal shock at $A_s = 20\ \text{cm}^2$ (14 marks)

Before the shock the flow is supersonic with $A_s/A^* = 20/10 = 2.0$. Supersonic root: $M_1 \approx 2.20$ (4 marks)

Rankine–Hugoniot downstream Mach: $M_2^2 = \frac{1 + \frac{\gamma-1}{2}M_1^2}{\gamma M_1^2 - \frac{\gamma-1}{2}} = \frac{1 + 0.2(4.84)}{1.4(4.84) - 0.2} = \frac{1.968}{6.576} = 0.2993$ $M_2 = 0.547$ (3 marks)

Static pressure ratio: $\frac{P_2}{P_1} = 1 + \frac{2\gamma}{\gamma+1}(M_1^2 - 1) = 1 + \frac{2.8}{2.4}(3.84) = 1 + 4.48 = 5.48$ (3 marks)

Stagnation-pressure ratio (total-pressure loss): $\frac{P_{0,2}}{P_{0,1}} = \left[\frac{\frac{\gamma+1}{2}M_1^2}{1+\frac{\gamma-1}{2}M_1^2}\right]^{\frac{\gamma}{\gamma-1}}\left[\frac{\gamma+1}{2\gamma M_1^2-(\gamma-1)}\right]^{\frac{1}{\gamma-1}}$ With $M_1 = 2.20$: first bracket $= \left(\dfrac{1.2\times4.84}{1.968}\right)^{3.5} = (2.951)^{3.5}=44.8$; second bracket $=\left(\dfrac{2.4}{6.576}\right)^{2.5}=(0.3650)^{2.5}=0.0805$; product $\approx 3.607$...

Correcting: $P_{0,2}/P_{0,1} = 44.8 \times 0.0805/?$ — proper formula gives $P_{0,2}/P_{0,1} \approx 0.628$ at $M_1=2.20$ (standard normal-shock table value). (4 marks)

Accept table value $P_{0,2}/P_{0,1} \approx 0.628$.

(d) Post-shock exit & expansion state (8 marks)

Across the shock, stagnation pressure drops but stagnation temperature is conserved; the entropy rise means the effective sonic area increases: $A^*_2 = A^* \cdot (P_{0,1}/P_{0,2}) = 10 \times (1/0.628) = 15.9\ \text{cm}^2$. (3 marks)

Downstream, subsonic isentropic flow uses $A_e/A^*_2 = 25/15.9 = 1.57$. The subsonic root gives $M_e \approx 0.40$, and $P_e = P_{0,2}/(1+0.2 M_e^2)^{3.5}$. With $P_{0,2}=0.628\ \text{MPa}$: $P_e \approx 0.628/1.117 = 0.562\ \text{MPa}$. (3 marks)

Compared with the design exit pressure of 64.4 kPa, the actual exit pressure is far higher → the flow has been compressed by the internal shock; relative to design (supersonic) operation the nozzle is over-expanded (back pressure too high for full supersonic expansion, forcing an internal shock). (2 marks)

Question 2

(a) Oblique shock, $\theta = 10^\circ$, $M_1=3$ (12 marks)

$\theta$–$\beta$–$M$ relation: $\tan\theta = 2\cot\beta\,\frac{M_1^2\sin^2\beta - 1}{M_1^2(\gamma+\cos 2\beta)+2}$ Solving the weak root gives $\beta \approx 27.4^\circ$. (4 marks)

Normal component upstream: $M_{n1} = M_1\sin\beta = 3\sin27.4^\circ = 3(0.4602)=1.381$. (2 marks)

Downstream normal Mach: $M_{n2}^2 = \frac{1+0.2 M_{n1}^2}{1.4 M_{n1}^2 - 0.2} = \frac{1+0.2(1.907)}{1.4(1.907)-0.2}=\frac{1.381}{2.470}=0.5593,\ M_{n2}=0.748$ $M_2 = \frac{M_{n2}}{\sin(\beta-\theta)} = \frac{0.748}{\sin17.4^\circ}=\frac{0.748}{0.2990}=2.50$ (3 marks)

Pressure: $\frac{P_2}{P_1}=1+\frac{2\gamma}{\gamma+1}(M_{n1}^2-1)=1+\frac{2.8}{2.4}(0.907)=1+1.058=2.058$ $P_2 = 2.058 \times 20 = 41.2\ \text{kPa}$. (3 marks)

(b) Prandtl–Meyer expansion of $20^\circ$ (10 marks)

$\nu(M)=\sqrt{\frac{\gamma+1}{\gamma-1}}\tan^{-1}\sqrt{\frac{\gamma-1}{\gamma+1}(M^2-1)}-\tan^{-1}\sqrt{M^2-1}$

$\nu(M_2=2.50)=39.12^\circ$. (3 marks)

After turning $20^\circ$: $\nu(M_3)=39.12+20=59.12^\circ$. Inverting gives $M_3 \approx 3.51$. (3 marks)

Isentropic (stagnation pressure constant through expansion; $P_{0,2}$ from behind the shock): $\frac{P_3}{P_2}=\left(\frac{1+0.2M_2^2}{1+0.2M_3^2}\right)^{3.5}=\left(\frac{1+0.2(6.25)}{1+0.2(12.32)}\right)^{3.5}=\left(\frac{2.25}{3.464}\right)^{3.5}=(0.6496)^{3.5}=0.2213$ $P_3 = 0.2213 \times 41.2 = 9.12\ \text{kPa}$. (4 marks)

(c) Python solver (8 marks)

import numpy as np
from scipy.optimize import brentq
 
def prandtl_meyer(M, gamma=1.4):
    if M < 1: raise ValueError("M must be >= 1")
    g = gamma
    term = np.sqrt((g+1)/(g-1))
    nu = (term*np.arctan(np.sqrt((g-1)/(g+1)*(M**2-1)))
          - np.arctan(np.sqrt(M**2-1)))
    return np.degrees(nu)
 
def theta_beta_M(beta_rad, M1, theta_rad, gamma=1.4):
    b, g = beta_rad, gamma
    num = M1**2*np.sin(b)**2 - 1
    den = M1**2*(g+np.cos(2*b)) + 2
    return 2/np.tan(b)*num/den - np.tan(theta_rad)
 
def solve_oblique(M1, theta_deg, gamma=1.4):
    th = np.radians(theta_deg)
    mu = np.arcsin(1/M1)              # Mach angle: lower bound
    # weak solution lies between Mach angle and the theta_max angle (<~ 65 deg).
    # Bracket just above mu up to ~65 deg and take the FIRST root -> weak.
    beta = brentq(theta_beta_M, mu+1e-4, np.radians(65),
                  args=(M1, th, gamma))
    return np.degrees(beta)
 
# solve_oblique(3.0, 10) -> ~27.4 ; prandtl_meyer(2.5) -> ~39.1

*(6 mar