Isentropic, variable-area flows

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:

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

Combined together, we have

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

we obtain a relationship between area change and density change:

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

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

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 (91), we see that $dp < 0$ results in either a positive or negative area change for the subsonic and supersonic regimes:

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

and

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

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:

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:

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

Putting the expressions for $\frac{p_1}{p_2}$ and $\frac{T_2}{T_1}$ back into Equation (96), 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:

Combining Equation (101) with (100) and (98), along with the stagnation relationships

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 (100) and (98) with the ideal gas law to get a density ratio:

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 (98) through (103).

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 (101) 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:

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:

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

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

Isentropic relations¶

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

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 (107) 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$):

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:

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

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:

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$.

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:

1. Find $M_2$ associated with the area ratio $A_2 / A_2^*$, then

2. Use that to find the stagnation pressure ratio $p_2 / p_{t2}$.

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

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

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$:

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

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.

Option 3: solving the full equation¶

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

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 (101) 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