Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

Prandtl–Meyer flows

Notebook Cell
%matplotlib inline
from matplotlib import pyplot as plt

import numpy as np
from scipy.optimize import root_scalar

# Pint gives us some helpful unit conversion
from pint import UnitRegistry
ureg = UnitRegistry()
Q_ = ureg.Quantity # We will use this to construct quantities (value + unit)
Notebook Cell
# these lines are only for helping improve the display
import matplotlib_inline.backend_inline
matplotlib_inline.backend_inline.set_matplotlib_formats('pdf', 'png')
plt.rcParams['figure.dpi']= 300
plt.rcParams['savefig.dpi'] = 300
plt.rcParams['mathtext.fontset'] = 'cm'

Prandtl–Meyer flows are isentropic turning flows that can involve compression or expansion. To describe them, we examine shocks at their weak limit.

Analysis of weak shocks

We can express the pressure change through a normal shock using Equation (143):

p2p1=2γγ+1M12γ1γ+1p2p11=Δpp1=2γγ+1(M121)  ,\begin{gather*} \frac{p_2}{p_1} = \frac{2\gamma}{\gamma+1} M_1^2 - \frac{\gamma-1}{\gamma+1} \\ \frac{p_2}{p_1} - 1 = \frac{\Delta p}{p_1} = \frac{2\gamma}{\gamma+1} (M_1^2 - 1) \;, \end{gather*}

so therefore through a normal shock Δp(M121) \Delta p \propto (M_1^2 - 1) .

We can also examine the entropy change through a normal shock by using one of our thermodynamic state relationships:

Δs12=cplogv2v1+cvlogp2p1s2s1R=γγ1logρ1ρ2+1γ1logp2p1  ,\begin{gather*} \Delta s_{1-2} = c_p \log \frac{v_2}{v_1} + c_v \log \frac{p_2}{p_1} \\ \frac{s_2 - s_1}{R} = \frac{\gamma}{\gamma-1} \log \frac{\rho_1}{\rho_2} + \frac{1}{\gamma-1} \log \frac{p_2}{p_1} \;, \end{gather*}

which, after considerable effort, we can express as

s2s1R=2γ(M121)33(γ+1)2+O(M121)4  .\frac{s_2 - s_1}{R} = \frac{2\gamma \left(M_1^2 - 1\right)^3}{3 (\gamma+1)^2} + \mathcal{O} \left( M_1^2 - 1 \right)^4 \;.

For weak shocks, where M11M_1 \rightarrow 1, we can neglect the higher-order terms, and so Δs(M121)3 \Delta s \propto (M_1^2 - 1)^3 .

For turning flows, the above analysis applies to the normal component of the Mach number: M1n=M1sinθM_{1n} = M_1 \sin \theta:

Δpp1=2γγ+1(M12sin2θ1)\frac{\Delta p}{p_1} = \frac{2\gamma}{\gamma+1} (M_1^2 \sin^2 \theta - 1)

and

ΔsR=2γ(M12sin2θ1)33(γ+1)2+O(M12sin2θ1)4  .\frac{\Delta s}{R} = \frac{2\gamma \left(M_1^2 \sin^2 \theta - 1\right)^3}{3 (\gamma+1)^2} + \mathcal{O} \left( M_1^2 \sin^2 \theta - 1 \right)^4 \;.

For very weak oblique shocks, we can make specific approximations:

  1. The deflection angle δ\delta is very small (δ1 \delta \ll 1 ), meaning we can use the small angle approximations to say tanδδ \tan \delta \approx \delta .

  2. The shock angle θ\theta approaches the Mach angle μ\mu, where sinμ=1/M\sin \mu = 1 / M.

Let’s apply these approximations to our oblique-shock relationship given by Equation (150):

tanδ=2cotθ[M12sin2θ1M12(γ+cos2θ)+2)δ(constants)(M12sin2θ1)δ(M12sin2θ1)  ,\begin{gather*} \tan \delta = 2 \cot \theta \left[ \frac{M_1^2 \sin^2 \theta - 1}{M_1^2 (\gamma + \cos 2\theta) + 2} \right) \\ \rightarrow \delta \approx (\text{constants}) \left( M_1^2 \sin^2 \theta - 1 \right) \\ \delta \propto \left( M_1^2 \sin^2 \theta - 1 \right) \;, \end{gather*}

since μ\mu is known for a given value of M1M_1. Therefore, we can express the changes in pressure and entropy using the deflection angle:

Δpp12γγ+1(constants)  δΔsR2γ3(γ+1)2[(constants)  δ]3\begin{gather*} \frac{\Delta p}{p_1} \approx \frac{2\gamma}{\gamma+1} (\text{constants}) \; \delta \\ \frac{\Delta s}{R} \approx \frac{2\gamma}{3 (\gamma+1)^2} \left[ (\text{constants}) \; \delta \right]^3 \end{gather*}

Now, we can say that for any weak oblique shock, for any initial conditions, that

Δpδ  .Δsδ3  .\begin{gather*} \Delta p \propto \delta \phantom{\;.} \\ \Delta s \propto \delta^3 \;. \end{gather*}

Isentropic turns

So far we have been analyzing a single weak shock, which deflects the flow by an infinitesimal angle. Now, consider a flow undergoing a turn over a finite angle, but comprised of nn small turns.

Diagram showing many small shock turns

Figure 4:A finite turning flow accomplished with many small turns from weak shocks.

With nn equal segments, the overall turning angle is δT=nδ \delta_T = n \delta. For a given δT\delta_T, increasing the number of segments nn will mean that each deflection δ\delta decreases. Each segment is a very weak oblique shock, and for each Δpδ \Delta p^{\prime} \propto \delta and Δsδ3 \Delta s^{\prime} \propto \delta^3 .

Then, if we look at the total change in pressure and entropy over the whole turn:

Δp=nΔpnδΔs=nΔsnδ3  ,\begin{gather*} \Delta p = \sum^n \Delta p^{\prime} \propto n \, \delta \\ \Delta s = \sum^n \Delta s^{\prime} \propto n \, \delta^3 \;, \end{gather*}

but δ=δT/n \delta = \delta_T / n , so

Δplimnn(δTn)=δTΔslimnn(δTn)3=0  .\begin{gather*} \Delta p \propto \lim_{n \rightarrow \infty} n \left( \frac{\delta_T}{n} \right) = \delta_T \\ \Delta s \propto \lim_{n \rightarrow \infty} n \left( \frac{\delta_T}{n} \right)^3 = 0 \;. \end{gather*}

Therefore, as the number of segments increases and nn \rightarrow \infty,

  1. The wall makes a smooth turn through a finite angle δT\delta_T,

  2. The shock waves become Mach waves,

  3. The Mach number changes continuously,

  4. There is a finite pressure change continuously through the turn, and

  5. There is no change in entropy—the process is isentropic.

So, we can use smooth turns to achieve isentropic compression with a concave turn, and isentropic expansion with a convex turn—smooth or sharp.

Analysis of Prandtl-Meyer flows

Knowing that Prandtl-Meyer expansion/compression is isentropic, we want to know how the Mach number changes for a given turning angle: ΔM=f(δT) \Delta M = f( \delta_T) . To derive this relationship, we can examine a single Mach wave to see how the flow changes instantaneously in an expansion fan.

Diagram of flow change through a single Mach wave

Figure 5:Diagram of how flow changes through a single Mach wave in a Prandtl-Meyer expansion turn.

We start with a supersonic flow at Mach number MM, encountering a expansion Mach wave at angle μ\mu. The flow turns through angle dνd \nu and accelerates to M+dMM + dM. Looking at the components of velocity, the tangential component VtV_t does not change through the expansion wave:

Vt=Vcosμ=(V+dV)cos(μ+dν)Vcosμ=(V+dV)(cosμcosdνsinμsindν)\begin{gather*} V_t = V \cos \mu = (V + dV) \cos (\mu + d \nu) \\ V \cos \mu = (V + dV) \left( \cos \mu \cos d \nu - \sin \mu \sin d \nu \right) \end{gather*}

Since a single Mach wave will only turn the flow by an infinitesimal angle, we can apply the small-angle approximation: dν1 d\nu \ll 1 , so cosdν1\cos d \nu \approx 1 and sindνdν\sin d\nu \approx d\nu. Therefore,

Vcosμ=(V+dV)(cosμdνsinμ)dν=cotμdVV\begin{gather*} V \cos \mu = (V + dV) \left( \cos \mu - d\nu \sin \mu \right) \\ d \nu = \cot \mu \frac{dV}{V} \end{gather*}

But, since we have that sinμ=1M\sin \mu = \frac{1}{M}, we know that cosμ=M21\cos \mu = \sqrt{M^2 - 1}. Then,

dν=M21dVVdν=M211+γ12M2dMM  .\begin{gather*} d \nu = \sqrt{M^2 - 1} \frac{dV}{V} \\ \rightarrow d \nu = \frac{\sqrt{M^2 - 1}}{1 + \frac{\gamma-1}{2} M^2} \frac{dM}{M} \;. \end{gather*}

This gives us dν=f(γ,M)d\nu = f(\gamma, M), which we can integrate to get ν+const=f(γ,M)\nu + \text{const} = f(\gamma, M). If we specify the constraint that ν=0\nu = 0 for M=1 M = 1, then the integration constant becomes zero, and the closed-form expression for ν\nu is

ν=(γ+1γ1)1/2tan1[γ1γ+1(M21)]1/2tan1(M21)1/2  .\nu = \left( \frac{\gamma+1}{\gamma-1} \right)^{1/2} \tan^{-1} \left[ \frac{\gamma-1}{\gamma+1} \left( M^2 - 1 \right) \right]^{1/2} - \tan^{-1} \left( M^2 - 1\right)^{1/2} \;.

This is the Prandtl–Meyer function, the angle measured from the flow direction where M=1M = 1, through which a flow turned isentropically reaches MM.

If we know the downstream Mach number, we can solve Equation (171) easily to get the associated turning angle. Finding the new Mach number from the angle requires interpolation or numerical solution.

def get_prandtl_meyer(mach, gamma=1.4):
    '''Evaluate Prandtl-Meyer function at given Mach number.
    
    Defined as the angle from the flow direction where Mach = 1 through which the 
    flow turned isentropically reaches the specified Mach number.
    '''
    return (
        np.sqrt((gamma + 1) / (gamma - 1)) *
        np.arctan(np.sqrt((gamma - 1)*(mach**2 - 1)/(gamma + 1))) -
        np.arctan(np.sqrt(mach**2 - 1))
        )

def solve_prandtl_meyer(mach, nu, gamma=1.4):
    '''Solve for unknown Mach number, given Prandtl-Meyer function (in radians).'''
    return (nu - get_prandtl_meyer(mach, gamma))

Example: turn from sonic flow

Consider a horizontal airflow moving at M=1M = 1, which is turned 28° via a sharp corner. A centered Prandtl–Meyer expansion fan isentropically expands the flow. What is the final Mach number?

In this case, we can apply the Prandtl–Meyer function directly: ν\nu = 28° from the flow direction where M=1M = 1. We can then solve for M2M_2:

gamma = 1.4
turn_angle = 28 * np.pi/180

# guesses for M2: 2.0 and 3.0
root = root_scalar(solve_prandtl_meyer, x0=2.0, x1=3.0, args=(turn_angle, gamma))
print(f'M2 = {root.root: .3f}')
M2 =  2.059

Example: turn from general supersonic flow

Now, consider a horizontal airflow moving at a Mach number of 2.059, that encounters a sharp turn downward of 12°. What is the final Mach number?

We cannot just solve Equation (171) because the initial Mach number is not 1.0.

However, we can imagine that to obtain the true starting Mach number of 2.059, we hypothetically started with a flow at Mach = 1.0 and expanded it through a turn of 28°. Then, adding an additional 12° will reach the same final Mach number.

In other words, modify the problem such that M0=1M1M2M_0 = 1 \rightarrow M_1 \rightarrow M_2, and ν2=Δν+ν1\nu_2 = \Delta \nu + \nu_1, where ν1\nu_1 is the Prandtl–Meyer function evaluated at the actual starting Mach number M1M_1, Δν\Delta \nu is the turning angle (12° in this case), and ν2\nu_2 is the Prandtl–Meyer function of the final Mach number M2M_2.

nu_1 = 28 * np.pi/180
delta_nu = 12 * np.pi/180

nu_2 = nu_1 + delta_nu
print(f'nu_2 = {nu_2*180/np.pi: .1f}°')

# guesses for M2: 2.0 and 3.0
root = root_scalar(solve_prandtl_meyer, x0=2.0, x1=3.0, args=(nu_2, gamma))
print(f'M_2 = {root.root: .3f}')
nu_2 =  40.0°
M_2 =  2.538

Therefore, for any Prandtl–Meyer flow, we can find the final Mach number from an arbitrary initial Mach number using

ν2=ν1+Δν  ,\nu_2 = \nu_1 + \Delta \nu \;,

where ν1=f(M1,γ)\nu_1 = f(M_1, \gamma), Δν\Delta \nu is the turning angle (positive for convex expanding flow, negative for concave compressing flow), and ν2=f(M2,γ)\nu_2 = f(M_2, \gamma).

Example: isentropic compression

Consider a Prandtl–Meyer compression problem: a horizontal airflow at Mach 2.40 encounters a smooth turn upward at 20°. What is the Mach number near the wall after the turn?

gamma = 1.4
M1 = 2.40

# turning angle is negative because compression
delta_nu = -20 * np.pi/180

nu_1 = get_prandtl_meyer(M1, gamma)
nu_2 = nu_1 + delta_nu

root = root_scalar(solve_prandtl_meyer, x0=1.5, x1=2.0, args=(nu_2, gamma))
print(f'M_2 = {root.root: .3f}')
M_2 =  1.664

However, what happens away from the wall? The flow needs to turn sharply, and an oblique shock forms with a deflection angle δ\delta of 20°.

def oblique_shock_delta(theta, mach, gamma=1.4):
    '''Calculate oblique shock deflection from shock angle and Mach'''
    return (
        np.arctan((2.0/np.tan(theta)) * (
            (mach**2 * np.sin(theta)**2 - 1) /
            (2 + mach**2 * (gamma + np.cos(2*theta)))
            ))
        )

def solve_oblique_theta(theta, delta, mach, gamma=1.4):
    '''Solve for oblique wave angle theta'''
    return (delta - oblique_shock_delta(theta, mach, gamma))

def get_mach_normal(gamma, mach1):
    '''Calculate Mach number after a normal shock'''
    return np.sqrt(
        (mach1**2 + 2/(gamma - 1))/(2 * gamma * mach1**2/(gamma - 1) - 1)
        )
gamma = 1.4
M1 = 2.40
delta = 20 * np.pi/180

root = root_scalar(
    solve_oblique_theta, x0=40*np.pi/180, x1=50*np.pi/180,
    args=(delta, M1, gamma)
    )
theta = root.root
print(f'Theta: {theta * 180/np.pi: .2f}°')

M1n = M1 * np.sin(theta)
M2n = get_mach_normal(gamma, M1n)
M2 = M2n / np.sin(theta - delta)
print(f'M2: {M2: .3f}')
Theta:  44.34°
M2:  1.569

Due to the oblique shock, the Mach number away from the wall is lower than it is near the wall, and the pressure will be different too! So, there will need to be another wave phenomenon occuring between the two regions.