Isentropic, variable-area flows¶

# Necessary modules to solve problems
import numpy as np
import pandas as pd
from scipy.optimize import root_scalar

%matplotlib inline
from matplotlib import pyplot as plt

# 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'

Varying-area adiabatic flows¶

Area change is one of the important factors that can adjust flow properties in compressible flow systems. (The others are friction and heat transfer, which we will cover briefly later.)

For now, let’s proceed to analyze flows with the following conditions/assumptions:

steady, one-dimensional flow
adiabatic: $ \delta q = 0 $, $ d s_e = 0$
no shaft work: $ \delta w_s = 0 $
no or negligible change in potential energy: $ dz = 0 $
no losses (i.e., reversible): $d s_i = 0$

As a result of the flow being adiabatic and reversible, it is also isentropic: $ds = 0$.

Our goal is now to see how changes in pressure, density, and velocity relate with changing area.

Starting with the energy equation, apply our conditions:

\[\begin{gather*} \delta q = \delta w_s + dh + \frac{dV^2}{2} + g dz \\ dh = -V dV \end{gather*}\]

We also have our thermodynamic relationships derived from Gibbs’ identities:

\[\begin{gather*} T ds = dh - \frac{dp}{\rho} \\ dh = \frac{dp}{\rho} \end{gather*}\]

Combined together, we have

\[ dV = - \frac{dp}{\rho V} \;. \]

Substituting this, and the speed of sound ($ dp = a^2 d\rho $) into the continuity equation:

\[\begin{align*} 0 &= \frac{d\rho}{\rho} + \frac{dA}{A} + \frac{dA}{A} \\ \frac{dp}{\rho} &= V^2 \left( \frac{d\rho}{\rho} + \frac{dA}{A} \right) \\ \frac{d\rho}{\rho} &= M^2 \left( \frac{d\rho}{\rho} + \frac{dA}{A} \right) \end{align*}\]

we obtain a relationship between area change and density change:

(9)¶\[ \frac{d\rho}{\rho} = \left( \frac{M^2}{1-M^2} \right) \frac{dA}{A} \;. \]

Substituting this back into the continuity equation, we obtain a relationship between area change and velocity change:

(10)¶\[ \frac{dV}{V} = -\left( \frac{1}{1-M^2} \right) \frac{dA}{A} \;. \]

And, finally, recalling that $ dV = -dp/\rho V $, we substitute that into Equation (10) to get a relationsihp between area change and pressure change:

(11)¶\[ dp = \rho V^2 \left( \frac{1}{1-M^2} \right) \frac{dA}{A} \;. \]

If we focus on situations where the pressure is decreasing ($ dp < 0 $), going from a high-pressure reservoir to a low-pressure receiver, we can examine how changes in the other properties must occur.

Examining Equation (11), we see that $ dp < 0 $ results in either a positive or negative area change for the subsonic and supersonic regimes:

\[ (-) = \frac{1}{1-M^2} \frac{dA}{A} \;, \]

so if $ M < 1 $ then $ dA < 0 $, and if $ M > 1 $ then $ dA > 0 $. Using Equations (9) and (10) we can see how changes in density and velocity occur in the different regimes as well:

\[\begin{gather*} M < 1: \quad \frac{d\rho}{\rho} = (+)(-) \rightarrow d\rho < 0 \\ M > 1: \quad \frac{d\rho}{\rho} = (-)(+) \rightarrow d\rho < 0 \end{gather*}\]

and

\[\begin{gather*} M < 1: \quad \frac{dV}{V} = -(+)(-) \rightarrow dV > 0 \\ M > 1: \quad \frac{dV}{V} = -(-)(+) \rightarrow dV > 0 \end{gather*}\]

In summary, for decreasing pressure ($ dp < 0 $), we have:

property	$ M < 1 $	$ M > 1 $
$ A $	↓	↑
$ \rho $	↓	↓
$ V $	↑	↑

This particular example shows how properties change in a nozzle, which converts pressure/enthalpy to kinetic energy. We can see that a subsonic nozzle has a converging shape (i.e., has a decreasing area) while a supersonic nozzle is diverging (with an increasing area).

In contrast, a diffuser converts kinetic energy into enthalpy/pressure, and is associated with increasing pressure ($ dp > 0 $). A subsonic diffuser is diverging, while a supersonic diffuser is converging.

For propulsion applications, when we need to accelerate a gas from low speed to supersonic speeds, we will need a converging-diverging nozzle. This will be the focus of deeper analysis later.

Equations for perfect gases¶

Using our governing equations, we can generate working equations that apply to more-general flows of perfect/ideal gases, where losses (i.e., irreversibilities) may be present.

The flow assumptions are:

steady, one-dimensional flow
adiabatic
no shaft work
perfect/ideal gas
no, or negligible, potential energy changes

Our goal is to find relations between properties at two points in the flow that are only a function of the gas, Mach numbers at both locations, and entropy change between the two locations: $ f \left( M_1, M_2, \gamma, \delta s \right) $.

Start with the continuity equation applied to a control volume, with flow entering at location 1 and leaving at location 2:

\[ \rho_1 A_1 V_1 = \rho_2 A_2 V_2 \]

(12)¶\[ \rightarrow \frac{A_2}{A_1} = \frac{\rho_1 V_1}{\rho_2 V_2} = \frac{p_1}{p_2} \frac{M_1}{M_2} \left( \frac{T_2}{T_1} \right)^{1/2} \;, \]

where we got the final relation by using the ideal gas law ($ p = \rho R T $), Mach number definition ($ V = M a $), and speed of sound expression in an ideal gas ($ a^2 = \gamma R T $).

We can simplify this expression even further by finding ways to express the pressure and temperature ratios as functions of the Mach numbers and gas properties.

For the temperature ratio, we can use conservation of energy and our expression for stagnation temperature:

\[\begin{gather*} h_{t1} + q = h_{t2} + w_s \\ \rightarrow T_{t1} = T_{t2} \\ T_t = T \left(1 + \frac{\gamma - 1}{2} M^2 \right) \\ \end{gather*}\]

(13)¶\[ \therefore \frac{T_2}{T_1} = \frac{1 + \frac{\gamma-1}{2} M_1^2}{1 + \frac{\gamma-1}{2} M_2^2} \]

For pressure, recall our relationship for the ratio of stagnation pressures between two locations:

\[\begin{gather*} \frac{p_{t2}}{p_{t1}} = e^{-\Delta s / R} \\ p_t = p \left(1 + \frac{\gamma - 1}{2} M^2 \right)^{\gamma / (\gamma-1)} \\ \frac{p_{t2}}{p_{t1}} = \frac{p_2}{p_1} \left[ \frac{1 + \frac{\gamma - 1}{2} M_2^2}{1 + \frac{\gamma - 1}{2} M_1^2} \right]^{\gamma / (\gamma-1)} = e^{-\Delta s / R} \end{gather*}\]

(14)¶\[ \therefore \frac{p_1}{p_2} = \left[ \frac{1 + \frac{\gamma - 1}{2} M_2^2}{1 + \frac{\gamma - 1}{2} M_1^2} \right]^{\gamma / (\gamma-1)} e^{\Delta s / R} \;. \]

Putting the expressions for $ \frac{p_1}{p_2} $ and $ \frac{T_2}{T_1} $ back into Equation (12), we can obtain a final expression for the area ratio between two locations as a function of the Mach numbers at the locations, the specific heat ratio of the gas, and any entropy increase between the locations:

(15)¶\[ \frac{A_2}{A_1} = \frac{M_1}{M_2} \left[ \frac{1 + \frac{\gamma-1}{2} M_2^2}{1 + \frac{\gamma-1}{2} M_1^2} \right]^{\frac{\gamma+1}{2(\gamma-1)}} e^{\Delta s/R} \;. \]

Combining Equation (15) with (14) and (13), along with the stagnation relationships

\[\begin{gather*} T_{t2} = T_{t1} \\ \frac{p_{t2}}{p_{t1}} = e^{-\Delta s / R} \;, \end{gather*}\]

we have relationships between properties at two locations in any steady, one-dimensional flow with varying area, as long as there is no heat transfer or shaft work. We could also combine Equations (14) and (13) with the ideal gas law to get a density ratio:

(16)¶\[ \frac{\rho_2}{\rho_1} = \left[ \frac{1 + \frac{\gamma - 1}{2} M_1^2}{1 + \frac{\gamma - 1}{2} M_2^2} \right]^{1 /(\gamma-1)} e^{-\Delta s / R} \;. \]

Sonic reference state¶

If we know the Mach numbers at two locations, it is pretty easy to find the associated ratios of area, temperature, pressure, and/or density using Equations (13) through (16).

However, if we instead know the area ratio or the properties at the two locations, plus one Mach number, and want to find the unknown Mach number, these equations are a bit tougher to solve. We can solve them directly using numerical methods, as you will see in a bit, but we can also take advantage of a helpful concept to ease the calculations.

Similar to the stagnation reference state, or the state associated with decelerating a fluid flow to rest and zero potential, isentropically, we can introduce the sonic reference state, indicated using a * superscript. This is the state associated with either decelerating or accelerating a fluid to sonic velocity (i.e., Mach = 1.0) by some process; currently, that will be through isentropic area change.

This reference state may not physically exist in the system, but could through appropriate area change. As a result, they are legitimate locations in a flow system, and so we can use all of our developed equations to consider flow from a real location to this reference location.

In particular, we can write Equation (15) between the reference locations associated with two real locations: $ A_1 $ and $ A_1^* $, then $ A_2 $ and $ A_2^* $. By definition, the Mach numbers are 1 at the two reference states. So, the area ratio equation becomes:

(17)¶\[ \frac{A_2^*}{A_1^*} = e^{\Delta s / R} \;. \]

This equation expresses how the area associated with reference state changes in a flow system in the presence of losses; for isentropic flow, the reference area $ A^* $ remains constant.

Recall our relationship for the stagnation pressure ratio between two locations:

\[ \frac{p_{t2}}{p_{t1}} = e^{-\Delta s / R} \;. \]

If we take the product of these two equations, we get

(18)¶\[ p_{t1} A_1^* = p_{t2} A_2^* \;, \]

which gives us a quantity that is conserved in any adiabatic flow.

Isentropic relations¶

For isentropic flows, where $\Delta s = 0$, Equation (15) reduces to

(19)¶\[ \frac{A_2}{A_1} = \frac{M_1}{M_2} \left( \frac{1+\frac{\gamma-1}{2} M_2^2}{1+\frac{\gamma-1}{2} M_1^2} \right)^{\frac{\gamma+1}{2(\gamma-1)}} \;. \]

If we have $M_1$ and $M_2$, and know the ideal gas (and therefore $\gamma$) then finding the area ratio between two sections is straightforward using the equation. Working backwards takes slightly more effort, since the equation cannot simply be inverted.

Instead, it either needs to be solved numerically using a root-finding algorithm, or we can take advantage of the * (sonic) reference state to provide an easy way to solve problems.

Let’s apply Equation (19) between any point in the flow and its sonic reference state, such that $A_2 \rightarrow A$ (so $M_2 \rightarrow M$) and $A_1 \rightarrow A^*$ (and $M_1 \rightarrow 1$):

(20)¶\[ \frac{A}{A^*} = \frac{1}{M} \left( \frac{1+\frac{\gamma-1}{2} M^2}{\frac{\gamma+1}{2}} \right)^{\frac{\gamma+1}{2(\gamma-1)}} \;, \]

which shows that $\frac{A}{A^*} = f(\gamma, M)$. For a given value of $\gamma$ (such as 1.4, used for air), it is easy to precalculate and tabulate the reference area ratio vs. values of Mach number. Similarly, the quantities $p/p_t$, $T/T_t$, and $pA/p_t A^*$ can be tabulated as well:

(21)¶\[ \frac{p}{p_t} = \left(1 + \frac{\gamma-1}{2} M^2 \right)^{-\gamma/(\gamma-1)} \]

(22)¶\[ \frac{T}{T_t} = \left(1 + \frac{\gamma-1}{2} M^2 \right)^{-1} \]

For example, let’s tabulate the values from $M = $ 0 to 0.3:

gamma = 1.4
machs = np.arange(0, 0.31, 0.01)
# hiding divide by zero warning
with np.errstate(divide='ignore'):
    area_ratios = (
        (1.0/machs)*(2*(1.0 + 0.5*(gamma-1)*machs**2)/(gamma+1))**(0.5*(gamma+1)/(gamma-1))
        )
pressure_ratios = (1.0/(1 + 0.5*(gamma-1)*machs**2))**(gamma/(gamma-1))
temperature_ratios = 1.0 / (1 + 0.5*(gamma-1)*machs**2)
df = pd.DataFrame({
    r'$M$': machs, r'$p/p_t$': pressure_ratios, r'$T/T_t$': temperature_ratios,
    r'$A/A^*$': area_ratios
    })
df.style.\
    hide_index().\
    format(formatter={(r'$M$'): "{:.2f}"})

$M$	$p/p_t$	$T/T_t$	$A/A^*$
0.00	1.000000	1.000000	inf
0.01	0.999930	0.999980	57.873843
0.02	0.999720	0.999920	28.942130
0.03	0.999370	0.999820	19.300542
0.04	0.998881	0.999680	14.481486
0.05	0.998252	0.999500	11.591444
0.06	0.997484	0.999281	9.665910
0.07	0.996578	0.999021	8.291525
0.08	0.995533	0.998722	7.261610
0.09	0.994351	0.998383	6.461342
0.10	0.993031	0.998004	5.821829
0.11	0.991576	0.997586	5.299230
0.12	0.989985	0.997128	4.864318
0.13	0.988259	0.996631	4.496859
0.14	0.986400	0.996095	4.182400
0.15	0.984408	0.995520	3.910343
0.16	0.982285	0.994906	3.672739
0.17	0.980030	0.994253	3.463509
0.18	0.977647	0.993562	3.277926
0.19	0.975135	0.992832	3.112259
0.20	0.972497	0.992063	2.963520
0.21	0.969733	0.991257	2.829294
0.22	0.966845	0.990413	2.707602
0.23	0.963835	0.989531	2.596812
0.24	0.960703	0.988611	2.495562
0.25	0.957453	0.987654	2.402710
0.26	0.954085	0.986660	2.317287
0.27	0.950600	0.985630	2.238471
0.28	0.947002	0.984562	2.165554
0.29	0.943291	0.983458	2.097928
0.30	0.939470	0.982318	2.035065

With tabulated data like this, we can then solve problems by constructing an appropriate property ratio based on known information and looking up the corresponding Mach number.

For example, given a gas ($\gamma$), the areas at two locations ($A_1$ and $A_2$), and the Mach number at one location ($M_1$), we can find the Mach number at the second location ($M_2$) by constructing the reference area ratio:

\[ \frac{A_2}{A_2^*} = \frac{A_2}{A_1} \frac{A_1}{A_1^*} \frac{A_1^*}{A_2^*} \;, \]

where $\frac{A_2}{A_2^*} = f(M_2)$, $\frac{A_2}{A_1}$ is known, $\frac{A_1}{A_1^*} = f(M_1)$, and $\frac{A_1^*}{A_2^*} = 1.0$ for isentropic flow between the two sections. So, we can find $M_2$.

Note

For a given area ratio $\frac{A}{A^*}$, there will always be two possible Mach numbers: a subsonic solution and a supersonic solution.

To determine which Mach number is the correct solution, you need to apply other knowledge and problem context. For example, if the initial Mach number is subsonic and the two sections are only connected by a diverging or converging area duct, then the second Mach number must also be subsonic.

Example: isentropic flow problem¶

Problem: Air flows isentropically through a duct ($\gamma = 1.4$) where the area is changing from point 1 to 2, with no heat transfer or shaft work. The area ratio is $\frac{A_2}{A_1} = 2.5$, the flow starts at $M_1 = 0.5$ and 4 bar. Find the Mach number and pressure at the second point in the duct.

We can solve this using the classical approach (pre-calculated isentropic tables) or a numerical approach; both follow the same general approach:

Find $M_2$ associated with the area ratio $A_2 / A_2^*$, then
Use that to find the stagnation pressure ratio $p_2 / p_{t2}$.

\[ \frac{A_2}{A_2^*} = \frac{A_2}{A_1} \frac{A_1}{A_1^*} \frac{A_1^*}{A_2^*} \;, \]

where $\frac{A_2}{A_1} = 2.5$ is given, we can find $\frac{A_1}{A_1^*}$ using

\[ \frac{A}{A^*} = \frac{1}{M} \left( \frac{1 + \frac{\gamma - 1}{2} M^2}{\frac{\gamma+1}{2}} \right)^{\frac{\gamma+1}{2(\gamma-1)}} \;, \]

(either by calculating or looking up in the $\gamma = 1.4$ table) and $\frac{A_1^*}{A_2^*} = 1$ because the flow is isentropic.

gamma = 1.4
mach_1 = 0.5

A2_A1 = 2.5
A1star_A2star = 1.0 # isentropic

A1_A1star = (1.0/mach_1) * (
    (1 + 0.5*(gamma-1)*mach_1**2) / ((gamma + 1)/2)
    )**((gamma+1) / (2*(gamma-1)))
print(f'A1/A1^* = {A1_A1star:.4f}')

A1/A1^* = 1.3398

A2_A2star = A2_A1 * A1_A1star * A1star_A2star
print(f'A2/A2star = {A2_A2star:.4f}')

A2/A2star = 3.3496

We can then find $M_2$, because $\frac{A_2}{A_2*} = f(M_2)$. Our options are to use the $\gamma = 1.4$ tables and interpolate, or solve the associated equation numerically.

Option 1: using tables¶

We can find in the tables that:

at $M=0.17$, $A/A^* = 3.46351$
at $M = 0.18$, $A/A^* = 3.27793$

and interpolate to find the precise $M_2$:

machs = np.array([0.17, 0.18])
areas = np.array([3.46351, 3.27793])
mach_2 = (
    machs[0] * (areas[1] - A2_A2star) + machs[1] * (A2_A2star - areas[0])
    ) / (areas[1] - areas[0])
print(f'M2 = {mach_2:.4f}')

M2 = 0.1761

This is probably sufficient, but we could get a more-accurate result by interpolating using more points and using the numpy.interp() function:

machs = np.array([0.15, 0.16, 0.17, 0.18, 0.19])
areas = np.array([3.91034, 3.67274, 3.46351, 3.27793, 3.11226])

mach_2 = np.interp(A2_A2star, areas[::-1], machs[::-1])
print(f'M2 = {mach_2:.4f}')

M2 = 0.1761

Note that we have to reverse the order of the values, since interp expects the x-values to be increasing. Also, we could easily generate these values ourselves for a different value of $\gamma$, but it is likely easier to just solve the equation directly in that case

Option 2: solving the equation¶

Alternately, we can solve the $\frac{A_2}{A_2*} = f(M_2, \gamma)$ equation directly using scipy.optimize.root_scalar. We need to give the function initial guesses x0 and x1 that are subsonic, to ensure we get a subsonic solution.

def area_function(mach, gamma, area_ratio):
    '''Function for area ratio, solving for M2'''
    return (
        area_ratio - (
            (1.0/mach) * ((1 + 0.5*(gamma-1)*mach**2) / 
            ((gamma + 1)/2))**((gamma+1) / (2*(gamma-1)))
            )
        )

sol = root_scalar(area_function, args=(gamma, A2_A2star), x0=0.1, x1=0.5)
print(f'M2 = {sol.root:.4f}')

M2 = 0.1760

Warning

Make sure you check the output from root-finding functions, because they can converge on the wrong solution even with an appropriate initial guess. You might need to adjust the initial guess if you get a solution that is not consistent with the physics of the problem.

Option 3: solving the full equation¶

As a third option, we could actually bypass the ratio approach and just numerically solve the full Equation (15) directly, using the known $ A_2/A_1 $, $ M_1 $, $ \gamma $, and $ \Delta s = 0 $!

def full_area_function(mach2, mach1, area_ratio, gamma=1.4, delta_s=0.0, R=287.0):
    '''Function for full area ratio equation, solving for M2'''
    # gamma, delta_s, and R have default values
    return (
        area_ratio - (mach1 / mach2) * (
            (1 + 0.5*(gamma-1)*mach2**2) / (1 + 0.5*(gamma-1)*mach1**2)
            )**((gamma+1) / (2*(gamma-1))) * np.exp(-delta_s / R)
        )

sol = root_scalar(full_area_function, args=(mach_1, A2_A1, gamma, 0.0), x0=0.1, x1=0.5)
print(f'M2 = {sol.root:.4f}')

M2 = 0.1760

This approach is the most direct.

Example: isentropic flow (using Matlab)¶

If you prefer to use Matlab, we can also solve the same example using the three different strategies.

gamma = 1.4;
mach_1 = 0.5;

A2_A1 = 2.5;
A1star_A2star = 1.0;

A1_A1star = (1.0/mach_1)*((1+0.5*(gamma-1)*mach_1^2)/((gamma+1)/2))^((gamma+1)/(2*(gamma-1)));
fprintf('A1/A1^* = %.4f\n', A1_A1star)

A2_A2star = A2_A1 * A1_A1star * A1star_A2star;
fprintf('A2/A2star = %.4f\n', A2_A2star)

A1/A1^* = 1.3398
A2/A2star = 3.3496

Option 1: using the tables¶

We find in the tables that:

at $M=0.17$, $A/A^* = 3.46351$
at $M = 0.18$, $A/A^* = 3.27793$

and linearly interpolate to find the precise $M_2$:

machs = [0.17, 0.18];
areas = [3.46351, 3.27793];
mach_2 = (machs(1)*(areas(2)-A2_A2star) + machs(2)*(A2_A2star-areas(1)))/(areas(2)-areas(1));
fprintf('M2 = %.4f\n', mach_2)

M2 = 0.1761

We can also use Matlab’s built-in interp1 function to interpolate with more data points:

machs = [0.15, 0.16, 0.17, 0.18, 0.19];
areas = [3.91034, 3.67274, 3.46351, 3.27793, 3.11226];

mach_2 = interp1(areas, machs, A2_A2star);
fprintf('M2 = %.4f\n', mach_2)

M2 = 0.1761

Option 2: solving the equation¶

Alternately, we can solve the $\frac{A_2}{A_2*} = f(M_2, \gamma)$ equation directly using the function fzero:

f = @(M) (A2_A2star - ((1.0/M)*((1+0.5*(gamma-1)*M^2) ...
    / ((gamma+1)/2))^((gamma+1)/(2*(gamma-1)))));

M2 = fzero(f, 0.1);
fprintf('M2 = %.4f\n', M2)

M2 = 0.1760

Note that, in this example, the values of $\gamma$ and $A_2 / A_2^*$ are hard-coded into the function.

Option 3: solving the full equation¶

As a third option, we can bypass the ratio approach and just numerically solve the full Equation (15) directly, using the known $ A_2/A_1 $, $ M_1 $, $ \gamma $, and $ \Delta s = 0 $!

To do this, we create two functions: one that depends on all the variables, and a parameterized function that only depends on $M_2$.

fun = @(M2, M1, A2_A1, gamma, delta_s, R) (A2_A1 - (M1/M2)*( ...
      (1 + 0.5*(gamma-1)*M2^2) / (1 + 0.5*(gamma-1)*M1^2) ...
            )^((gamma+1)/(2*(gamma-1)))*exp(-delta_s/R));

% parameterized function
f = @(M2) fun(M2, mach_1, A2_A1, gamma, 0.0, 287);
M2 = fzero(f, 0.1);
fprintf('M2 = %.4f\n', M2)

M2 = 0.1760

Gas Dynamics notes

Isentropic, variable-area flows¶

Varying-area adiabatic flows¶

Equations for perfect gases¶

Sonic reference state¶

Isentropic relations¶

Example: isentropic flow problem¶

Option 1: using tables¶

Option 2: solving the equation¶

Option 3: solving the full equation¶

Example: isentropic flow (using Matlab)¶

Option 1: using the tables¶

Option 2: solving the equation¶

Option 3: solving the full equation¶