Flowpedia

• Incompressible flow
• Compressible flow

Incompressible flow subpages ordered by category

Thermodynamics

Specific heat

For thermally perfect and calorically perfect gases

$\begin{array}{cc}& C_p=\frac{dh}{dT}\\ & C_v=\frac{de}{dT}\end{array}$

(Eq. 1)

From the definition of enthalpy and the equation of state $p=\rho RT$

$h=e+\frac{p}{\rho }=e+RT$

(Eq. 2)

Differentiate (Eq. 2) with respect to temperature gives

$\frac{dh}{dT}=\frac{de}{dT}+\frac{d(RT)}{dT}$

(Eq. 3)

Inserting the specific heats gives

$C_p=C_v+R$

(Eq. 4)

Dividing (Eq. 4) by $C_v$ gives

$\frac{C_p}{C_v}=1+\frac{R}{C_v}$

(Eq. 5)

Introducing the ratio of specific heats defined as

$\gamma =\frac{C_p}{C_v}$

(Eq. 6)

Now, inserting (Eq. 6) in Eqn. \ref{eq:specificheat:c} gives

$C_v=\frac{R}{\gamma -1}$

(Eq. 7)

In the same way, dividing (Eq. 4) with $C_p$ gives

$1=\frac{C_v}{C_p}+\frac{R}{C_p}=\frac{1}{\gamma }+\frac{R}{C_p}$

(Eq. 8)

and thus

$C_p=\frac{\gamma R}{\gamma -1}$

(Eq. 9)

Isentropic relations

First law of thermodynamics

First law of thermodynamics:

$de=\delta q-\delta w$

(Eq. 10)

For a reversible process: $\delta w=pd(1/\rho )$ and $\delta q=Tds$

$de=Tds-pd\left(\frac{1}{\rho }\right)$

(Eq. 11)

Enthalpy is defined as: $h=e+p/\rho$ and thus

$dh=de+pd\left(\frac{1}{\rho }\right)+\left(\frac{1}{\rho }\right)dp$

(Eq. 12)

Eliminate $de$ in (Eq. 11) using (Eq. 12)

$Tds=dh-\cancel{pd\left(\frac{1}{\rho }\right)}-\left(\frac{1}{\rho }\right)dp+\cancel{pd\left(\frac{1}{\rho }\right)}$

(Eq. 13)

$ds=\frac{dh}{T}-\frac{dp}{\rho T}$

(Eq. 14)

Using $dh=C_{p}T$ and the equation of state $p=\rho RT$, we get

$ds=C_p\frac{dT}{T}-R\frac{dp}{p}$

(Eq. 15)

Integrating (Eq. 15) gives

$s_2-s_1=\underset{1}{\overset{2}{\int }}{C}_p\frac{dT}{T}-R\ln \left(\frac{p_2}{p_1}\right)$

(Eq. 16)

For a calorically perfect gas, $C_p$ is constant (not a function of temperature) and can be moved out from the integral and thus

$s_2-s_1=C_p\ln \left(\frac{T_2}{T_1}\right)-R\ln \left(\frac{p_2}{p_1}\right)$

(Eq. 17)

An alternative form of (Eq. 17) is obtained by using $de=C_{v}dT$ in (Eq. 11), which gives

$s_2-s_1=\underset{1}{\overset{2}{\int }}{C}_v\frac{dT}{T}-R\ln \left(\frac{{\rho }_2}{{\rho }_1}\right)$

(Eq. 18)

Again, for a calorically perfect gas, we get

$s_2-s_1=C_v\ln \left(\frac{T_2}{T_1}\right)-R\ln \left(\frac{{\rho }_2}{{\rho }_1}\right)$

(Eq. 19)

Isentropic Relations

Adiabatic and reversible processes, i.e., isentropic processes implies $ds=0$ and thus (Eq. 17) reduces to

$\frac{C_p}{R}\ln \left(\frac{T_2}{T_1}\right)=\ln \left(\frac{p_2}{p_1}\right)$

(Eq. 20)

$$\frac{C_p}{R}=\frac{\gamma }{\gamma -1}$$

$$\frac{\gamma }{\gamma -1}\ln \left(\frac{T_2}{T_1}\right)=\ln \left(\frac{p_2}{p_1}\right)\Rightarrow$$

$\frac{p_2}{p_1}={\left(\frac{T_2}{T_1}\right)}^{\gamma /(\gamma -1)}$

(Eq. 21)

In the same way, (Eq. 19) gives

$\frac{{\rho }_2}{{\rho }_1}={\left(\frac{T_2}{T_1}\right)}^{1/(\gamma -1)}$

(Eq. 22)

Eqn. (Eq. 21) and Eqn. (Eq. 22) constitutes the isentropic relations

$\frac{p_2}{p_1}={\left(\frac{{\rho }_2}{{\rho }_1}\right)}^{\gamma }={\left(\frac{T_2}{T_1}\right)}^{\gamma /(\gamma -1)}$

(Eq. 23)

Thermodynamic processes

$ds=C_v\frac{dT}{T}+R\frac{d\nu }{\nu }$

(Eq. 24)

$d\nu =\frac{\nu }{R}ds-C_v\frac{\nu }{RT}dT=\frac{\nu }{R}ds-\frac{C_v}{p}dT$

(Eq. 25)

for an isentropic process ($ds=0$), $d\nu <0$ for positive values of $dT$.

$ds=C_p\frac{dT}{T}-R\frac{dp}{p}$

(Eq. 26)

$dp=-\frac{p}{R}ds+C_p\frac{p}{RT}dT=-\frac{p}{R}ds+C_p\rho dT$

(Eq. 27)

for an isentropic process ($ds=0$), $dp>0$ for positive values of $dT$.

Since $\nu$ decreases with temperature and pressure increases with temperature for an isentropic process, we can see from (Eq. 25) that $d\nu$ will be greater at lower temperatures and thus isochores (lines of constant specific volume) will be closely spaced at low temperatures and more sparse at higher temperatures. For an isochore $dv=0$ which implies

$0=\frac{\nu }{R}(ds-C_v\frac{dT}{T})\Rightarrow \frac{dT}{ds}=\frac{T}{C_v}$

(Eq. 28)

and thus we can see that the slope of an isochore in a $T-s$-diagram is positive and that the slope increases with temperature.

In analogy, we can see that an isobar ($dp=0$) leads to the following relation

$0=\frac{p}{R}(C_p\frac{dT}{T}-ds)\Rightarrow \frac{dT}{ds}=\frac{T}{C_p}$

(Eq. 29)

and consequently isobars will also have a positive slope that increases with temperature in a $T-s$-diagram. Moreover, isobars are less steep than ischores as $C_p>C_v$.

Governing equations

Governing equations on integral form

The governing equations stems from mass conservation, conservation of momentum and conservation of energy

The Continuity Equation

"Mass can be neither created nor destroyed, which implies that mass is conserved"

The net massflow into the control volume $\Omega$ in Fig. \ref{fig:generic:cv} is obtained by integrating mass flux over the control volume surface $\partial \Omega$

$-\underset{\partial \Omega }{\iint }\rho 𝐯\cdot 𝐧dS$

(Eq. 30)

Now, let's consider a small infinitesimal volume $dV$ inside $\Omega$. The mass of $dV$ is $\rho dV$. Thus, the mass enclosed within $\Omega$ can be calculated as

$\underset{\Omega }{\iiint }\rho dV$

(Eq. 31)

The rate of change of mass within $\Omega$ is obtained as

$\frac{d}{dt}\underset{\Omega }{\iiint }\rho dV$

(Eq. 32)

Mass is conserved, which means that the rate of change of mass within $\Omega$ must equal the net flux over the control volume surface.

$\frac{d}{dt}\underset{\Omega }{\iiint }\rho dV=-\underset{\partial \Omega }{\iint }\rho 𝐯\cdot 𝐧dS$

(Eq. 33)

or

$\frac{d}{dt}\underset{\Omega }{\iiint }\rho dV+\underset{\partial \Omega }{\iint }\rho 𝐯\cdot 𝐧dS=0$

(Eq. 34)

which is the integral form of the continuity equation.

The Momentum Equation

"The time rate of change of momentum of a body equals the net force exerted on it"

$\frac{d}{dt}(m𝐯)=𝐅$

(Eq. 35)

What type of forces do we have?

• Body forces acting on the fluid inside $\Omega$
  • gravitation
  • electromagnetic forces
  • Coriolis forces
• Surface forces: pressure forces and shear forces

Body forces inside $\Omega$:

$\underset{\Omega }{\iiint }\rho 𝐟dV$

(Eq. 36)

Surface force on $\partial \Omega$:

$-\underset{\partial \Omega }{\iint }p𝐧dS$

(Eq. 37)

Since we are considering inviscid flow, there are no shear forces and thus we have the net force as

$𝐅=\underset{\Omega }{\iiint }\rho 𝐟dV-\underset{\partial \Omega }{\iint }p𝐧dS$

(Eq. 38)

The fluid flowing through $\Omega$ will carry momentum and the net flow of momentum out from $\Omega$ is calculated as

$\underset{\partial \Omega }{\iint }(\rho 𝐯\cdot 𝐧dS)𝐯=\underset{\partial \Omega }{\iint }(\rho 𝐯\cdot 𝐧)𝐯dS$

(Eq. 39)

Integrated momentum inside $\Omega$

$\underset{\Omega }{\iiint }\rho 𝐯dV$

(Eq. 40)

Rate of change of momentum due to unsteady effects inside $\Omega$

$\frac{d}{dt}\underset{\Omega }{\iiint }\rho 𝐯dV$

(Eq. 41)

Combining the rate of change of momentum, the net momentum flux and the net forces we get

$\frac{d}{dt}\underset{\Omega }{\iiint }\rho 𝐯dV+\underset{\partial \Omega }{\iint }(\rho 𝐯\cdot 𝐧)𝐯dS=\underset{\Omega }{\iiint }\rho 𝐟dV-\underset{\partial \Omega }{\iint }p𝐧dS$

(Eq. 42)

combining the surface integrals, we get

$\frac{d}{dt}\underset{\Omega }{\iiint }\rho 𝐯dV+\underset{\partial \Omega }{\iint }[(\rho 𝐯\cdot 𝐧)𝐯+p𝐧]dS=\underset{\Omega }{\iiint }\rho 𝐟dV$

(Eq. 43)

which is the momentum equation on integral form.

The Energy Equation

"Energy can be neither created nor destroyed; it can only change in form"

$$E_1+E_2=E_3$$

$E_1$ Rate of heat added to the fluid in $\Omega$ from the surroundings

heat transfer

radiation

$E_2$ Rate of work done on the fluid in $\Omega$

$E_3$ Rate of change of energy of the fluid as it flows through $\Omega$

$E_1=\underset{\Omega }{\iiint }\dot{q}\rho dV$

(Eq. 44)

where $\dot{q}$ is the rate of heat added per unit mass

The rate of work done on the fluid in $\Omega$ due to pressure forces is obtained from the pressure force term in the momentum equation.

$E_{2_{pressure}}=-\underset{\partial \Omega }{\iint }(p𝐧dS)\cdot 𝐯=-\underset{\partial \Omega }{\iint }p𝐯\cdot 𝐧dS$

(Eq. 45)

The rate of work done on the fluid in $\Omega$ due to body forces is

$E_{2_{bodyforces}}=\underset{\Omega }{\iiint }(\rho 𝐟dV)\cdot 𝐯=\underset{\Omega }{\iiint }\rho 𝐟\cdot 𝐯dV$

(Eq. 46)

$E_2=E_{2_{pressure}}+E_{2_{bodyforces}}=-\underset{\partial \Omega }{\iint }p𝐯\cdot 𝐧dS+\underset{\Omega }{\iiint }\rho 𝐟\cdot 𝐯dV$

(Eq. 47)

The energy of the fluid per unit mass is the sum of internal energy $e$ (molecular energy) and the kinetic energy $V^2/2$ and the net energy flux over the control volume surface is calculated by the following integral

$\underset{\partial \Omega }{\iint }(\rho 𝐯\cdot 𝐧dS)(e+\frac{V^2}{2})$

(Eq. 48)

Analogous to mass and momentum, the total amount of energy of the fluid in $\Omega$ is calculated as

$\underset{\Omega }{\iiint }\rho (e+\frac{V^2}{2})dV$

(Eq. 49)

The time rate of change of the energy of the fluid in $\Omega$ is obtained as

$\frac{d}{dt}\underset{\Omega }{\iiint }\rho (e+\frac{V^2}{2})dV$

(Eq. 50)

Now, $E_3$ is obtained as the sum of the time rate of change of energy of the fluid in $\Omega$ and the net flux of energy carried by fluid passing the control volume surface.

$E_3=\frac{d}{dt}\underset{\Omega }{\iiint }\rho (e+\frac{V^2}{2})dV+\underset{\partial \Omega }{\iint }(\rho 𝐯\cdot 𝐧dS)(e+\frac{V^2}{2})$

(Eq. 51)

With all elements of the energy equation defined, we are now ready to finally compile the full equation

$\frac{d}{dt}\underset{\Omega }{\iiint }\rho (e+\frac{V^2}{2})dV+\underset{\partial \Omega }{\iint }[\rho (e+\frac{V^2}{2})(𝐯\cdot 𝐧)+p𝐯\cdot 𝐧]dS=$ $\underset{\Omega }{\iiint }\rho 𝐟\cdot 𝐯dV+\underset{\Omega }{\iiint }\dot{q}\rho dV$

(Eq. 52)

The surface integral in the energy equation may be rewritten as

$\underset{\partial \Omega }{\iint }[\rho (e+\frac{V^2}{2})(𝐯\cdot 𝐧)+p𝐯\cdot 𝐧]dS=$ $\underset{\partial \Omega }{\iint }\rho [e+\frac{p}{\rho }+\frac{V^2}{2}](𝐯\cdot 𝐧)dS$

(Eq. 53)

and with the definition of enthalpy $h=e+p/\rho$, we get

$\underset{\partial \Omega }{\iint }\rho [h+\frac{V^2}{2}](𝐯\cdot 𝐧)dS$

(Eq. 54)

Furthermore, introducing total internal energy $e_o$ and total enthalpy $h_o$ defined as

$e_o=e+\frac{1}{2}{V}^2$

(Eq. 55)

and

$h_o=h+\frac{1}{2}{V}^2$

(Eq. 56)

the energy equation is written as

$\frac{d}{dt}\underset{\Omega }{\iiint }\rho e_{o}dV+\underset{\partial \Omega }{\iint }\rho h_o(𝐯\cdot 𝐧)dS=$ $\underset{\Omega }{\iiint }\rho 𝐟\cdot 𝐯dV+\underset{\Omega }{\iiint }\dot{q}\rho dV$

(Eq. 57)

Summary

The integral form of the governing equations for inviscid compressible flow has been derived

Continuity:

$\frac{d}{dt}\underset{\Omega }{\iiint }\rho dV+\underset{\partial \Omega }{\iint }\rho 𝐯\cdot 𝐧dS=0$

Momentum:

$\frac{d}{dt}\underset{\Omega }{\iiint }\rho 𝐯dV+\underset{\partial \Omega }{\iint }[(\rho 𝐯\cdot 𝐧)𝐯+p𝐧]dS=\underset{\Omega }{\iiint }\rho 𝐟dV$

Energy:

$\frac{d}{dt}\underset{\Omega }{\iiint }\rho e_{o}dV+\underset{\partial \Omega }{\iint }\rho h_o(𝐯\cdot 𝐧)dS=$ $\underset{\Omega }{\iiint }\rho 𝐟\cdot 𝐯dV+\underset{\Omega }{\iiint }\dot{q}\rho dV$

Governing equations on differential form

The Differential Equations on Conservation Form

Conservation of Mass

The continuity equation on integral form reads

$\frac{d}{dt}\underset{\Omega }{\iiint }\rho dV+\underset{\partial \Omega }{\iint }\rho 𝐯\cdot 𝐧dS=0$

Apply Gauss's divergence theorem on the surface integral gives

$\underset{\partial \Omega }{\iint }\rho 𝐯\cdot 𝐧dS=\underset{\Omega }{\iiint }\mathrm{\nabla }\cdot (\rho 𝐯)dV$

(Eq. 58)

Also, if $\Omega$ is a fixed control volume

$\frac{d}{dt}\underset{\Omega }{\iiint }\rho dV=\underset{\Omega }{\iiint }\frac{\partial \rho }{\partial t}dV$

(Eq. 59)

The continuity equation can now be written as a single volume integral.

$\underset{\Omega }{\iiint }[\frac{\partial \rho }{\partial t}+\mathrm{\nabla }\cdot (\rho 𝐯)]dV=0$

(Eq. 60)

$\Omega$ is an arbitrary control volume and thus

$\frac{\partial \rho }{\partial t}+\mathrm{\nabla }\cdot (\rho 𝐯)=0$

(Eq. 61)

which is the continuity equation on partial differential form.

Conservation of Momentum

The momentum equation on integral form reads

$\frac{d}{dt}\underset{\Omega }{\iiint }\rho 𝐯dV+\underset{\partial \Omega }{\iint }[(\rho 𝐯\cdot 𝐧)𝐯+p𝐧]dS=\underset{\Omega }{\iiint }\rho 𝐟dV$

As for the continuity equation, the surface integral terms are rewritten as volume integrals using Gauss's divergence theorem.

$\underset{\partial \Omega }{\iint }(\rho 𝐯\cdot 𝐧)𝐯dS=\underset{\Omega }{\iiint }\mathrm{\nabla }\cdot (\rho 𝐯𝐯)dV$

(Eq. 62)

$\underset{\partial \Omega }{\iint }p𝐧dS=\underset{\Omega }{\iiint }\mathrm{\nabla }pdV$

(Eq. 63)

Also, if $\Omega$ is a fixed control volume

$\frac{d}{dt}\underset{\Omega }{\iiint }\rho 𝐯dV=\underset{\Omega }{\iiint }\frac{\partial }{\partial t}(\rho 𝐯)dV$

(Eq. 64)

The momentum equation can now be written as one single volume integral

$\underset{\Omega }{\iiint }[\frac{\partial }{\partial t}(\rho 𝐯)+\mathrm{\nabla }\cdot (\rho 𝐯𝐯)+\mathrm{\nabla }p-\rho 𝐟]dV=0$

(Eq. 65)

$\Omega$ is an arbitrary control volume and thus

$\frac{\partial }{\partial t}(\rho 𝐯)+\mathrm{\nabla }\cdot (\rho 𝐯𝐯)+\mathrm{\nabla }p=\rho 𝐟$

(Eq. 66)

which is the momentum equation on partial differential form

Conservation of Energy

The energy equation on integral form reads

$\frac{d}{dt}\underset{\Omega }{\iiint }\rho e_{o}dV+\underset{\partial \Omega }{\iint }\rho h_o(𝐯\cdot 𝐧)dS=$ $\underset{\Omega }{\iiint }\rho 𝐟\cdot 𝐯dV+\underset{\Omega }{\iiint }\dot{q}\rho dV$

Gauss's divergence theorem applied to the surface integral term in the energy equation gives

$\underset{\partial \Omega }{\iint }\rho h_o(𝐯\cdot 𝐧)dS=\underset{\Omega }{\iiint }\mathrm{\nabla }\cdot (\rho h_o𝐯)dV$

(Eq. 67)

Fixed control volume

$\frac{d}{dt}\underset{\Omega }{\iiint }\rho e_{o}dV=\underset{\Omega }{\iiint }\frac{\partial }{\partial t}(\rho e_o)dV$

(Eq. 68)

The energy equation can now be written as

$\underset{\Omega }{\iiint }[\frac{\partial }{\partial t}(\rho e_o)+\mathrm{\nabla }\cdot (\rho h_o𝐯)-\rho 𝐟\cdot 𝐯-\dot{q}\rho ]dV=0$

(Eq. 69)

$\Omega$ is an arbitrary control volume and thus

$\frac{\partial }{\partial t}(\rho e_o)+\mathrm{\nabla }\cdot (\rho h_o𝐯)=\rho 𝐟\cdot 𝐯+\dot{q}\rho$

(Eq. 70)

which is the energy equation on partial differential form

Summary

The governing equations for compressible inviscid flow on partial differential form:

Continuity:

$\frac{\partial \rho }{\partial t}+\mathrm{\nabla }\cdot (\rho 𝐯)=0$

Momentum:

$\frac{\partial }{\partial t}(\rho 𝐯)+\mathrm{\nabla }\cdot (\rho 𝐯𝐯)+\mathrm{\nabla }p=\rho 𝐟$

Energy:

$\frac{\partial }{\partial t}(\rho e_o)+\mathrm{\nabla }\cdot (\rho h_o𝐯)=\rho 𝐟\cdot 𝐯+\dot{q}\rho$

The Differential Equations on Non-Conservation Form

The Substantial Derivative

The substantial derivative operator is defined as

$\frac{D}{Dt}=\frac{\partial }{\partial t}+𝐯\cdot \mathrm{\nabla }$

(Eq. 71)

where the first term of the right hand side is the local derivative and the second term is the convective derivative.

Conservation of Mass

If we apply the substantial derivative operator to density we get

$\frac{D\rho }{Dt}=\frac{\partial \rho }{\partial t}+𝐯\cdot \mathrm{\nabla }\rho$

(Eq. 72)

From before we have the continuity equation on differential form as

$\frac{\partial \rho }{\partial t}+\mathrm{\nabla }\cdot (\rho 𝐯)=0$

(Eq. 73)

which can be rewritten as

$\frac{\partial \rho }{\partial t}+\rho (\mathrm{\nabla }\cdot 𝐯)+𝐯\cdot \mathrm{\nabla }\rho =0$

(Eq. 74)

and thus

$\frac{D\rho }{Dt}+\rho (\mathrm{\nabla }\cdot 𝐯)=0$

(Eq. 75)

Eq. 75 says that the mass of a fluid element with a fixed set of fluid particles is constant as the element moves in space.

Conservation of Momentum

We start from the momentum equation on differential form derived above

$\frac{\partial }{\partial t}(\rho 𝐯)+\mathrm{\nabla }\cdot (\rho 𝐯𝐯)+\mathrm{\nabla }p=\rho 𝐟$

(Eq. 76)

Expanding the first and the second terms gives

$\rho \frac{\partial 𝐯}{\partial t}+𝐯\frac{\partial \rho }{\partial t}+\rho 𝐯\cdot \mathrm{\nabla }𝐯+𝐯(\mathrm{\nabla }\cdot \rho 𝐯)+\mathrm{\nabla }p=\rho 𝐟$

(Eq. 77)

Collecting terms, we can identify the substantial derivative operator applied to the velocity vector and the continuity equation.

$\rho \underset{=\frac{D𝐯}{Dt}}{\underbrace{[\frac{\partial 𝐯}{\partial t}+𝐯\cdot \mathrm{\nabla }𝐯]}}+𝐯\underset{=0}{\underbrace{[\frac{\partial \rho }{\partial t}+\mathrm{\nabla }\cdot \rho 𝐯]}}+\mathrm{\nabla }p=\rho 𝐟$

(Eq. 78)

which gives us the non-conservation form of the momentum equation

$\frac{D𝐯}{Dt}+\frac{1}{\rho }\mathrm{\nabla }p=𝐟$

(Eq. 79)

Conservation of Energy

The last equation on non-conservation differential form is the energy equation. We start by rewriting the energy equation on differential form (Eq. 70), repeated here for convenience

$\frac{\partial }{\partial t}(\rho e_o)+\mathrm{\nabla }\cdot (\rho h_o𝐯)=\rho 𝐟\cdot 𝐯+\dot{q}\rho$

Total enthalpy, $h_o$, is replaced with total energy, $e_o$

$h_o=e_o+\frac{p}{\rho }$

(Eq. 80)

which gives

$\frac{\partial }{\partial t}(\rho e_o)+\mathrm{\nabla }\cdot (\rho e_o𝐯)+\mathrm{\nabla }\cdot (p𝐯)=\rho 𝐟\cdot 𝐯+\dot{q}\rho$

(Eq. 81)

Expanding the two first terms as

$\rho \frac{\partial e_o}{\partial t}+e_o\frac{\partial \rho }{\partial t}+\rho 𝐯\cdot \mathrm{\nabla }{e}_o+e_o\mathrm{\nabla }\cdot (\rho 𝐯)+\mathrm{\nabla }\cdot (p𝐯)=$ $=\rho 𝐟\cdot 𝐯+\dot{q}\rho$

(Eq. 82)

Collecting terms, we can identify the substantial derivative operator applied on total energy, $D{e}_o/Dt$ and the continuity equation

$\rho \underset{=\frac{D{e}_o}{Dt}}{\underbrace{[\frac{\partial e_o}{\partial t}+𝐯\cdot \mathrm{\nabla }{e}_o]}}+e_o\underset{=0}{\underbrace{[\frac{\partial \rho }{\partial t}+\mathrm{\nabla }\cdot (\rho 𝐯)]}}+\mathrm{\nabla }\cdot (p𝐯)=\rho 𝐟\cdot 𝐯+\dot{q}\rho$

(Eq. 83)

and thus we end up with the energy equation on non-conservation differential form

$\rho \frac{D{e}_o}{Dt}+\mathrm{\nabla }\cdot (p𝐯)=\rho 𝐟\cdot 𝐯+\dot{q}\rho$

(Eq. 84)

Summary

Continuity:

$\frac{D\rho }{Dt}+\rho (\mathrm{\nabla }\cdot 𝐯)=0$

Momentum:

$\frac{D𝐯}{Dt}+\frac{1}{\rho }\mathrm{\nabla }p=𝐟$

Energy:

$\rho \frac{D{e}_o}{Dt}+\mathrm{\nabla }\cdot (p𝐯)=\rho 𝐟\cdot 𝐯+\dot{q}\rho$

Alternative Forms of the Energy Equation

Internal Energy Formulation

Total internal energy is defined as

$e_o=e+\frac{1}{2}𝐯\cdot 𝐯$

(Eq. 85)

Inserted in Eq. 84, this gives

$\rho \frac{De}{Dt}+\rho 𝐯\cdot \frac{D𝐯}{Dt}+\mathrm{\nabla }\cdot (p𝐯)=\rho 𝐟\cdot 𝐯+\dot{q}\rho$

(Eq. 86)

Now, let's replace the substantial derivative $D𝐯/Dt$ using the momentum equation on non-conservation form (Eq. 79).

$\rho \frac{De}{Dt}-𝐯\cdot \mathrm{\nabla }p+\cancel{\rho 𝐟\cdot 𝐯}+\mathrm{\nabla }\cdot (p𝐯)=\cancel{\rho 𝐟\cdot 𝐯}+\dot{q}\rho$

(Eq. 87)

Now, expand the term $\mathrm{\nabla }\cdot (p𝐯)$ gives

$\rho \frac{De}{Dt}\cancel{-𝐯\cdot \mathrm{\nabla }p}+\cancel{𝐯\cdot \mathrm{\nabla }p}+p(\mathrm{\nabla }\cdot 𝐯)=\dot{q}\rho \Rightarrow$ $\Rightarrow \rho \frac{De}{Dt}+p(\mathrm{\nabla }\cdot 𝐯)=\dot{q}\rho$

(Eq. 88)

Divide by $\rho$

$\frac{De}{Dt}+\frac{p}{\rho }(\mathrm{\nabla }\cdot 𝐯)=\dot{q}$

(Eq. 89)

Conservation of mass gives

$\frac{D\rho }{Dt}+\rho (\mathrm{\nabla }\cdot 𝐯)=0\Rightarrow \mathrm{\nabla }\cdot 𝐯=-\frac{1}{\rho }\frac{D\rho }{Dt}$

(Eq. 90)

Insert in Eq. 89

$\frac{De}{Dt}-\frac{p}{{\rho }^2}\frac{D\rho }{Dt}=\dot{q}\Rightarrow \frac{De}{Dt}+p\frac{D}{Dt}\left(\frac{1}{\rho }\right)=\dot{q}$

(Eq. 91)

$\frac{De}{Dt}+p\frac{D\nu }{Dt}=\dot{q}$

(Eq. 92)

Compare with the first law of thermodynamics: $de=\delta q-\delta w$

Enthalpy Formulation

$h=e+\frac{p}{\rho }\Rightarrow \frac{Dh}{Dt}=\frac{De}{Dt}+\frac{1}{\rho }\frac{Dp}{Dt}+p\frac{D}{Dt}\left(\frac{1}{\rho }\right)$

(Eq. 93)

with $De/Dt$ from Eq. 89

$\frac{Dh}{Dt}=\dot{q}-\cancel{p\frac{D}{Dt}\left(\frac{1}{\rho }\right)}+\frac{1}{\rho }\frac{Dp}{Dt}+\cancel{p\frac{D}{Dt}\left(\frac{1}{\rho }\right)}$

(Eq. 94)

$\frac{Dh}{Dt}=\dot{q}+\frac{1}{\rho }\frac{Dp}{Dt}$

(Eq. 95)

Total Enthalpy Formulation

$h_o=h+\frac{1}{2}𝐯𝐯\Rightarrow \frac{D{h}_o}{Dt}=\frac{Dh}{Dt}+𝐯\cdot \frac{D𝐯}{Dt}$

(Eq. 96)

From the momentum equation (Eq. 79)

$\frac{D𝐯}{Dt}=𝐟-\frac{1}{\rho }\mathrm{\nabla }p$

(Eq. 97)

which gives

$\frac{D{h}_o}{Dt}=\frac{Dh}{Dt}+𝐯\cdot 𝐟-\frac{1}{\rho }𝐯\cdot \mathrm{\nabla }p$

(Eq. 98)

Inserting $Dh/Dt$ from Eq. 95 gives

$\frac{D{h}_o}{Dt}=\dot{q}+\frac{1}{\rho }\frac{Dp}{Dt}+𝐯\cdot 𝐟-\frac{1}{\rho }𝐯\cdot \mathrm{\nabla }p=$ $=\frac{1}{\rho }[\frac{Dp}{Dt}-𝐯\cdot \mathrm{\nabla }p]+\dot{q}+𝐯\cdot 𝐟$

(Eq. 99)

The substantial derivative operator applied to pressure

$\frac{Dp}{Dt}=\frac{\partial p}{\partial t}+𝐯\cdot \mathrm{\nabla }p$

(Eq. 100)

and thus

$\frac{Dp}{Dt}-𝐯\cdot \mathrm{\nabla }p=\frac{\partial p}{\partial t}$

(Eq. 101)

which gives

$\frac{D{h}_o}{Dt}=\frac{1}{\rho }\frac{\partial p}{\partial t}+\dot{q}+𝐯\cdot 𝐟$

(Eq. 102)

If we assume adiabatic flow without body forces

$\frac{D{h}_o}{Dt}=\frac{1}{\rho }\frac{\partial p}{\partial t}$

(Eq. 103)

If we further assume the flow to be steady state we get

$\frac{D{h}_o}{Dt}=0$

(Eq. 104)

This means that in a steady-state adiabatic flow without body forces, total enthalpy is constant along a streamline.

The entropy equation

From the second law of thermodynamics

$\frac{De}{Dt}=T\frac{Ds}{Dt}-p\frac{D}{Dt}\left(\frac{1}{\rho }\right)$

(Eq. 105)

From the energy equation on differential non-conservation form internal energy formulation

$\frac{De}{Dt}=\dot{q}-\frac{p}{\rho }(\mathrm{\nabla }\cdot 𝐯)$

(Eq. 106)

The continuity equation on differential non-conservation form

$\frac{D\rho }{Dt}+\rho (\mathrm{\nabla }\cdot 𝐯)=0\Rightarrow \mathrm{\nabla }\cdot 𝐯=-\frac{1}{\rho }\frac{D\rho }{Dt}$

(Eq. 107)

and thus

$\frac{De}{Dt}=\dot{q}+\frac{p}{{\rho }^2}\frac{D\rho }{Dt}$

(Eq. 108)

$\frac{D\rho }{Dt}=-\frac{1}{{\nu }^2}\frac{D\nu }{Dt}$

(Eq. 109)

$\rho \frac{De}{Dt}=\rho \dot{q}-\frac{p}{\rho {\nu }^2}\frac{D\nu }{Dt}=\rho \dot{q}-\rho p\frac{D\nu }{Dt}$

(Eq. 110)

$\rho [\frac{De}{Dt}+p\frac{D\nu }{Dt}-\dot{q}]=0\Rightarrow \frac{De}{Dt}=\dot{q}-p\frac{D\nu }{Dt}$

(Eq. 111)

Insert $De/Dt$ in Eqn. \ref{eq:second:law}

$\dot{q}-p\frac{D}{Dt}\left(\frac{1}{\rho }\right)=T\frac{Ds}{Dt}-p\frac{D}{Dt}\left(\frac{1}{\rho }\right)\Rightarrow$

(Eq. 112)

$T\frac{Ds}{Dt}=-\dot{q}$

(Eq. 113)

Adiabatic flow:

$T\frac{Ds}{Dt}=0$

(Eq. 114)

In an adiabatic, steady-state, inviscid flow, entropy is constant along a streamline.

Crocco's equation

The momentum equation without body forces

$\rho \frac{D𝐯}{Dt}=-\mathrm{\nabla }p$

(Eq. 115)

Expanding the substantial derivative

$\rho \frac{\partial 𝐯}{\partial t}+\rho 𝐯\cdot \mathrm{\nabla }𝐯=-\mathrm{\nabla }p$

(Eq. 116)

The first and second law of thermodynamics gives

$T\mathrm{\nabla }s=\mathrm{\nabla }h-\frac{\mathrm{\nabla }p}{\rho }$

(Eq. 117)

Insert $\mathrm{\nabla }p$ from the momentum equation

$T\mathrm{\nabla }s=\mathrm{\nabla }h+\frac{\partial 𝐯}{\partial t}+𝐯\cdot \mathrm{\nabla }𝐯$

(Eq. 118)

Definition of total enthalpy ($h_o$)

$h_o=h+\frac{1}{2}𝐯\cdot 𝐯\Rightarrow \mathrm{\nabla }h=\mathrm{\nabla }{h}_o-\mathrm{\nabla }(\frac{1}{2}𝐯\cdot 𝐯)$

(Eq. 119)

The last term can be rewritten as

$\mathrm{\nabla }(\frac{1}{2}𝐯\cdot 𝐯)=𝐯\times (\mathrm{\nabla }\times 𝐯)+𝐯\cdot \mathrm{\nabla }𝐯$

(Eq. 120)

which gives

$\mathrm{\nabla }h=\mathrm{\nabla }{h}_o-𝐯\times (\mathrm{\nabla }\times 𝐯)-𝐯\cdot \mathrm{\nabla }𝐯$

(Eq. 121)

Insert $\mathrm{\nabla }h$ in the entropy equation gives

$T\mathrm{\nabla }s=\mathrm{\nabla }{h}_o-𝐯\times (\mathrm{\nabla }\times 𝐯)-\cancel{𝐯\cdot \mathrm{\nabla }𝐯}+\frac{\partial 𝐯}{\partial t}+\cancel{𝐯\cdot \mathrm{\nabla }𝐯}$

(Eq. 122)

$T\mathrm{\nabla }s=\mathrm{\nabla }{h}_o-𝐯\times (\mathrm{\nabla }\times 𝐯)+\frac{\partial 𝐯}{\partial t}$

(Eq. 123)

One-dimensional inviscid flow

Reference flow states

Stagnation Flow Properties

$h_o=h+\frac{1}{2}{u}^2\Rightarrow Cp{T}_o=C_{p}T+\frac{1}{2}{u}^2\Rightarrow \frac{T_o}{T}=1+\frac{1}{2C_p}{M}^2\gamma RT=1+\frac{\gamma -1}{2}{M}^2$

(Eq. 124)

Acoustic waves

In Fig. \ref{fig:soundwave}, station 1 represents the flow state ahead of the sound wave and station 2 the flow state behind the sound wave. Set up the continuity equation for one-dimensional flows between 1 and 2. If we could change frame of reference and follow the sound wave, we would see fluid approaching the wave with the propagation speed of the wave, $a$, and behind the wave, the fluid would have a slightly modified speed, $a+da$. There would also be a slight in all other flow properties. Let's apply the one-dimensional continuity equation between station 1 and station 2.

${\rho }_1{u}_1={\rho }_2{u}_2$

(Eq. 125)

$\rho a=(\rho +d\rho )(a+da)$

(Eq. 126)

$\cancel{\rho a}=\cancel{\rho a}+\rho da+ad\rho +\underset{\sim 0}{\underbrace{d\rho da}}\Rightarrow$

(Eq. 127)

$a=-\rho \frac{da}{d\rho }$}

(Eq. 128)

The one-dimensional momentum equation between station 1 and station 2 gives

${\rho }_1{u}_1^2+p_1={\rho }_2{u}_2^2+p_2$

(Eq. 129)

$\rho a^2+p=(\rho +d\rho )(a+da{)}^2+(p+dp)$

(Eq. 130)

$\cancel{\rho a^2}+\cancel{p}=\cancel{\rho a^2}+2\rho ada+\underset{\sim 0}{\underbrace{\rho d{a}^2}}+d\rho a^2+\underset{\sim 0}{\underbrace{2d\rho ada}}+\underset{\sim 0}{\underbrace{d\rho d{a}^2}}+\cancel{p}+dp\Rightarrow$

(Eq. 131)

$dp=-2\rho ada-d\rho a^2\Rightarrow$

(Eq. 132)

$da=-\frac{dp+d\rho a^2}{2\rho a}=-\frac{d\rho }{2a\rho }(\frac{dp}{d\rho }+a^2)\Rightarrow$

(Eq. 133)

$\frac{da}{d\rho }=-\frac{1}{2a\rho }(\frac{dp}{d\rho }+a^2)$

(Eq. 134)

Eq. 134 in Eq. 128 gives

$a=\frac{1}{2a}(\frac{dp}{d\rho }+a^2)\Rightarrow$

(Eq. 135)

$a^2=\frac{dp}{d\rho }$

(Eq. 136)

Sound wave:

• gradients are small
• irreversible (dissipative effects are negligible)
• no heat addition

Thus, the change of flow properties as the sound wave passes can be assumed to be an isentropic process

$a^2={\left(\frac{dp}{d\rho }\right)}_s$

(Eq. 137)

$a=\sqrt{{\left(\frac{dp}{d\rho }\right)}_s}=\sqrt{\frac{1}{\rho {\tau }_s}}$

(Eq. 138)

where ${\tau }_s$ is the compressibility of the gas. Eq. 137 is valid for all gases. It can be seen from the equation, that truly incompressible flow (${\tau }_s=0$) would imply infinite speed of sound.

Since the process is isentropic, we can use the isentropic relations if we also assume the gas to be calorically perfect

$\frac{p_2}{p_1}={\left(\frac{{\rho }_2}{{\rho }_1}\right)}^{\gamma }\Rightarrow p=C{\rho }^{\gamma }$

(Eq. 139)

$a^2={\left(\frac{dp}{d\rho }\right)}_s=\gamma C{\rho }^{\gamma -1}=\gamma \underset{=p}{\underbrace{\left[C{\rho }^{\gamma }\right]}}{\rho }^{-1}=\frac{\gamma p}{\rho }\Rightarrow$

(Eq. 140)

$a=\sqrt{\frac{\gamma p}{\rho }}$

(Eq. 141)

or

$a=\sqrt{\gamma RT}$

(Eq. 142)

From the relation above, it is obvious that the local speed of sound is related to the temperature of the flow, which in turn is a measure of the motion of elementary particles (atoms and/or molecules) of the fluid at a specific location. This stems from the fact that sound waves are propagated via interaction of these elementary particles. Since information in a flow is propagated via molecular interaction the relation between the speed at which this information is conveyed and the speed of the flow has important physical implications. Figure~\ref{fig:speed:of:sound} compares three sound wave patterns generated by a a beacon. In the left picture, the sound transmitter is stationary and thus the acoustic waves are centered around the transmitter. In the middle image, the transmitter is moving to the left at a speed less than the speed of sound and thus the transmitter will always be within all sound wave circles but it will be off-centered with a bias in the direction that the transmitter is moving. In the right image the transmitter is moving faster than the speed of sound and thus it will always be located outside of all acoustic waves. In a supersonic flow, no information can travel upstream and therefore there is no way for the flow to adjust to downstream obstacles. This is compensated for by the introduction of shocks in the flow. Over a shock flow properties changes discontinuity. An example is given in Figure~\ref{fig:supersonic:flow}.

Shock waves

The starting point is to set up the governing equations for one-dimensional steady compressible flow over a control volume enclosing the normal shock (Fig. \ref{fig:shock:cv}).

continuity:

${\rho }_1{u}_1={\rho }_2{u}_2$

(Eq. 143)

momentum:

${\rho }_1{u}_1^2+p_1={\rho }_2{u}_2^2+p_2$

(Eq. 144)

energy:

$h_1+\frac{1}{2}{u}_1^2=h_2+\frac{1}{2}{u}_2^2$

(Eq. 145)

Divide the momentum equation by ${\rho }_1{u}_1$

$\frac{1}{{\rho }_1{u}_1}({\rho }_1{u}_1^2+p_1)=\frac{1}{{\rho }_1{u}_1}({\rho }_2{u}_2^2+p_2)=\{{\rho }_1{u}_1={\rho }_2{u}_2\}=\frac{1}{{\rho }_2{u}_2}({\rho }_2{u}_2^2+p_2)\Rightarrow$

(Eq. 146)

$\frac{p_1}{{\rho }_1{u}_1}-\frac{p_2}{{\rho }_2{u}_2}=u_2-u_1$

(Eq. 147)

For a calorically perfect gas $a=\sqrt{\gamma p/\rho }$, which if implemented in Eqn. \ref{eq:governing:mom:b} gives

$\frac{a_1^2}{\gamma u_1}-\frac{a_2^2}{\gamma u_2}=u_2-u_1$

(Eq. 148)

The energy equation (Eqn. \ref{eq:governing:energy}) with $h=C_{p}T$

$C_p{T}_1+\frac{1}{2}{u}_1^2=C_p{T}_2+\frac{1}{2}{u}_2^2$

(Eq. 149)

Replacing $C_p$ with $\gamma R/(\gamma -1)$ gives

$\frac{\gamma R{T}_1}{\gamma -1}+\frac{1}{2}{u}_1^2=\frac{\gamma R{T}_2}{\gamma -1}+\frac{1}{2}{u}_2^2$

(Eq. 150)

With $a=\sqrt{\gamma RT}$ this becomes

$\frac{a_1^2}{\gamma -1}+\frac{1}{2}{u}_1^2=\frac{a_2^2}{\gamma -1}+\frac{1}{2}{u}_2^2$

(Eq. 151)

Eqn. \ref{eq:governing:energy:d} can be set up between any two points in the flow. Specifically, we can use the relation to relate the flow velocity, $u$, and speed of sound, $a$, in any point to the corresponding flow properties at sonic conditions ($u=a=a^{*}$).

$\frac{a^2}{\gamma -1}+\frac{1}{2}{u}^2=\frac{\gamma +1}{2(\gamma -1)}{a^{*}}^2$

(Eq. 152)

If Eqn. \ref{eq:governing:energy:e} is evaluated in locations 1 and 2, we get

$\begin{array}{cc}& a_1^2=\frac{\gamma +1}{2}{a^{*}}^2-\frac{\gamma -1}{2}{u}_1^2\\ & a_2^2=\frac{\gamma +1}{2}{a^{*}}^2-\frac{\gamma -1}{2}{u}_2^2\end{array}$

(Eq. 153)

Since the change in flow conditions over the shock is adiabatic (no heat is added inside the shock), critical properties will be constant over the shock. Especially $a^{*}$ will be constant.

Eqn. \ref{eq:governing:energy:f} inserted in \ref{eq:governing:mom:c} gives\\

$\frac{1}{\gamma u_1}(\frac{\gamma +1}{2}{a^{*}}^2-\frac{\gamma -1}{2}{u}_1^2)-\frac{1}{\gamma u_2}(\frac{\gamma +1}{2}{a^{*}}^2-\frac{\gamma -1}{2}{u}_2^2)=u_2-u_1\Rightarrow$

(Eq. 154)

$\left(\frac{\gamma +1}{2\gamma }\right){a^{*}}^2(\frac{1}{u_1}-\frac{1}{u_2})=\left(\frac{\gamma +1}{2\gamma }\right)(u_2-u_1)\Rightarrow$

(Eq. 155)

${a^{*}}^2(\frac{1}{u_1}-\frac{1}{u_2})=(u_2-u_1)\Rightarrow$

(Eq. 156)

${a^{*}}^2(\frac{u_2}{u_1{u}_2}-\frac{u_1}{u_1{u}_2})=(u_2-u_1)\Rightarrow$

(Eq. 157)

$\frac{1}{u_1{u}_2}{a^{*}}^2(u_2-u_1)=(u_2-u_1)\Rightarrow$

(Eq. 158)

${a^{*}}^2=u_1{u}_2$

(Eq. 159)

Eqn. \ref{eq:prandtl} is sometimes referred to as the Prandtl relation. Divide the Prandtl relation by ${a^{*}}^2$ on both sides gives

$1=\frac{u_1}{a^{*}}\frac{u_2}{a^{*}}=M_1^{*}{M}_2^{*}$

(Eq. 160)

or

$M_2^{*}=\frac{1}{M_1^{*}}$

(Eq. 161)

The relation between $M^{*}$ and $M$ is given by

${M^{*}}^2=\frac{(\gamma +1)M^2}{2+(\gamma -1)M^2}$

(Eq. 162)

from which is can be seen that $M^{*}$ will follow the Mach number $M$ in the sense that

• $M=1\Rightarrow M^{*}=1$
• $M<1\Rightarrow M^{*}<1$
• $M>1\Rightarrow M^{*}>1$

Eqn. \ref{eq:MachStar} inserted in Eqn. \ref{eq:NormalMach} gives

$\frac{(\gamma +1)M_1^2}{2+(\gamma -1)M_1^2}=\frac{2+(\gamma -1)M_2^2}{(\gamma +1)M_2^2}$

(Eq. 163)

$M_2^2=\frac{1+[(\gamma -1)/2]M_1^2}{\gamma M_1^2-(\gamma -1)/2}$

(Eq. 164)

The Mach number relations above effectively show that if the Mach number upstream of the shock is greater than one, the downstream Mach number must be less than one and vice versa. We can also see that a sonic upstream flow $M_1=1.0$ gives sonic flow downstream of the shock. So, apparently the relation as such holds for both supersonic and subsonic upstream flow mathematically. The question is if it is also physically correct. For a supersonic upstream flow we will get a discontinuous compression and if the flow upstream of the control volume is subsonic we will instead get a discontinuous expansion inside the control volume but, again, is this physically correct? We will get the answer by analyzing the entropy change over the control volume.

Analyzing the energy equation and the second law of thermodynamics shows that there is a direct relation between entropy increase and total pressure drop.

$s_2-s_1=C_p\ln \frac{T_2}{T_1}-R\ln \frac{p_2}{p_1}$

(Eq. 165)

$s_2-s_1=C_p\ln \frac{T_2}{T_{o_2}}\frac{T_{o_1}}{T_1}\frac{T_{o_2}}{T_{o_1}}-R\ln \frac{p_2}{p_{o_2}}\frac{p_{o_1}}{p_1}\frac{p_{o_2}}{p_{o_1}}$

(Eq. 166)

using the isentropic relations we get

$s_2-s_1=C_p\ln \frac{T_{o_2}}{T_{o_1}}-R\ln \frac{p_{o_2}}{p_{o_1}}$

(Eq. 167)

and since the process is adiabatic and thus $T_{o_2}=T_{o_1}$ the change in entropy is directly related to the change in total pressure as

$s_2-s_1=-R\ln \frac{p_{o_2}}{p_{o_1}}$

(Eq. 168)

or

$\frac{p_{o_2}}{p_{o_1}}=e^{-(s_2-s_1)/R}$

(Eq. 169)

Figure~\ref{fig:shock:entropy} shows the entropy change over a normal shock. As can be seen in the figure, a subsonic upstream Mach number leads to a reduction of entropy, which once and for all rules out all such solutions as non-physical and thus the question about the upstream conditions can now be considered answered. This in turn implies that the Mach number downstream of a normal shock will always be subsonic, which can be seen in Fig~\ref{fig:shock:downstream:Mach} below.

By rewriting the right-hand side of Eqn.\ref{eq:NormalMach:b}, it is easy to realize that the downstream Mach number $M_2$ approaches a finite value for large values of the upstream Mach number, $M_1$.

${M_2^2|}_{M_1\to \mathrm{\infty }}={\frac{2/M_1^2+(\gamma -1)}{2\gamma -(\gamma -1)/M_1^2}|}_{M_1\to \mathrm{\infty }}=\frac{\gamma -1}{2\gamma }$

(Eq. 170)

Normal-shock relations

Rewriting the continuity equation (Eqn. \ref{eq:governing:cont})

$\frac{{\rho }_2}{{\rho }_1}=\frac{u_1}{u_2}=\frac{u_1^2}{u_1{u}_2}=\{{a^{*}}^2=u_1{u}_2\}=\frac{u_1^2}{{a^{*}}^2}={M_1^{*}}^2$

(Eq. 171)

Eqn. \ref{eq:MachStar} in Eqn. \ref{eq:Normal:density:a} gives

$\frac{{\rho }_2}{{\rho }_1}=\frac{(\gamma +1)M_1^2}{2+(\gamma -1)M_1^2}$

(Eq. 172)

To get a corresponding relation for the pressure ratio over the shock, we go back to the momentum equation (Eqn. \ref{eq:governing:mom})

$p_2-p_1={\rho }_1{u}_1^2-{\rho }_2{u}_2^2=\{{\rho }_1{u}_1={\rho }_2{u}_1\}={\rho }_1{u}_1(u_1-u_2)={\rho }_1{u}_1^2(1-\frac{u_2}{u_1})$

(Eq. 173)

$\frac{p_2-p_1}{p_1}=\frac{{\rho }_1{u}_1^2}{p_1}(1-\frac{u_2}{u_1})=\{a_1=\sqrt{\frac{\gamma p_1}{{\rho }_1}}\}=\gamma \frac{u_1^2}{a_1^2}(1-\frac{u_2}{u_1})=\gamma M_1^2(1-\frac{u_2}{u_1})$

(Eq. 174)

$\frac{p_2}{p_1}-1=\gamma M_1^2(1-\frac{u_2}{u_1})=\{\frac{u_2}{u_1}=\frac{{\rho }_1}{{\rho }_2}\}=\gamma M_1^2(1-\frac{2+(\gamma -1)M_1^2}{(\gamma +1)M_1^2})$

(Eq. 175)

$\frac{p_2}{p_1}=1+\frac{2\gamma }{\gamma +1}(M_1^2-1)$

(Eq. 176)

Figure~\ref{fig:shock:pressure:ratio} shows that the pressure must increase over the shock due to the fact that, based on the discussion above, the upstream Mach number must be greater than one and thus the shock is a discontinuous compression process.

The temperature ratio over the shock can be obtained using the already derived relations for pressure ratio and density ratio together with the equation of state $p=\rho RT$

$\frac{T_2}{T_1}=\left(\frac{p_2}{p_1}\right)\left(\frac{{\rho }_1}{{\rho }_2}\right)$

(Eq. 177)

$\frac{T_2}{T_1}=[1+\frac{2\gamma }{\gamma +1}(M_1^2-1)]\left[\frac{(\gamma +1)M_1^2}{2+(\gamma -1)M_1^2}\right]$

(Eq. 178)

Figure~\ref{fig:normal:shock:relations} below shows how different flow properties change over a normal shock as a function of upstream Mach number.

Now, one question remains. How come that we by analyzing the control volume using the upstream and downstream states get the normal shock relations. There is no way that the governing equations could have known about the fact that we assumed that there would be a shock inside of the control volume, or is it? The answer is that we have assumed that there will be a change in flow properties from upstream to downstream. We have further assumed that the flow is adiabatic (we are using the adiabatic energy equation) so there is no heat exchange. We are, however, allowing for irreversibilities in the flow. The only way to accomplish a change in flow properties under those constraints is a formation of a normal shock (a discontinuity in flow properties - a sudden flow compression) between station 1 and station 2.

The Hugoniot equation

The Hugoniot equation is an alternative normal shock relation based on thermodynamic quantities only. It is derived from the governing equations and relates the change in energy to the change in pressure and specific volume. The starting point of the derivation of the Hugoniot equation is the governing equations (Eqns~\ref{eq:governing:cont} - \ref{eq:governing:energy}).

The continuity equation is rewritten and inserted into the momentum equation

$u_1=\left(\frac{{\rho }_2}{{\rho }_1}\right)u_2$

(Eq. 179)

Replace $u_1$ in Eqn. \ref{eq:governing:mom} using Eqn. \ref{eq:governing:cont:b}

${\rho }_1{\left(\frac{{\rho }_2}{{\rho }_1}\right)}^2{u}_2^2+p_1={\rho }_2{u}_2^2+p_2$

(Eq. 180)

$u_2^2({\rho }_1{\left(\frac{{\rho }_2}{{\rho }_1}\right)}^2-{\rho }_2)=(p_2-p_1)$

(Eq. 181)

$u_2^2\left(\left(\frac{{\rho }_2}{{\rho }_1}\right)({\rho }_2-{\rho }_1)\right)=(p_2-p_1)$

(Eq. 182)

$u_2^2=\left(\frac{{\rho }_1}{{\rho }_2}\right)\frac{p_2-p_1}{{\rho }_2-{\rho }_1}$

(Eq. 183)

Eqn. \ref{eq:governing:cont:b} and \ref{eq:governing:mom:b} gives

$u_1^2=\left(\frac{{\rho }_2}{{\rho }_1}\right)\frac{p_2-p_1}{{\rho }_2-{\rho }_1}$

(Eq. 184)

Eqn. \ref{eq:governing:mom:b} and Eqn. \ref{eq:governing:mom:c} inserted in the energy equation (Eqn. \ref{eq:governing:energy}) gives

$h_1+\frac{1}{2}\left(\frac{{\rho }_2}{{\rho }_1}\right)\left(\frac{p_2-p_1}{{\rho }_2-{\rho }_1}\right)=h_2+\frac{1}{2}\left(\frac{{\rho }_1}{{\rho }_2}\right)\left(\frac{p_2-p_1}{{\rho }_2-{\rho }_1}\right)$

(Eq. 185)

$h_2-h_1=\frac{p_2-p_1}{2}[\left(\frac{{\rho }_2}{{\rho }_1}\right)\left(\frac{1}{{\rho }_2-{\rho }_1}\right)-\left(\frac{{\rho }_1}{{\rho }_2}\right)\left(\frac{1}{{\rho }_2-{\rho }_1}\right)]$

(Eq. 186)

$h_2-h_1=\frac{p_2-p_1}{2}\left[\frac{{\rho }_2^2-{\rho }_1^2}{{\rho }_1{\rho }_2({\rho }_2-{\rho }_1)}\right]=\frac{p_2-p_1}{2}\left[\frac{{\rho }_2+{\rho }_1}{{\rho }_1{\rho }_2}\right]$

(Eq. 187)

$h_2-h_1=\frac{p_2-p_1}{2}(\frac{1}{{\rho }_1}+\frac{1}{{\rho }_2})$

(Eq. 188)

Now, replacing the enthalpies with internal energies using $h=e+p/\rho$ gives

$e_2-e_1=\frac{p_1}{{\rho }_1}-\frac{p_2}{{\rho }_2}+\frac{p_2-p_1}{2}(\frac{1}{{\rho }_1}+\frac{1}{{\rho }_2})$

(Eq. 189)

which after some rewriting becomes the Hugoniot equation

$e_2-e_1=\frac{p_2+p_1}{2}(\frac{1}{{\rho }_1}-\frac{1}{{\rho }_2})=\frac{p_2+p_1}{2}({\nu }_1-{\nu }_2)$

(Eq. 190)

To give an idea about how the normal shock relates to an isentropic compression (a flow compression process without losses) the change in flow density as a function of pressure ratio is compared in Figure~\ref{fig:normal:shock:compression:vs:isentropic}. One can see that the normal-shock compression is more effective but less efficient than the corresponding isentropic process.

Introducing $C$ as the massflow per unit area (which is a constant)

${\rho }_1{u}_1={\rho }_2{u}_2=C$

(Eq. 191)

Inserted into the momentum equation this gives

$p_1+\frac{C^2}{{\rho }_1}=p_2+\frac{C^2}{{\rho }_2}$

(Eq. 192)

or

$\frac{p_2-p_1}{{\nu }_2-{\nu }_1}=-C^2$

(Eq. 193)

which implies that all possible solutions to the governing equations must be located on a line in $p\nu$-space (the so-called Rayleigh line). If we add the Hugoniot relation to this we will find that there are two possible solutions, the upstream condition and the condition downstream of the normal shock and the flow cannot be in any of the intermediate stages. The normal-process is a so-called wave solution to the governing equations where the flow state must jump directly from one flow state to another without passing the intermediate conditions. If we add heat or friction to the problem we will instead get continuous solutions as we will see in the following sections. Figures \ref{fig:shock:pv} and \ref{fig:shock:Ts} shows a normal shock process in a $p\nu$- and $Ts$-diagram, respectively. Note that the flow passes the characteristic conditions as it is going through the shock, which means that the flow goes from supersonic to subsonic.

One-dimensional flow with heta addition

Flow-station relations

The aim is to derive relations for pressure ratio and temperature ratio as a function of Mach numbers. We will do that starting from the momentum equation.

$p_2-p_1={\rho }_1{u}_1^2-{\rho }_2{u}_2^2$

(Eq. 194)

Assuming calorically perfect gas

$\rho u^2=\rho a^2{M}^2=\rho \frac{\gamma p}{\rho }{M}^2=\gamma p{M}^2$

(Eq. 195)

which inserted in Eqn. \ref{eq:governing:mom} gives

$p_2-p_1=\gamma p_1{M}_1^2-\gamma p_2{M}_2^2$

(Eq. 196)

$p_2(1+\gamma M_2^2)=p_1(1+\gamma M_1^2)$

(Eq. 197)

and thus

$\frac{p_2}{p_1}=\frac{1+\gamma M_1^2}{1+\gamma M_2^2}$

(Eq. 198)

From the equation of state $p=\rho RT$, we get

$\frac{T_2}{T_1}=\frac{p_2}{{\rho }_{2}R}\frac{{\rho }_{1}R}{p_1}=\frac{p_2}{p_1}\frac{{\rho }_1}{{\rho }_2}$

(Eq. 199)

Using the continuity equation, we can get ${\rho }_1/{\rho }_2$

${\rho }_1{u}_1={\rho }_2{u}_2\Rightarrow \frac{{\rho }_1}{{\rho }_2}=\frac{u_2}{u_1}$

(Eq. 200)

Inserted in Eqn. \ref{eq:tr:a} gives

$\frac{T_2}{T_1}=\frac{p_2}{p_1}\frac{u_2}{u_1}$

(Eq. 201)

$\frac{u_2}{u_1}=\frac{M_2{a}_2}{M_1{a}_1}=\frac{M_2}{M_1}\frac{\sqrt{\gamma R{T}_2}}{\sqrt{\gamma R{T}_1}}=\frac{M_2}{M_1}\sqrt{\frac{T_2}{T_1}}$

(Eq. 202)

Eqn. \ref{eq:tr:c} in Eqn. \ref{eq:tr:b} gives

$\sqrt{\frac{T_2}{T_1}}=\frac{p_2}{p_1}\frac{M_2}{M_1}$

(Eq. 203)

With $p_2/p_1$ from Eqn. \ref{eq:governing:mom:b}, Eqn \ref{eq:tr:d} becomes

$\frac{T_2}{T_1}={\left(\frac{1+\gamma M_1^2}{1+\gamma M_2^2}\right)}^2{\left(\frac{M_2}{M_1}\right)}^2$

(Eq. 204)

Differential Relations

The equations presented in the previous section gives us the flow state after heat addition but since the heat addition, unlike the normal shock, is a continuous process, it is of interest to study the the heat addition from start to end. In order to do so we will now derive differential relations starting from the governing equations on differential form. We will start with converting the integral equation for conservation of mass for one-dimensional flows to differential form.

${\rho }_1{u}_1={\rho }_2{u}_2=const\Rightarrow d(\rho u)=0$

(Eq. 205)

$d(\rho u)=\rho du+ud\rho =0$

(Eq. 206)

Divide by $\rho u$ gives

$\frac{d\rho }{\rho }=\frac{du}{u}$

(Eq. 207)

The integral form of the conservation of momentum equation for one-dimensional flows is converted to differential form as follows.

$p_1+{\rho }_1{u}_1^2=p_2+{\rho }_2{u}_2^2=const\Rightarrow d(p+\rho u^2)=0$

(Eq. 208)

$dp+\rho udu+u\underset{=0}{\underbrace{d(\rho u)}}=0\Rightarrow dp=-\rho udu$

(Eq. 209)

with $\rho =\frac{p}{RT}$ and $u^2=M^2{a}^2=M^2\gamma RT$ in Eqn.~\ref{eq:governing:mom:diff:b}, we get

$dp=-\frac{p}{RT}{u}^2\frac{du}{u}=-\frac{p}{RT}{M}^2\gamma RT\frac{du}{u}\Rightarrow$

(Eq. 210)

$\frac{dp}{p}=-\gamma M^2\frac{du}{u}$

(Eq. 211)

which gives the relative change in pressure, $dp/p$, as a function of the relative change in flow velocity, $du/u$. The next equation to derive is an equation that describes the relative change in temperature, $dT/T$, as a function of the relative change in flow velocity, $du/u$. The starting point is the equation of state (the gas law).

$p=\rho RT\Rightarrow dp=R(\rho dT+Td\rho )\Rightarrow dT=\frac{1}{R\rho }dp-\frac{T}{\rho }d\rho$

(Eq. 212)

Divide by $T$

$\frac{dT}{T}=\frac{1}{\rho RT}dp-\frac{1}{\rho }d\rho$

(Eq. 213)

substitute $dp$ from Eqn.~\ref{eq:governing:mom:diff:c} and $d\rho$ from Eqn.~\ref{eq:governing:cont:diff:b} gives

$\frac{dT}{T}=(1-\gamma M^2)\frac{du}{u}$

(Eq. 214)

The entropy equation reads

$ds=C_v\frac{dp}{p}-C_p\frac{d\rho }{\rho }$

(Eq. 215)

which after substituting $dp$ from Eqn.~\ref{eq:governing:mom:diff:c} and $d\rho$ from Eqn.~\ref{eq:governing:cont:diff:b} becomes

$ds=C_v\gamma (1-M^2)\frac{du}{u}$

(Eq. 216)

From the definition of total temperature $T_o$ we get

$T_o=T+\frac{u^2}{2C_p}\Rightarrow d{T}_o=dT+\frac{1}{C_p}udu=dT+\frac{\gamma -1}{\gamma RT}T{u}^2\frac{du}{u}\Rightarrow$

(Eq. 217)

$d{T}_o=dT+(\gamma -1)M^{2}T\frac{du}{u}$

(Eq. 218)

Inserting $dT$ from Eqn~\ref{eq:governing:temp:diff:c} in Eqn~\ref{eq:governing:To:diff:a} we get

$d{T}_o=(1-\gamma M^2)T\frac{du}{u}+(\gamma -1)M^{2}T\frac{du}{u}$

(Eq. 219)

or

$d{T}_o=(1-M^2)T\frac{du}{u}$

(Eq. 220)

Dividing Eqn.~\ref{eq:governing:To:diff:b} by $T_o$ and using

$T_o=T(1+\frac{\gamma -1}{2}{M}^2)$

(Eq. 221)

we get

$\frac{d{T}_o}{T_o}=\frac{1-M^2}{1+\frac{\gamma -1}{2}{M}^2}\frac{du}{u}$

(Eq. 222)

Finally, we will derive a differential relation that describes the change in Mach number.

$M=\frac{u}{\sqrt{\gamma RT}}\Rightarrow dM=\frac{1}{\sqrt{\gamma R}}(T^{1/2}du+ud(T^{-1/2}))=\frac{du}{\sqrt{\gamma RT}}-\frac{u}{2\sqrt{\gamma R}}{T}^{-3/2}dT\Rightarrow$

(Eq. 223)

$dM=\frac{1}{\sqrt{\gamma RT}}\frac{du}{u}-\frac{1}{2}\frac{u}{\sqrt{\gamma RT}}\frac{dT}{T}=M\frac{du}{u}-\frac{M}{2}\frac{dT}{T}$

(Eq. 224)

Inserting $dT/T$ from Eqn.~\ref{eq:governing:temp:diff:c}, we get

$\frac{dM}{M}=\frac{1+\gamma M^2}{2}\frac{du}{u}$

(Eq. 225)

All the derived differential relations are expressed as functions of <math<du/u</math> but it would be more convenient to relate the changes in flow properties to the added heat or the change in total temperature, which can be related to the added heat through the energy equation.

$d{T}_o=\frac{\delta q}{C_p}$

(Eq. 226)

From Eqn.~\ref{eq:governing:To:diff:c}, we get

$\frac{du}{u}=\left(\frac{1+\frac{\gamma -1}{2}{M}^2}{1-M^2}\right)\frac{d{T}_o}{T_o}$

(Eq. 227)

Now, we can substitute $du/u$ in all the above relations using Eqn.~\ref{eq:governing:du:diff:final}, we get the following relations

$\frac{d\rho }{\rho }=-\left(\frac{1+\frac{\gamma -1}{2}{M}^2}{1-M^2}\right)\frac{d{T}_o}{T_o}$

(Eq. 228)

$\frac{dp}{p}=\gamma M^2\left(\frac{1+\frac{\gamma -1}{2}{M}^2}{1-M^2}\right)\frac{d{T}_o}{T_o}$

(Eq. 229)

$\frac{dT}{T}=(1-\gamma M^2)\left(\frac{1+\frac{\gamma -1}{2}{M}^2}{1-M^2}\right)\frac{d{T}_o}{T_o}$

(Eq. 230)

$\frac{dT}{T}=\left(\frac{1+\gamma M^2}{2}\right)\left(\frac{1+\frac{\gamma -1}{2}{M}^2}{1-M^2}\right)\frac{d{T}_o}{T_o}$

(Eq. 231)

$\frac{dT}{T}=C_v\gamma (1+\frac{\gamma -1}{2}{M}^2)\frac{d{T}_o}{T_o}$

(Eq. 232)

Heat Addition Process

With the differential relations in place, we can now study the continuous change in flow quantities from the initial flow state to the flow state after the heat addition process by dividing the total amount of heat added to the flow, $q$, into small portions, $\delta q$, and calculate the change in flow properties for each of these heat additions, see Figure~\ref{fig:dq}.

Let's first examine the temperature change by rewriting Eqn.~\ref{eq:governing:dT:diff:final} as

$dT=\frac{1-\gamma M^2}{1-M^2}d{T}_o\iff \frac{dT}{d{T}_o}=\frac{1-\gamma M^2}{1-M^2}$

(Eq. 233)

which is equivalent to

$\frac{dh}{\delta q}=\frac{1-\gamma M^2}{1-M^2}$

(Eq. 234)

Form Eqn.~\ref{eq:governing:dT:diff:mod:a} we can make the following observation

$\frac{dT}{d{T}_o}=0\Rightarrow \gamma M^2=1\Rightarrow M=\sqrt{1/\gamma }$

(Eq. 235)

which means that the maximum temperature will be reached when the Mach number is $\sqrt{1/\gamma }$. Since $\gamma$ is a number greater than one for all gases, this implies that the maximum temperature can only be reached if the flow is subsonic. For air, this the maximum temperature will be reached at $M=0.845$.

If we evaluate Eqn.~\ref{eq:governing:dT:diff:mod:a} for sonic flow ($M=1$), we see that the derivative becomes infinite.

$|M|\to 1.0\Rightarrow \frac{dT}{d{T}_o}\to \pm \mathrm{\infty }$

(Eq. 236)

Now, by specifying an initial subsonic flow state and dividing the heat addition corresponding to choked flow, $q^{\ast }$, into small portions $\delta q$, one can perform integration as indicated in Figure~\ref{fig:dq}. The result is presented in the in Figure~\ref{fig:TS:closeup}. The subsonic process corresponds to the upper line. As heat is added the Mach number is increased and at $M={\gamma }^{-1/2}$ the maximum temperature is reached. Adding more heat will reduce the temperature and increase the Mach number until sonic conditions are reached ($M=1.0$). As can be seen in Figure~\ref{fig:TS:closeup}, the lean of the subsonic branch of the Rayleigh line is lower than the isobars (gray lines), which means the increasing heat will reduce pressure. The lower part of the blue line in Figure~\ref{fig:TS:closeup} is the supersonic branch of the Rayleigh line, which is obtained in the same way starting from a supersonic flow condition. A flow state resulting in the same sonic conditions as for the subsonic case is calculated and used as a starting state. The corresponding $q^\ast$ is calculated and the same calculation of consecutive flow states in a step-wise manner is performed. As can be seen in Figure~\ref{fig:TS:closeup}, the lean of the supersonic part of the Rayleigh curve is steeper than the isobars (gray lines), which means that pressure increases as heat is added to the flow. As we saw from Eqn.~\ref{eq:governing:dT:diff:mod:b}, $dT/d{T}_o$ becomes infinite when the flow approaches the sonic the sonic state. After the sonic state is reached, further heat addition is impossible without changing the upstream flow conditions. This will be made clearer in the next section.

Using the differential relations above, we can get a good picture of the development of flow variables as heat is continuously added to the flow (see Figure~\ref{fig:rayleigh:trends}).

Rayleigh Line

The continuity equation for steady-state, one-dimensional flow reads

${\rho }_1{u}_1={\rho }_2{u}_2=C$

(Eq. 237)

where $C$ is the massflow per square meter (massflow divided by area). Inserted in the momentum equation we get

$p_1+\frac{C^2}{{\rho }_1}=p_2+\frac{C^2}{{\rho }_2}\iff p_1+{\nu }_1{C}^2=p_2+{\nu }_2{C}^2\Rightarrow \frac{p_2-p_1}{{\nu }_2-{\nu }_1}=-C^2$

(Eq. 238)

Eqn.~\ref{eq:governing:mom:b} tells us that any solution to the governing flow equations must lie along a line (a so-called Rayleigh line) in a $p\nu$-diagram. In Figure~\ref{fig:PV}, 1 corresponds to the flow state before heat addition and states 2 and 3 corresponds to the flow state after heat is added. If the flow in state 1 is subsonic, adding heat will change the flow state following the Rayleigh line to the right, i.e. towards flow state 2. If the initial flow state instead is supersonic, heat addition will move the flow state towards state 3.

Now we know in which direction we will move along the Rayleigh curve when heat is added but in order to find the flow state after heat addition we need to add the energy equation to the problem. If we draw a curve corresponding to the energy equation including the heat addition in the same $p\nu$-diagram, the intersection of this curve and the Rayleigh line corresponds to the downstream flow state (the flow state that fulfils the continuity, momentum, and energy equations). To be able to do this we will rewrite the energy equation such that it can be represented by a line in the $p\nu$-diagram.

The energy equation for one-dimensional flow with heat addition reads

$h_1+\frac{1}{2}{u}_1^2+q=h_2+\frac{1}{2}{u}_2^2$

(Eq. 239)

Inserting the constant $C$ from above (the massflow per $m^2$) and and and $h=C_{p}T$, we get

$\frac{\gamma R}{\gamma -1}{T}_1+\frac{1}{2}{C}^2{\nu }_1^2+q=\frac{\gamma R}{\gamma -1}+\frac{1}{2}{C}^2{\nu }_2^2$

(Eq. 240)

which may be rewritten as

$\frac{p_2}{p_1}=(\frac{{\nu }_2}{{\nu }_1}-\frac{\gamma +1}{\gamma -1}-\frac{2q}{R{T}_1}){(1-\frac{\gamma +1}{\gamma -1}\frac{{\nu }_2}{{\nu }_1})}^{-1}$

(Eq. 241)

As you can see in the examples above (Figures~\ref{fig:TSPV:b} and \ref{fig:TSPV:d}), sonic conditions are reached when the Rayleigh line is tangent to the curve representing the energy equation in the $p\nu$-diagram. Adding more heat would move the energy equation line upwards and thus there can not be any solution after reaching this state unless the upstream conditions are changed such that the energy line intersects the Rayleigh line after further heat addition. Let's have a second look at the equations and see if it is possible to verify that the case where the Rayleigh line is a tangent to the energy-equation curve is in fact the sonic state.

Starting from Eqn.~\ref{eq:governing:energy:b}, it is easy to see that for any point along the energy equation curve the flow state may be expressed as a function of the initial flow state and the added heat $q$ as

$\frac{\gamma }{\gamma -1}p\nu +\frac{1}{2}{C}^2{\nu }^2=\frac{\gamma }{\gamma -1}{p}_1{\nu }_1+\frac{1}{2}{C}^2{\nu }_1^2+q=D$

(Eq. 242)

where $D$ is a constant.

Now, let's different the Eqn.\ref{eq:governing:energy:d} with respect to $\nu$

$\frac{\gamma }{\gamma -1}(\nu \frac{dp}{d\nu }+p)+C^2\nu =0\Rightarrow \frac{dp}{d\nu }=-\frac{\gamma }{\gamma -1}{C}^2-\frac{p}{\nu }$

(Eq. 243)

The Rayleigh line is a tangent to the energy equation curve when $dp/d\nu =-C^2$ and thus

$\frac{C^2}{\gamma }=\frac{p}{\nu }$

(Eq. 244)

By definition $C=\rho u$ and $\nu =1/\rho$, which inserted in Eqn.~\ref{eq:governing:energy:f} gives

$u=\sqrt{\frac{\gamma p}{\rho }}=a$

(Eq. 245)

Thermal Choking

When the heat addition reaches $q^\ast$ the flow becomes sonic and the flow is said to thermally choked. Thermal choking is illustrated in Figure~\ref{fig:TSPV:d}, where the curve representing the energy equation (the blue line in the $p\nu$-diagram) is tangent to the Rayleigh line and if more heat is added the blue line will move to the right of the Rayleigh line and thus there are no solutions for $q>q^{\ast }$. So what happens if more heat is added to the flow after thermal choking is reached. The answer is different if the flow is subsonic or supersonic. For a subsonic flow, the upstream flow will be adjusted such that the slope of the Rayleigh line changes and the energy equation curve becomes tangent to the Rayleigh line. This means that the massflow per unit area ($C$) is reduced and $q^{\ast }$ is increased such that $q^{\ast }$ equals the heat added to the flow. Note that the upstream total conditions will not be changed in this process (see Figure~\ref{fig:thermal:choking:sub}).

$M_{1^{\prime }}=f(q^{\ast })$ $T_{1^{\prime }}=f(T_o,M_{1^{\prime }})$ $p_{1^{\prime }}=f(p_o,M_{1^{\prime }})$ ${\rho }_{1^{\prime }}=f(p_{1^{\prime }},T_{1^{\prime }})$ $a_{1^{\prime }}=f(T_{1^{\prime }})$ $u_{1^{\prime }}=M_{1^{\prime }}{a}_{1^{\prime }}$

In a choked supersonic flow, there is no possibility for pressure waves to travel upstream in the flow and thus the upstream flow conditions can not be changed as in the subsonic case. Moreover, since a normal shock is an adiabatic process (a jump between two points on the same Rayleigh line), the total temperature is not changed over a chock. From before we have

$\frac{T_o}{T_o^{\ast }}=\frac{(\gamma +1)M_1^2}{(1+\gamma M_1^2{)}^2}(2+(\gamma -1)M_1^2)$

(Eq. 246)

Inserting the normal shock relation

$M_2^2=\frac{2+(\gamma -1)M_1^2}{2\gamma M_1^2-(\gamma -1)}$

(Eq. 247)

one can show that

$\frac{T_o}{T_o^{\ast }}=\frac{(\gamma +1)M_2^2}{(1+\gamma M_2^2{)}^2}(2+(\gamma -1)M_2^2)$

(Eq. 248)

and thus $T_o^{\ast }$ is not changed by the normal shock and consequently $q^{\ast }$ is unchanged if there is a normal shock between station 1 and 2. So, it is not possible to change the upstream static flow conditions and a normal shock will not make it possible to add more heat. The only possible solution is a normal shock upstream of station 1 and thus subsonic flow through the heat addition process.

One-dimensional flow with friction

Flow-station data

The starting point is the governing equations for one-dimensional steady-state flow

Continuity

$${\rho }_1{u}_1={\rho }_2{u}_2$$

Momentum

${\rho }_1{u}_1^2+p_1-{\overline{\tau }}_w\frac{bL}{A}={\rho }_2{u}_2^2+p_2$

(Eq. 249)

where ${\overline{\tau }}_w$ is the average wall-shear stress

${\overline{\tau }}_w=\frac{1}{L}\underset{0}{\overset{L}{\int }}{\tau }_{w}dx$

(Eq. 250)

$b$ is the tube perimeter, and $L$ is the tube length. For circular cross sections

$\frac{bL}{A}=\{A=\frac{\pi D^2}{4},b=\pi D\}=\frac{4L}{D}$

(Eq. 251)

and thus

${\rho }_1{u}_1^2+p_1-\frac{4}{D}\underset{0}{\overset{L}{\int }}{\tau }_{w}dx={\rho }_2{u}_2^2+p_2$

(Eq. 252)

Energy

$h_1+\frac{1}{2}{u}_1^2=h_2+\frac{1}{2}{u}_2^2$

(Eq. 253)

Differential Form

In order to remove the integral term in the momentum equation, the governing equations are written in differential form

Continuity

${\rho }_1{u}_1={\rho }_2{u}_2=const\Rightarrow$

(Eq. 254)

$\frac{d}{dx}(\rho u)=0$

(Eq. 255)

Momentum

$({\rho }_2{u}_2^2+p_2-{\rho }_1{u}_1^2+p_1)=-\frac{4}{D}\underset{0}{\overset{L}{\int }}{\tau }_{w}dx\Rightarrow$

(Eq. 256)

$\frac{d}{dx}(\rho u^2+p)=-\frac{4}{D}{\tau }_w$

(Eq. 257)

$\frac{d}{dx}(\rho u^2+p)=\rho u\frac{du}{dx}+u\frac{d}{dx}(\rho u)+\frac{dp}{dx}=\{\frac{d}{dx}(\rho u)=0\}=\rho u\frac{du}{dx}+\frac{dp}{dx}$

(Eq. 258)

$\rho u\frac{du}{dx}+\frac{dp}{dx}=-\frac{4}{D}{\tau }_w$

(Eq. 259)

The wall shear stress is often approximated using a shear-stress factor, $f$, according to

${\tau }_w=f\frac{1}{2}\rho u^2$

(Eq. 260)

and thus

$\rho u\frac{du}{dx}+\frac{dp}{dx}=-\frac{2}{D}f\rho u^2$

(Eq. 261)

Energy

$h_1+\frac{1}{2}{u}_1^2=h_2+\frac{1}{2}{u}_2^2=const$

(Eq. 262)

$h_{o_1}=h_{o_2}=const$

(Eq. 263)

$\frac{d}{dx}{h}_o=0$

(Eq. 264)

Summary

continuity:

$\frac{d}{dx}(\rho u)=0$

(Eq. 265)

momentum:

$\rho u\frac{du}{dx}+\frac{dp}{dx}=-\frac{2}{D}f\rho u^2$

(Eq. 266)

energy:

$\frac{d}{dx}{h}_o=0$

(Eq. 267)

From chapter 3.9 we have the following expression for the momentum equation for one-dimensional flow with friction (equation (3.95))

$dp+\rho udu=-\frac{1}{2}\rho u^2\frac{4fdx}{D}$

(Eq. 268)

For cases dealing with calorically perfect gas, (3.95) can be recast completely in terms of Mach number using the following relations

• speed of sound: $a^2=\gamma p/\rho$
• the definition of Mach number: $M^2=u^2/a^2$
• the ideal gas law for thermally perfect gas: $p=\rho RT$
• the continuity equation: $\rho u=const$
• the energy equation: $C_{p}T+u^2/2=const$

Continuity equation

We start with the continuity equation which for one-dimensional steady flows reads

$\rho u=const$

(Eq. 269)

Differentiating (\ref{eqn:cont:a}) gives

$d(\rho u)=0.\iff \rho du+ud\rho =0.$

(Eq. 270)

If $u\ne 0.$ we can divide by $\rho u$ which gives us

$\frac{du}{u}+\frac{d\rho }{\rho }=0.$

(Eq. 271)

Now, if we divide and multiply the first term in (\ref{eqn:cont:b}) by $2u$ and use the chain rule for derivatives we get

$\frac{d(u^2)}{2u^2}+\frac{d\rho }{\rho }=0.$

(Eq. 272)

Energy equation

For an adiabatic one-dimensional flow we have that

$C_{p}T+\frac{u^2}{2}=const$

(Eq. 273)

If we differentiate (\ref{eqn:ttot:a}) we get

$C_{p}dT+\frac{1}{2}d(u^2)=0.$

(Eq. 274)

We replace $C_p$ with $\gamma R/(\gamma -1)$ and multiply and divide the first term with $T$ which gives us

$\frac{\gamma RT}{(\gamma -1)}\frac{dT}{T}+\frac{1}{2}d(u^2)=0.$

(Eq. 275)

Now, divide by $\gamma RT/(\gamma -1)$ and multiply and divide the second term by $u^2$ gives

$\frac{dT}{T}+\frac{(\gamma -1)}{2}{M}^2\frac{d(u^2)}{u^2}=0.$

(Eq. 276)

We want to remove the $dT/T$-term in (\ref{eqn:ttot}). From the definition of Mach number we have that

$a^2{M}^2=u^2$

(Eq. 277)

which we can rewrite using the expression for speed of sound $(a^2=\gamma RT)$ according to

$\gamma RT{M}^2=u^2$

(Eq. 278)

Differentiating (\ref{eqn:Mach:b}) gives us

$\gamma R{M}^{2}dT+\gamma RTd(M^2)=d(u^2)$

(Eq. 279)

Now, if we divide (\ref{eqn:Mach:c}) by $\gamma RT{M}^2$ and use $a^2=\gamma RT$ and $a^2{M}^2=u^2$ we get

$\frac{dT}{T}+\frac{d(M^2)}{M^2}=\frac{d(u^2)}{u^2}$

(Eq. 280)

Equation (\ref{eqn:Mach}) may now be used to replace the $dT/T$-term in equation (\ref{eqn:ttot})

$-\frac{d(M^2)}{M^2}+\frac{d(u^2)}{u^2}+\frac{(\gamma -1)}{2}{M}^2\frac{d(u^2)}{u^2}=0.$

(Eq. 281)

which can be rewritten according to

$\frac{d(u^2)}{u^2}={[1+\frac{(\gamma -1)}{2}{M}^2]}^{-1}\frac{d(M^2)}{M^2}$

(Eq. 282)

Using the chain rule for derivatives, the last term may be rewritten according to

$\frac{d(M^2)}{M^2}=2M\frac{dM}{M^2}=2\frac{dM}{M}$

(Eq. 283)

which gives

$\frac{d(u^2)}{u^2}=2{[1+\frac{(\gamma -1)}{2}{M}^2]}^{-1}\frac{dM}{M}$

(Eq. 284)

The ideal gas law

For a perfect gas the ideal gas law reads

$p=\rho RT$

(Eq. 285)

Differentiating (\ref{eqn:gaslaw:a}) gives:

$dp=\rho RdT+RTd\rho$

(Eq. 286)

If $p\ne 0.$, we can divide (\ref{eqn:gaslaw:b}) by $p$ which gives

$\frac{dp}{p}=\frac{dT}{T}+\frac{d\rho }{\rho }$

(Eq. 287)

which can be rearranged according to

$[\frac{dp}{p}-\frac{d\rho }{\rho }]=\frac{dT}{T}$

(Eq. 288)

Now, inserting $dT/T$ from equation (\ref{eqn:ttot}) gives

$[\frac{dp}{p}-\frac{d\rho }{\rho }]+\frac{(\gamma -1)}{2}{M}^2\frac{d(u^2)}{u^2}=0.$

(Eq. 289)

The $d\rho /\rho$-term can be replaced using equation (\ref{eqn:cont})

$\frac{dp}{p}+\frac{d(u^2)}{2u^2}+\frac{(\gamma -1)}{2}{M}^2\frac{d(u^2)}{u^2}=0.$

(Eq. 290)

Collect terms and rewrite gives

$\frac{dp}{p}+\left[\frac{1+(\gamma -1)M^2}{2}\right]\frac{d(u^2)}{u^2}=0.$

(Eq. 291)

Momentum equation

By combining the above derived relations and the momentum equation on the form given by (3.95), we can get an expression where the friction force is a function of Mach number only

For convenience equation (3.95) is written again here

$dp+\rho udu=-\frac{1}{2}\rho u^2\frac{4fdx}{D}$

(Eq. 292)

if $u\ne 0.$, we can divide by $0.5\rho u^2$ which gives

$2\frac{dp}{\rho u^2}+2\frac{\rho udu}{\rho u^2}=-\frac{4fdx}{D}$

(Eq. 293)

using $M^2=u^2/a^2$, $a^2=\gamma p/\rho$ and the chain rule in (\ref{eqn:mom:a}) gives

$\frac{2}{\gamma M^2}\frac{dp}{p}+\frac{d(u^2)}{u^2}=-\frac{4fdx}{D}$

(Eq. 294)

From equation (\ref{eqn:gaslaw}) we can get a relation that expresses the pressure derivative term, $dp/p$, in terms of Mach number and $d(u^2)/u^2$. Inserting this in (\ref{eqn:mom:b}) gives

$\frac{2}{\gamma M^2}\{-\left[\frac{1+(\gamma -1)M^2}{2}\right]\frac{d(u^2)}{u^2}\}+\frac{d(u^2)}{u^2}=-\frac{4fdx}{D}$

(Eq. 295)

collecting terms and rearranging gives

$\frac{M^2-1}{\gamma M^2}\frac{d(u^2)}{u^2}=\frac{4fdx}{D}$

(Eq. 296)

if we now use equation (\ref{eqn:ttot:Mach}) to get rid of the $d(u^2)/u^2$-term we end up with the following expression

$\frac{4fdx}{D}=\frac{2}{\gamma M^2}(1-M^2){[1+\frac{(\gamma -1)}{2}{M}^2]}^{-1}\frac{dM}{M}$

(Eq. 297)

Differential Relations

In analogy with the heat addition process discussed in the previous section, one-dimensional flow with heat addition is a continuous process. We will derive the differential relations for one-dimensional flow with friction, which will lead to trends for supersonic and supersonic flow with friction.

The continuity equation gives

$d(\rho u)=ud\rho +\rho du\Rightarrow$

(Eq. 298)

$\frac{d\rho }{\rho }=-\frac{du}{u}$

(Eq. 299)

The addition of friction does not affect total temperature and thus the total temperature is constant

$T_o=T+\frac{u^2}{2C_p}=const$

(Eq. 300)

differentiating gives

$d{T}_o=dT+\frac{1}{Cp}udu=0$

(Eq. 301)

with $u=M\sqrt{\gamma RT}$, we get

$\frac{dT}{T}=-(\gamma -1)M^2\frac{du}{u}$

(Eq. 302)

A differential relation for pressure can be obtained from the ideal gas relation

$p=\rho RT\Rightarrow dp=R(Td\rho +\rho dT)\Rightarrow$

(Eq. 303)

$\frac{dp}{p}=-(1+(\gamma -1)M^2)\frac{du}{u}$

(Eq. 304)

The entropy increase can be obtained from

$ds=C_v\frac{dp}{p}-C_p\frac{d\rho }{\rho }$

(Eq. 305)

and thus

$ds=-R(1-M^2)\frac{du}{u}$

(Eq. 306)

Finally, a relation describing the change in Mach number can be obtained from

$M=\frac{u}{\sqrt{\gamma RT}}\Rightarrow dM=M\frac{du}{u}-\frac{M}{2}\frac{dT}{T}$

(Eq. 307)

which can be rewritten as

$\frac{dM}{M}=(1+(\gamma -1)M^2)\frac{du}{u}$

(Eq. 308)

Eqns.~\ref{eq:fanno:drho} - \ref{eq:fanno:dM} are expressed as functions of $du$ and in order to get a direct relation to the addition of friction caused by the increase in pipe length $dx$, the equations are rewritten so that all variable changes are functions of the entropy increase $ds$.

$\frac{d\rho }{\rho }=-\frac{1}{R(1-M^2)}ds$

(Eq. 309)

$\frac{dT}{T}=-(\gamma -1)M^2\frac{1}{R(1-M^2)}ds$

(Eq. 310)

$\frac{dp}{p}=-(1+(\gamma -1)M^2)\frac{1}{R(1-M^2)}ds$

(Eq. 311)

$\frac{dM}{M}=(1+\frac{(\gamma -1)}{2}{M}^2)\frac{1}{R(1-M^2)}ds$

(Eq. 312)

$\frac{du}{u}=\frac{1}{R(1-M^2)}ds$

(Eq. 313)

A relation for the change in total pressure can be obtained from

$ds=C_p\frac{d{T}_o}{T_o}-R\frac{d{p}_o}{p_o}$

(Eq. 314)

Since total temperature is constant the relation above gives

$\frac{d{p}_o}{p_o}=-\frac{ds}{R}$

(Eq. 315)

Using the differential relations above, we can get a good picture of the development of flow variables as friction is continuously added to the flow (see Figure~\ref{fig:fanno:trends}).

Friction Choking

Figure~\ref{fig:friction:Ts} shows the Fanno flow process in a $Ts$-diagram. The dashed line represents the sonic temperature, which means that the flow states along the process line above the dashed line are subsonic flow states and the part of the line below the dashed line represents supersonic flow states. In both subsonic and supersonic flow addition of friction leads to a change in temperature in the direction towards the sonic temperature, i.e. the flow approaches sonic conditions ($M=1$). When the length of the pipe through which the fluid flows is equal to the length at which the flow is sonic, the flow is choked (friction choking) and further pipe length cannot be added without a change in the flow conditions. For an initially subsonic flow, a pipe longer than $L^{\ast }$, the change in flow conditions is analogous to the what happens for addition of heat to a subsonic flow that has reached sonic state discussed in the previous section. The inlet conditions will change such that the massflow is reduced without changing the inlet total conditions such as the pipe length is equal to $L^{\ast }$ for the new inlet conditions.

$M_{1^{\prime }}=f(L^{\ast })$ $T_{1^{\prime }}=f(T_o,M_{1^{\prime }})$ $p_{1^{\prime }}=f(p_o,M_{1^{\prime }})$ ${\rho }_{1^{\prime }}=f(p_{1^{\prime }},T_{1^{\prime }})$ $a_{1^{\prime }}=f(T_{1^{\prime }})$ $u_{1^{\prime }}=M_{1^{\prime }}{a}_{1^{\prime }}$

For a choked supersonic flow, addition of more friction (increasing the length of the pipe such that $L>L^{\ast }$) may lead to the generation of a shock inside the pipe. In contrast to the one-dimensional flow with heat addition where a shock does not change $q^{\ast }$, $L^{\ast }$ is increased over a shock. The internal shock will be generated in an axial location such that $L^{\ast }$ downstream of the shock equals the remaining pipe length at the shock location (see Figure~\ref{fig:friction:choking:sup}). As more length is added to the pipe, the shock will move further and further upstream in the pipe until it stands at the pipe entrance. If the pipe is longer than $L^{\ast }$ after o shock standing at the inlet, the shock will move to the upstream system and the pipe flow will be subsonic and the massflow will be adjusted such that $L=L^{\ast }$ according to the process described for subsonic choking above.

From prvevious derivations, we know that $L^{\ast }$ is a function of mach number according to

$\frac{4\overline{f}{L}^{\ast }}{D}=\frac{1-M^2}{\gamma M^2}+\left(\frac{\gamma +1}{2\gamma }\right)\ln \left(\frac{(\gamma +1)M^2}{2+(\gamma -1)M^2}\right)$

(Eq. 316)

by dividing both the numerator and denominator in the fractions by $M^2$ it is easy to see that the choking length (Figure~\ref{fig:friction:factor}) approaches a finite length for great Mach numbers and thus the upper limit for the choking length $L_1^{\ast }$ is given by

${\frac{4\overline{f}{L}_1^{\ast }}{D}(M_1)|}_{M_1\to \mathrm{\infty }}=-\frac{1}{\gamma }+\left(\frac{\gamma +1}{2\gamma }\right)\ln \left(\frac{\gamma +1}{\gamma -1}\right)$

(Eq. 317)

From the normal shock relations we know that the downstream Mach number approaches the finite value $\sqrt{(\gamma -1)/2\gamma }$ large Mach numbers and thus the choking length downstream the shock is limited to

${\frac{4\overline{f}{L}_2^{\ast }}{D}(M_2)|}_{M_1\to \mathrm{\infty }}=\left(\frac{\gamma +1}{\gamma (\gamma -1)}\right)+\left(\frac{\gamma +1}{2\gamma }\right)\ln \left(\frac{(\gamma +1)(\gamma -1)}{4\gamma +(\gamma -1{)}^2}\right)$

(Eq. 318)

From the relations above we get

${(\frac{4\overline{f}{L}_2^{\ast }}{D}(M_2)-\frac{4\overline{f}{L}_1^{\ast }}{D}(M_1))|}_{M_1\to \mathrm{\infty }}=\left(\frac{2}{\gamma -1}\right)+\left(\frac{\gamma +1}{2\gamma }\right)\ln \left[\left(\frac{(\gamma +1)(\gamma -1)}{4\gamma +(\gamma -1{)}^2}\right)\left(\frac{\gamma -1}{\gamma +1}\right)\right]$

(Eq. 319)

Figure~\ref{fig:friction:factor:shock} shows the development of choking length $L_1^{\ast }$ in a supersonic flow as a function of Mach number in relation to the corresponding choking length $L_2^{\ast }$ downstream of a normal shock generated at the same Mach number. As can be seen from the figure, a normal shock will always increase the choking length.

Oblique shocks and expansion waves

Oblique shocks

The Shock Polar

The shock polar is a graphical representation of all possible flow deflection angles for a given Mach number. The shock polar is generated by plotting the normalized axial velocity component downstream of an oblique shock versus the normalized vertical velocity component. In essence the ratio of vertical and axial velocity components gives the deflection angle. Figure~\ref{fig:shock:polar:a} shows a set of shock polars generated for different upstream Mach numbers (indicated in the figure). The vertical and axial velocity components are normalized by the characteristic speed of sound (the speed of sound at sonic conditions). Each of the shock polars have two solutions where the vertical velocity component is zero, i.e. zero-deflection solutions. The zero-deflection solution furthest to the right represents the Mach wave solution and the solution to the left represents a normal shock. This is easy to realize as these are the only possible solutions that would result in zero flow deflection and the Mach wave solution is located where $V_x/a^{\ast }$ is greater than one, which means that the flow is supersonic on the downstream side, whereas the normal shock results in a subsonic flow on the downstream side ($V_x/a^{\ast }<1.0$). Figures~\ref{fig:shock:polar:b}-{fig:shock:polar:f} shows only the shock polar corresponding to an upstream mach number of 2.5. In Figure~\ref{fig:shock:polar:b}, a unit half-circle is added to indicate which parts of the shock polar that represents supersonic solutions and which parts that represent supersonic solutions. Everything that falls inside of the unit circle represents subsonic solutions. i.e. the downstream Mach number is subsonic and all solutions that are outside of the circle represents solutions for which the flow is supersonic. Figure~\ref{fig:shock:polar:c} shows how the maximum possible flow deflection relates to the shock polar for a specific upstream Mach number. The line in the figure represents the flow deflection and when the line is tangent to the shock polar, the maximum flow deflection is reached. Further increase of the flow deflection angle leads to that the flow deflection line is outside of the shock polar and thus there are no possible solutions for angles greater than ${\theta }_{max}$. For the maximum deflection angle there is only one possible solution since the flow deflection line is a tangent to the shock polar. For all angles smaller than ${\theta }_{max}$, there are, however, two possible solutions (see Figure~\ref{fig:shock:polar:d}). In most cases there is one supersonic solution (the weak solution) and one subsonic solution (the strong solution) as indicated in Figure~\ref{fig:shock:polar:d}). However, for angles close to ${\theta }_{max}$ both solutions may fall inside of the unit circle and thus both solutions will be subsonic. This is in line with what we saw for the $\theta$-$\beta$-Mach relation earlier (see section~\ref{sec:theta:beta:Mach}). Finally, we will have a look at how the shock angle ($\beta$) is related to the shock polar. Figure~\ref{fig:shock:polar:e} shows how the shock angle ($\beta$) is found for the weak solution of a given upstream mach number and a given flow deflection ($\theta$). Draw a line starting at the Mach wave solution going through the the downstream solution (in this case the weak solution). Now, make another line starting at the origin that is perpendicular to the first line. The angle of the second line is the shock angle ($\beta$). In analogy, figure~\ref{fig:shock:polar:f} shows how to obtain the shock angle for the strong solution.

Shock Reflection at a Solid Boundary

What happens when an oblique shock reaches a solid wall? To sort this out we will analyze the schematic flow situation illustrated in Figure~\ref{fig:regular:reflection}. The figure shows a supersonic flow through a channel where there is a sudden bend of the lower wall leading to the generation of an oblique shock in order to deflect the flow such that it follows the wall downstream of the corner. The shock angle ${\beta }_a$ is a function of the upstream Mach number $M_1$ and the flow deflection $\theta$. A bit further downstream the shock will reach the upper wall of the channel. The question now is what will happen at the point reaches the upper wall. As indicated in the figure the shock will deflect but in what way will it deflect. Will it be a specular reflection, i.e. will the angle of the shock be the same but in the other direction? To answer this question let's sort out why the shock will reflect. After the first shock, the flow will be deflected the angle $\theta$, which means that, at the upper wall, the flow must be deflected again such that it follows the direction of the upper wall. Hence, the deflection angle will again be $\theta$ but in the opposite direction. So, the deflection angle is the same, does that mean that the shock angle will be the same? The answer is no, but why? Passing the first shock, the Mach number is reduced and thus the $\theta$-$\beta$-Mach relation will give us another shock angle ${\beta }_2$ for the second shock. Hence, the shock is not reflected specularly since ${\beta }_1\ne {\beta }_2$.

The situation discussed in the previous paragraph assumed that the flow deflection at the upper wall was less than the maximum flow deflection possible for the Mach number ahead of the reflection (downstream of the first shock). If, however, the flow deflection that must take place exceeds the maximum possible flow deflection angle, it will not be possible to generate an oblique shock that fulfills the requirements. Instead, a normal shock will be generated at the upper wall (a so-called Mach reflection) that will be gradually converted into an oblique shock. This situation is depicted in Figure~\ref{fig:Mach:reflection}. The slip line indicated in the figure is a consequence of the fact that the flow the goes through the shock system experiences different entropy increases depending on if the flow passes a single stronger shock that generate higher losses or a set of weaker oblique shocks that will generate less losses.

Expansion waves

Prandtl-Meyer Expansion Waves

A single Mach wave has a insignificant effect on the flow passing it but an expansion region constitutes an infinite number of Mach waves and the integrated effect is significant. The net turning of the flow by a single Mach wave is depicted schematically in Fig. \ref{fig:machwave}. It can be shown geometrically that

$d\theta =\sqrt{M^2-1}\frac{dV}{V}$

(Eq. 320)

Since Eqn. \ref{eq:mach:turning} is derived from the flow turning geometry with the assumption that the net flow tuning is small, it is valid for all gas models.

To get the integrated effect of all Mach waves in the expansion region, we integrate Eqn. \ref{eq:mach:turning} over the expansion region

$\underset{{\theta }_1}{\overset{{\theta }_2}{\int }}d\theta =\underset{M_1}{\overset{M_2}{\int }}\sqrt{M^2-1}\frac{dV}{V}$

(Eq. 321)

To be able to do the integration, we need to rewrite it

$V=Ma\Rightarrow \ln V=\ln M+\ln a$

(Eq. 322)

Differentiate to get

$\frac{dV}{V}=\frac{dM}{M}+\frac{da}{a}$

(Eq. 323)

Each Mach wave is isentropic and thus the expansion is an isentropic process, which means that we can use the adiabatic energy equation

$\frac{T_o}{T}=1+\frac{\gamma -1}{2}{M}^2$

(Eq. 324)

For a calorically perfect gas $a=\sqrt{\gamma RT}$ and $a_o=\sqrt{\gamma R{T}_o}$ and thus

$\frac{T_o}{T}={\left(\frac{a_o}{a}\right)}^2\Rightarrow$

(Eq. 325)

${\left(\frac{a_o}{a}\right)}^2=1+\frac{\gamma -1}{2}{M}^2$

(Eq. 326)

Solve for $a$ gives

$a=a_o{(1+\frac{\gamma -1}{2}{M}^2)}^{-1/2}$

(Eq. 327)

Differentiate Eqn \ref{eq:adiatbatic:energy:b} to get

$da=a_o(-\frac{1}{2})\left(\frac{\gamma -1}{2}\right)2M{(1+\frac{\gamma -1}{2}{M}^2)}^{-3/2}dM$

(Eq. 328)

Eqn. \ref{eq:adiatbatic:energy:b} in Eqn. \ref{eq:adiatbatic:energy:c} gives

$\frac{da}{a}=-\left(\frac{\gamma -1}{2}\right){(1+\frac{\gamma -1}{2}{M}^2)}^{-1}MdM$

(Eq. 329)

From Eqn. \ref{eq:mach:turning:c}, we have

$\frac{dV}{V}=\frac{dM}{M}+\frac{da}{a}$

(Eq. 330)

With $da/a$ from Eqn. \ref{eq:adiatbatic:energy:d}, we get

$\frac{dV}{V}=\frac{dM}{M}-\left(\frac{\gamma -1}{2}\right){(1+\frac{\gamma -1}{2}{M}^2)}^{-1}MdM$

(Eq. 331)

$\frac{dV}{V}=dM\left[\frac{(1+\frac{\gamma -1}{2}{M}^2)-\left(\frac{\gamma -1}{2}\right)M^2}{(1+\frac{\gamma -1}{2}{M}^2)M}\right]={(1+\frac{\gamma -1}{2}{M}^2)}^{-1}\frac{dM}{M}$

(Eq. 332)

Now, insert $dV/V$ in Eqn. \ref{eq:mach:turning:b} to get

$\underset{{\theta }_1}{\overset{{\theta }_2}{\int }}d\theta =\underset{M_1}{\overset{M_2}{\int }}\sqrt{M^2-1}{(1+\frac{\gamma -1}{2}{M}^2)}^{-1}\frac{dM}{M}$

(Eq. 333)

The integral on the right hand side of Eqn. \ref{eq:mach:turning:c} is the Prandtl-Meyer function, which is usually denoted $\nu$. The Prandtl-Meyer function evaluated for Mach number $M$ becomes

$\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}$

(Eq. 334)

and thus the net turning of the flow can be calculated as

${\theta }_2-{\theta }_1=\nu (M_2)-\nu (M_1)$

(Eq. 335)

Solving Problems using the Prandtl Meyer Function

A typical problem is one where we know the net flow turning and the upstream flow conditions and want to calculate the flow conditions downstream of the expansion region. An example of such a problem is given in Fig. \ref{fig:expansion:corner}.

A problem of that type can be solved as follows:

1. Calculate $\nu M_1$ using Eqn. \ref{eq:prandtl:meyer} or tabulated values
2. Calculate $\nu (M_2)$ as $\nu (M_2)={\theta }_2-{\theta }_1+\nu (M_1)$
3. Calculate $M_2$ from the known $\nu M_2$ using Eqn. \ref{eq:prandtl:meyer} or tabulated values

The aim is to derive relations of temperature, pressure and density over an expansion wave. Total temperature upstream and downstream of the expansion wave is calculated as

$\frac{T_{o_1}}{T_1}=1+\frac{\gamma -1}{2}{M}_1^2$

(Eq. 336)

$\frac{T_{o_2}}{T_2}=1+\frac{\gamma -1}{2}{M}_2^2$

(Eq. 337)

The temperature ratio over the expansion wave may now be calculated as

$\frac{T_2}{T_1}=\frac{T_2}{T_{o_2}}\frac{T_{o_1}}{T_1}=\frac{1+\frac{\gamma -1}{2}{M}_1^2}{1+\frac{\gamma -1}{2}{M}_2^2}=\frac{2+(\gamma -1)M_1^2}{2+(\gamma -1)M_2^2}$

(Eq. 338)

The expansion is isentropic and thus total temperature and both total pressure are unaffected by the expansion. Therefore, $T_{o_1}=T_{o_2}$ and thus

$\frac{T_2}{T_1}=\frac{2+(\gamma -1)M_1^2}{2+(\gamma -1)M_2^2}$

(Eq. 339)

The pressure and density ratios can be obtained from Eqn. \ref{eq:tr} using the isentropic relations

$\frac{p_2}{p_1}={\left[\frac{2+(\gamma -1)M_1^2}{2+(\gamma -1)M_2^2}\right]}^{\gamma /(\gamma -1)}$

(Eq. 340)

$\frac{{\rho }_2}{{\rho }_1}={\left[\frac{2+(\gamma -1)M_1^2}{2+(\gamma -1)M_2^2}\right]}^{1/(\gamma -1)}$

(Eq. 341)

Quasi-one-dimensional flow

The Q1D equations

Governing Equations

In the following quasi-one-dimensional flow will be assumed. That means that the cross-section is allowed to vary smoothly but flow quantities varies in one direction only. The equations that are derived will thus describe one-dimensional flow in axisymmetric tubes. Let's assume flow in the $x$-direction, which means that all flow quantities and the cross-section area will vary with the axial coordinate $x$.

$A=A(x),\rho =\rho (x),u=u(x),p=p(x),\mathrm{...}$

(Eq. 342)

We will further assume steady-state flow, which means that unsteady terms will be zero.

The equations are derived with the starting point in the governing flow equations on integral form

Continuity Equation

Applying the integral form of the continuity equation on the quasi-one-dimensional flow control volume (Fig. \ref{fig:cv}) gives

$\underset{=0}{\underbrace{\frac{d}{dt}\underset{\Omega }{\iiint }\rho dV}}+\underset{\partial \Omega }{\iint }\rho 𝐯\cdot 𝐧dS=0$

(Eq. 343)

$\underset{\partial \Omega }{\iint }\rho 𝐯\cdot 𝐧dS=-{\rho }_1{u}_1{A}_1+{\rho }_2{u}_2{A}_2$

(Eq. 344)

${\rho }_1{u}_1{A}_1={\rho }_2{u}_2{A}_2$

(Eq. 345)

Momentum Equation

Applying the integral form of the momentum equation on the quasi-one-dimensional flow control volume (Fig. \ref{fig:cv}) gives

$\underset{=0}{\underbrace{\frac{d}{dt}\underset{\Omega }{\iiint }\rho 𝐯dV}}+\underset{\partial \Omega }{\iint }[\rho (𝐯\cdot 𝐧)𝐯+p𝐧]dS=0$

(Eq. 346)

$\underset{\partial \Omega }{\iint }\rho (𝐯\cdot 𝐧)𝐯dS=-{\rho }_1{u}_1^2{A}_1+{\rho }_2{u}_2^2{A}_2$

(Eq. 347)

$\underset{\partial \Omega }{\iint }p𝐧dS=-p_1{A}_1+p_2{A}_2-\underset{A_1}{\overset{A_2}{\int }}pdA$

(Eq. 348)

collecting terms

$({\rho }_1{u}_1^2+p_1)A_1+\underset{A_1}{\overset{A_2}{\int }}pdA=({\rho }_2{u}_2^2+p_2)A_2$

(Eq. 349)

Energy Equation

Applying the integral form of the energy equation on the quasi-one-dimensional flow control volume (Fig. \ref{fig:cv}) gives

$\underset{=0}{\underbrace{\frac{d}{dt}\underset{\Omega }{\iiint }\rho e_{o}dV}}+\underset{\partial \Omega }{\iint }[\rho h_o(𝐯\cdot 𝐧)]dS=0$

(Eq. 350)

$\underset{\partial \Omega }{\iint }[\rho h_o(𝐯\cdot 𝐧)]dS=-{\rho }_1{u}_1{h}_{o_1}{A}_1+{\rho }_2{u}_2{h}_{o_2}{A}_2$

(Eq. 351)

${\rho }_1{u}_1{h}_{o_1}{A}_1={\rho }_2{u}_2{h}_{o_2}{A}_2$

(Eq. 352)

Now, using the continuity equation ${\rho }_1{u}_1{A}_1={\rho }_2{u}_2{A}_2$ gives

$h_{o_1}=h_{o_2}$

(Eq. 353)

Differential Form

The integral term appearing the momentum equation is undesired and therefore the governing equations are converted to differential form.

The continuity equation (Eqn. \ref{eq:governing:cont:b}) is rewritten in differential form as

${\rho }_1{u}_1{A}_1={\rho }_2{u}_2{A}_2=const$

(Eq. 354)

$d(\rho uA)=0$

(Eq. 355)

The momentum equation (Eqn. \ref{eq:governing:mom:b}) is rewritten in differential form as

$({\rho }_1{u}_1^2+p_1)A_1+\underset{A_1}{\overset{A_2}{\int }}pdA=({\rho }_2{u}_2^2+p_2)A_2\Rightarrow d[(\rho u^2+p)A]=pdA$

(Eq. 356)

$d(\rho u^{2}A)+d(pA)=pdA$

(Eq. 357)

$ud(\rho uA)+\rho uAdu+Adp+\cancel{pdA}=\cancel{pdA}$

(Eq. 358)

From the continuity equation we have $d(\rho uA)$ and thus

$\rho u\cancel{A}du+\cancel{A}dp=0\Rightarrow$

(Eq. 359)

$dp=-\rho udu$

(Eq. 360)

which is the momentum equation on differential form. Also referred to as Euler's equation. Finally, the energy equation (Eqn. \ref{eq:governing:cont:b}) is rewritten in differential form as

$h_{o_1}=h_{o_2}=const\Rightarrow d{h}_o=0$

(Eq. 361)

$h_o=h+\frac{1}{2}{u}^2\Rightarrow dh+\frac{1}{2}d(u^2)=0$

(Eq. 362)

$dh+udu=0$

(Eq. 363)

Summary

Continuity:

$d(\rho uA)=0$

Momentum:

$dp=-\rho udu$

Energy:

$dh+udu=0$

The equations are valid for:

• quasi-one-dimensional flow
• steady state
• all gas models (no gas model assumptions made)
• inviscid flow

It should be noted that equations are exact but they are applied to a physical model that is approximate, i.e., the approximation that flow quantities varies in one dimension with a varying cross-section area. In reality, a variation of cross-section area would imply flow in three dimensions.

Area-velocity relation

The Area-Velocity Relation

Starting point - the continuity equation (Eqn. \ref{eq:governing:cont}):

$d(\rho uA)=0\Rightarrow \rho udA+\rho Adu+uAd\rho =0$

(Eq. 364)

divide by $\rho uA$ gives

$\frac{d\rho }{\rho }+\frac{du}{u}+\frac{dA}{A}=0$

(Eq. 365)

As the name suggests, the area-velocity relation is a relation including the area and the flow velocity. Therefore, the next step is to replace the density terms.

This can be achieved using the momentum equation (Eqn. \ref{eq:governing:mom})

$dp=-\rho udu\iff \frac{dp}{\rho }=-udu$

(Eq. 366)

$\frac{dp}{\rho }=\frac{dp}{d\rho }\frac{d\rho }{\rho }=-udu$

(Eq. 367)

If we assume adiabatic and reversible flow processes, i.e., isentropic flow

$\frac{dp}{d\rho }={\left(\frac{dp}{d\rho }\right)}_s=a^2\Rightarrow a^2\frac{d\rho }{\rho }=-udu$

(Eq. 368)

$a^2\frac{d\rho }{\rho }=-udu=-u^2\frac{du}{u}$

(Eq. 369)

$\frac{d\rho }{\rho }=-M^2\frac{du}{u}$

(Eq. 370)

Eqn. \ref{eq:governing:mom:b} inserted in Eqn. \ref{eq:governing:cont:b} gives

$-M^2\frac{du}{u}+\frac{du}{u}+\frac{dA}{A}=0$

(Eq. 371)

or

$\frac{dA}{A}=(M^2-1)\frac{du}{u}$

(Eq. 372)

which is the area-velocity relation.

From the area-velocity relation (Eqn. \ref{eq:governing:av}), we can learn that in a subsonic flow, the flow will accelerate if the cross-section area is decreased and decelerate if the cross-section area is increased. It can also be seen that for supersonic flow, the relation between flow velocity and cross-section area will be the opposite of that for subsonic flows, see Fig. \ref{fig:areavelocity}. For sonic flow, $M=1$, the relation shows that $dA=0$, which means that sonic flow can only occur at a cross-section area maximum or minimum. From the subsonic versus supersonic flow discussion, it can be understood that sonic flow at the minimum cross section area is the only valid option (see Fig. \ref{fig:sonic}).

Area-Mach relation

The Area-Mach-Number Relation

Starting point - the continuity equation (Eqn. \ref{eq:governing:cont}):

$d(\rho uA)=0\Rightarrow \rho uA=const$

(Eq. 373)

This applies everywhere in the nozzle and therefore the sonic conditions can be used as a reference

$\rho uA={\rho }^{*}{u}^{*}{A}^{*}=\{u^{*}=a^{*}\}={\rho }^{*}{a}^{*}{A}^{*}$

(Eq. 374)

divide by $\rho u{A}^{*}$ gives

$\frac{{\rho }^{*}}{\rho }\frac{a^{*}}{u}=\frac{A}{A^{*}}$

(Eq. 375)

$a^{*}/u=1/M^{*}$ but ${\rho }^{*}/\rho$ is unknown

$\frac{{\rho }^{*}}{\rho }=\frac{{\rho }^{*}}{{\rho }_o}\frac{{\rho }_o}{\rho }$

(Eq. 376)

and thus

$\frac{{\rho }^{*}}{{\rho }_o}\frac{{\rho }_o}{\rho }\frac{1}{M^{*}}=\frac{A}{A^{*}}$

(Eq. 377)

Using the isentropic relations, we get

$\frac{{\rho }^{*}}{{\rho }_o}=\frac{1}{{[\frac{1}{2}(\gamma -1)]}^{1/(\gamma -1)}}$

(Eq. 378)

$\frac{{\rho }_o}{\rho }={[1+\frac{1}{2}(\gamma +1)M^2]}^{1/(\gamma -1)}$

(Eq. 379)

Eqns. \ref{eq:rho:a} and \ref{eq:rho:b} in Eqn. \ref{eq:areamach:a} gives

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

(Eq. 380)

What remains now is to replace $M^{*}$

${M^{*}}^2=\frac{u^2}{{a^{*}}^2}=\frac{u^2}{a^2}\frac{a^2}{{a^{*}}^2}=\frac{u^2}{a^2}\frac{a^2}{a_o^2}\frac{a_o^2}{{a^{*}}^2}=M^2\frac{a^2}{a_o^2}\frac{a_o^2}{{a^{*}}^2}$

(Eq. 381)

For a calorically perfect gas $a=\sqrt{\gamma RT}$, which gives

$\frac{a^2}{a_o^2}=\frac{T}{T_o}={[1+\frac{1}{2}(\gamma -1)M^2]}^{-1}$

(Eq. 382)

$\frac{a_o^2}{{a^{*}}^2}=\frac{T_o}{T^{*}}=\frac{1}{2}(\gamma +1)$

(Eq. 383)

Eqns. \ref{eq:a:a} and \ref{eq:a:b} in Eqn. \ref{eq:mstar:a} gives

${M^{*}}^2=\frac{(\gamma +1)M^2}{2+(\gamma -1)M^2}$

(Eq. 384)

Now, rewrite Eqn. \ref{eq:areamach:b} as

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

(Eq. 385)

and insert ${M^{*}}^2$ from Eqn. \ref{eq:mstar:b}

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

(Eq. 386)

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

(Eq. 387)

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

(Eq. 388)

which is the area-Mach-number relation.

For a nozzle flow, the area-Mach-number relation gives the Mach number, $M$, at any location inside the nozzle as a function of the ratio between the local cross-section area, $A$, and the throat area at choked conditions, $A^{*}$.

$M=f\left(\frac{A}{A^{*}}\right)$

(Eq. 389)

Due to the assumptions made in the derivation, the area-Mach-number relation is only valid for isentropic flows of calorically perfect gases. This means that it cannot be used throughout the divergent part of a convergent-divergent nozzle in case there is a shock within the nozzle. It can, however, be used both upstream and downstream of the shock. Note that $A^{*}$ will change over the shock.

Choked flow

Geometric Choking

For steady-state nozzle flow, the massflow is obtained as

$\dot{m}=\rho uA=const$

(Eq. 390)

Eqn. \ref{eq:massflow:a} can be evaluated at any location inside the nozzle and if evaluated at sonic conditions we get

$\dot{m}={\rho }^{*}{u}^{*}{A}^{*}$

(Eq. 391)

By definition $u^{*}=a^{*}$ and thus

$\dot{m}={\rho }^{*}{a}^{*}{A}^{*}$

(Eq. 392)

${\rho }^{*}$ and $a^{*}$ can be obtained using the ratios ${\rho }^{*}/{\rho }_o$ and $a^{*}/a_o$

${\rho }^{*}=\left(\frac{{\rho }^{*}}{{\rho }_o}\right){\rho }_o=\frac{p_o}{R{T}_o}{\left(\frac{2}{\gamma +1}\right)}^{1/(\gamma -1)}$

(Eq. 393)

$a^{*}=\left(\frac{a^{*}}{a_o}\right)a_o=a_o{\left(\frac{2}{\gamma +1}\right)}^{1/2}=\sqrt{\gamma R{T}_o}{\left(\frac{2}{\gamma +1}\right)}^{1/2}$

(Eq. 394)

Eqns. \ref{eq:as} and \ref{eq:rhos} in Eqn. \ref{eq:massflow:b} gives

$\dot{m}=\frac{p_o}{R{T}_o}{\left(\frac{2}{\gamma +1}\right)}^{1/(\gamma -1)}\sqrt{\gamma R{T}_o}{\left(\frac{2}{\gamma +1}\right)}^{1/2}{A}^{*}$

(Eq. 395)

which can be rewritten as

$\dot{m}=\frac{p_o{A}^{*}}{\sqrt{T_o}}\sqrt{\frac{\gamma }{R}{\left(\frac{2}{\gamma +1}\right)}^{(\gamma +1)/(\gamma -1)}}$

(Eq. 396)

Eqn. \ref{eq:massflow:c} valid for:

• quasi-one-dimensional flow
• steady state
• inviscid flow
• calorically perfect gas

It should be noted that the choked massflow can be calculated using Eqn. \ref{eq:massflow:c} even for cases with shocks downstream of the throat.

Nozzle flow

add description of nozzle flows here...

Diffusers

Add description and examples here...

Unsteady waves

Moving shock waves

Moving Normal Shock Waves

The starting point is the governing equations for stationary normal shocks (repeated here for convenience).

${\rho }_1{u}_1={\rho }_2{u}_2$

(Eq. 397)

${\rho }_1{u}_1^2+p_1={\rho }_2{u}_2^2+p_2$

(Eq. 398)

$h_1+\frac{1}{2}{u}_1^2=h_2+\frac{1}{2}{u}_2^2$

(Eq. 399)

Shock moving to the right with the constant speed $W$ into a gas that is standing still. Moving with the shock, we would see a gas velocity ahead of the shock $u_1=W$, and the gas behind the shock moves to the right with the velocity $u_2=W-u_p$. Now, let's insert $u_1$ and $u_2$ in the stationary shock relations \ref{eq:stationary:cont} - \ref{eq:stationary:energy}.

${\rho }_{1}W={\rho }_2(W-u_p)$

(Eq. 400)

${\rho }_1{W}^2+p_1={\rho }_2(W-u_p{)}^2+p_2$

(Eq. 401)

$h_1+\frac{1}{2}{W}^2=h_2+\frac{1}{2}(W-u_p{)}^2$

(Eq. 402)

Rewriting Eqn. \ref{eq:unsteady:cont}

$(W-u_p)=W\frac{{\rho }_1}{{\rho }_2}$

(Eq. 403)

Inserting Eqn. \ref{eq:unsteady:cont:mod} in Eqn. \ref{eq:unsteady:mom} gives

$p_1+{\rho }_1{W}^2=p_2+{\rho }_2{W}^2{\left(\frac{{\rho }_1}{{\rho }_2}\right)}^2\Rightarrow p_2-p_1={\rho }_1{W}^2(1-\frac{{\rho }_1}{{\rho }_2})$

(Eq. 404)

$W^2=\frac{p_2-p_1}{{\rho }_2-{\rho }_1}\left(\frac{{\rho }_2}{{\rho }_1}\right)$

(Eq. 405)

From the continuity equation \ref{eq:unsteady:cont}, we get

$W=(W-u_p)\left(\frac{{\rho }_2}{{\rho }_1}\right)$

(Eq. 406)

Inserting Eqn. \ref{eq:unsteady:cont:modb} in Eqn. \ref{eq:unsteady:mom:mod} gives

$(W-u_p{)}^2=\frac{p_2-p_1}{{\rho }_2-{\rho }_1}\left(\frac{{\rho }_1}{{\rho }_2}\right)$

(Eq. 407)

Now, let's insert Eqns. \ref{eq:unsteady:mom:mod} and \ref{eq:unsteady:mom:modb} in the energy equation (Eqn. \ref{eq:unsteady:energy}).

$h_1+\frac{1}{2}\left[\frac{p_2-p_1}{{\rho }_2-{\rho }_1}\left(\frac{{\rho }_2}{{\rho }_1}\right)\right]=h_2+\frac{1}{2}\left[\frac{p_2-p_1}{{\rho }_2-{\rho }_1}\left(\frac{{\rho }_1}{{\rho }_2}\right)\right]$

(Eq. 408)

$h=e+\frac{p}{\rho }$

(Eq. 409)

$e_1+\frac{p_1}{{\rho }_1}+\frac{1}{2}\left[\frac{p_2-p_1}{{\rho }_2-{\rho }_1}\left(\frac{{\rho }_2}{{\rho }_1}\right)\right]=e_2+\frac{p_2}{{\rho }_2}+\frac{1}{2}\left[\frac{p_2-p_1}{{\rho }_2-{\rho }_1}\left(\frac{{\rho }_1}{{\rho }_2}\right)\right]$

(Eq. 410)

which can be rewritten as

$e_2-e_1=\frac{p_1+p_2}{2}(\frac{1}{{\rho }_1}-\frac{1}{{\rho }_2})$

(Eq. 411)

Eqn \ref{eq:unsteady:hugonoit} is the same Hugoniot equation as we get for a stationary normal shock. The Hugoniot equation is a relation of thermodynamic properties over a shock. As the shock in the unsteady case is moving with a constant velocity, the frame of reference moving with the shock is an inertial frame and thus the same physical relations apply in the moving shock case as in the stationary shock case. The fact that the Hugoniot relation does not include any velocities or Mach numbers but only thermodynamic properties, the relation will be unchanged for a moving shock.

Moving Shock Relations

For a calorically perfect gas we have $e=C_{v}T$. Inserted in the Hugoniot relation above this gives

$C_v(T_2-T_1)=\frac{p_1+p_2}{2}({\nu }_1-{\nu }_2)$

(Eq. 412)

where $\nu =1/\rho$

Now, using the ideal gas law $T=p\nu /R$ and $C_v/R=1/(\gamma -1)$ gives

$\left(\frac{1}{\gamma -1}\right)(p_2{\nu }_2-p_1{\nu }_1)=\frac{p_1+p_2}{2}({\nu }_1-{\nu }_2)$ $\iff$ $p_2(\frac{{\nu }_2}{\gamma -1}-\frac{{\nu }_1-{\nu }_2}{2})=p_1(\frac{{\nu }_1}{\gamma -1}+\frac{{\nu }_1-{\nu }_2}{2})$

(Eq. 413)

From this result, we can derive a relation for the pressure ratio over the shock as a function of density ratio

$\frac{p_2}{p_1}=\frac{\left(\frac{\gamma +1}{\gamma -1}\right)\left(\frac{{\nu }_1}{{\nu }_2}\right)-1}{\left(\frac{\gamma +1}{\gamma -1}\right)-\left(\frac{{\nu }_1}{{\nu }_2}\right)}$

(Eq. 414)

$\nu =RT/p$ and thus

$\frac{{\nu }_1}{{\nu }_2}=\frac{T_1}{T_2}\frac{p_2}{p_1}$

(Eq. 415)

Eqn. \ref{eq:unsteady:density:ratio} in Eqn. \ref{eq:unsteady:hugonoit:c} gives

$\frac{p_2}{p_1}=\frac{\left(\frac{\gamma +1}{\gamma -1}\right)\left(\frac{T_1}{T_2}\frac{p_2}{p_1}\right)-1}{\left(\frac{\gamma +1}{\gamma -1}\right)-\left(\frac{T_1}{T_2}\frac{p_2}{p_1}\right)}$

(Eq. 416)

Now, we can get a relation for calculation of the temperature ratio over the moving shock as function of the shock pressure ratio

$\frac{T_2}{T_1}=\frac{p_2}{p_1}\left[\frac{\left(\frac{\gamma +1}{\gamma -1}\right)+\left(\frac{p_2}{p_1}\right)}{1+\left(\frac{\gamma +1}{\gamma -1}\right)\left(\frac{p_2}{p_1}\right)}\right]$

(Eq. 417)

Once again using the ideal gas law

$\frac{{\rho }_2}{{\rho }_1}=\frac{\left(\frac{\gamma +1}{\gamma -1}\right)+\left(\frac{p_2}{p_1}\right)}{1+\left(\frac{\gamma +1}{\gamma -1}\right)\left(\frac{p_2}{p_1}\right)}$

(Eq. 418)

Going back to the momentum equation

$p_2-p_1={\rho }_1{W}^2(1-\frac{{\rho }_1}{{\rho }_2})=\{W=M_s{a}_1\}={\rho }_1{M}_s^2{a}_1^2(1-\frac{{\rho }_1}{{\rho }_2})$

(Eq. 419)

with $a_1^2=\gamma p_1/{\rho }_1$, we get

$\frac{p_2}{p_1}=\gamma M_s^2(1-\frac{{\rho }_1}{{\rho }_2})+1$

(Eq. 420)

From the normal shock relations, we have

$\frac{{\rho }_1}{{\rho }_2}=\frac{2+(\gamma -1)M_s^2}{(\gamma +1)M_s^2}$

(Eq. 421)

Eqn. \ref{eq:unsteady:Mach:b} in \ref{eq:unsteady:Mach:a} gives

$\frac{p_2}{p_1}=1+\left(\frac{2\gamma }{\gamma +1}\right)(M_s^2-1)$

(Eq. 422)

or

$M_s=\sqrt{\left(\frac{\gamma +1}{2\gamma }\right)(\frac{p_2}{p_1}-1)+1}$

(Eq. 423)

Eqn. \ref{eq:unsteady:Mach} with $M_s=W/a_1$

$W=a_1\sqrt{\left(\frac{\gamma +1}{2\gamma }\right)(\frac{p_2}{p_1}-1)+1}$

(Eq. 424)

Induced Flow Behind Moving Shock

Let's try to find a relation for calculation of the induced velocity behind the moving shock. Once again, the starting point is the continuity equation for moving shocks (Eqn. \ref{eq:unsteady:cont}) repeated here for convenience

${\rho }_{1}W={\rho }_2(W-u_p)$

(Eq. 425)

The induced velocity appears on the right side of the continuity equation

$W({\rho }_1-{\rho }_2)=-{\rho }_2{u}_p$

(Eq. 426)

$u_p=W(1-\frac{{\rho }_1}{{\rho }_2})$

(Eq. 427)

From before we have a relation for $W$ as a function of pressure ratio and one for ${\rho }_1/{\rho }_2$, also as a function of pressure ratio.

Eqn. \ref{eq:unsteady:up:a} togheter with Eqns. \ref{eq:unsteady:W} and \ref{eq:unsteady:density:ratio} gives

$u_p=a_1\underset{I}{\underbrace{\sqrt{\left(\frac{\gamma +1}{2\gamma }\right)(\frac{p_2}{p_1}-1)+1}}}\underset{II}{\underbrace{[1-\frac{\left(\frac{\gamma +1}{\gamma -1}\right)+\left(\frac{p_2}{p_1}\right)}{1+\left(\frac{\gamma +1}{\gamma -1}\right)\left(\frac{p_2}{p_1}\right)}]}}$

(Eq. 428)

The equation subsets I and II can be rewritten as:

Term I:

$\sqrt{\left(\frac{\gamma +1}{2\gamma }\right)(\frac{p_2}{p_1}-1)+1}=\sqrt{\frac{\gamma +1}{2\gamma }[\left(\frac{p_2}{p_1}\right)+\left(\frac{\gamma -1}{\gamma +1}\right)]}$

(Eq. 429)

Term II:

$[1-\frac{\left(\frac{\gamma +1}{\gamma -1}\right)+\left(\frac{p_2}{p_1}\right)}{1+\left(\frac{\gamma +1}{\gamma -1}\right)\left(\frac{p_2}{p_1}\right)}]=\frac{1}{\gamma }(\frac{p_2}{p_1}-1)\frac{\left(\frac{2\gamma }{\gamma +1}\right)}{\left(\frac{\gamma -1}{\gamma +1}\right)+\left(\frac{p_2}{p_1}\right)}$

(Eq. 430)

the rewritten terms I and II implemented, Eqn. \ref{eq:unsteady:up:b} becomes

$u_p=\frac{a_1}{\gamma }(\frac{p_2}{p_1}-1)\sqrt{\frac{\left(\frac{2\gamma }{\gamma +1}\right)}{\left(\frac{\gamma -1}{\gamma +1}\right)+\left(\frac{p_2}{p_1}\right)}}$

(Eq. 431)

Since the region behind the moving shock is region 2, the induced flow Mach number is obtained as

$M_p=\frac{u_p}{a_2}=\frac{u_p}{a_1}\frac{a_1}{a_2}=\frac{u_p}{a_1}\sqrt{\frac{\gamma R{T}_1}{\gamma R{T}_2}}=\frac{u_p}{a_1}\sqrt{\frac{T_1}{T_2}}$

(Eq. 432)

With $up/a_1$ from Eqn. \ref{eq:unsteady:up} and $T_1/T_2$ from Eqn. \ref{eq:unsteady:temperature:ratio}

$M_p=\frac{1}{\gamma }(\frac{p_2}{p_1}-1){\left(\frac{\left(\frac{2\gamma }{\gamma +1}\right)}{\left(\frac{\gamma -1}{\gamma +1}\right)+\left(\frac{p_2}{p_1}\right)}\right)}^{1/2}{\left(\frac{1+\left(\frac{\gamma +1}{\gamma -1}\right)\left(\frac{p_2}{p_1}\right)}{\left(\frac{\gamma +1}{\gamma -1}\right)\left(\frac{p_2}{p_1}\right)+{\left(\frac{p_2}{p_1}\right)}^2}\right)}^{1/2}$

(Eq. 433)

There is a theoretical upper limit for the induced Mach number $M_p$

$\underset{p_2/p_1\to \mathrm{\infty }}{\lim }{M}_p\left(\frac{p_2}{p_1}\right)=\sqrt{\frac{2}{\gamma (\gamma -1)}}$

(Eq. 434)

As can be seen, at the upper limit the induced Mach number is a function of $\gamma$ and for air ($\gamma =1.4$) we get

$\underset{p_2/p_1\to \mathrm{\infty }}{\lim }{M}_p\left(\frac{p_2}{p_1}\right)\simeq 1.89$

(Eq. 435)

Shock Wave Reflection

When the incident shock wave reaches the wall, a shock propagating in the opposite direction is generated with a shock strength such that the velocity of the induced flow behind the incident shock is reduced to zero. The flow can not go through the wall and thus the velocity must be zero in the vicinity of the wall. The properties of the incident shock wave are directly related to the pressure ratio over the shock wave. Therefore, it would be convenient to have a relation between the reflected shock wave and incident shock wave.

The Incident Shock Wave

The pressure ratio over the incident shock in Fig.~\ref{fig:reflection} can be obtained as

$\frac{p_2}{p_1}=1+\frac{2\gamma }{\gamma +1}(M_s^2-1)$

(Eq. 436)

where $M_s$ is the wave Mach number, which is calculated as

$M_s=\frac{W}{a_1}$

(Eq. 437)

In Eqn.~\ref{eq:incident:Mach:def}, $W$ is the speed with which the incident shock wave travels into region 1 and $a_1$ is the speed of sound in region 1 (see Fig.~\ref{fig:reflection}).

Solving Eqn.~\ref{eq:incident:pr} for $M_s$, we get

$M_s=\sqrt{\frac{\gamma +1}{2\gamma }(\frac{p_2}{p_1}-1)+1}$

(Eq. 438)

Anderson derives the relations for calculation of the ratio $T_2/T_1$

$\frac{T_2}{T_1}=\frac{p_2}{p_1}\left(\frac{\frac{\gamma +1}{\gamma -1}+\frac{p_2}{p_1}}{1+\frac{\gamma +1}{\gamma -1}\frac{p_2}{p_1}}\right)$

(Eq. 439)

From Eqn.~\ref{eq:incident:tr} it is easy to get the corresponding relation for ${\rho }_2/{\rho }_1$

$\frac{{\rho }_2}{{\rho }_1}=\frac{1+\frac{\gamma +1}{\gamma -1}\frac{p_2}{p_1}}{\frac{\gamma +1}{\gamma -1}+\frac{p_2}{p_1}}$

(Eq. 440)

Anderson also shows how to obtain the induced velocity, $u_p$, behind the incident shock wave, {\emph{i.e.}} the velocity in region 2 (see Fig.~\ref{fig:reflection}).

$u_p=W(1-\frac{{\rho }_1}{{\rho }_2})=M_s{a}_1(1-\frac{{\rho }_1}{{\rho }_2})$

(Eq. 441)

The Reflected Shock Wave

The pressure ratio over the reflected shock can be obtained from Eqn.~\ref{eq:incident:pr} by analogy

$\frac{p_5}{p_2}=1+\frac{2\gamma }{\gamma +1}(M_r^2-1)$

(Eq. 442)

where $M_r$ is the Mach number of the reflected shock wave defined as

$M_r=\frac{W_r+u_p}{a_2}$

(Eq. 443)

where $W_r$ is the speed of the reflected shock wave and $a_2$ is the speed of sound in region 2 (see Fig.~\ref{fig:reflection}).

Solving Eqn.~\ref{eq:reflected:pr} for $M_r$ gives

$M_r=\sqrt{\frac{\gamma +1}{2\gamma }(\frac{p_5}{p_2}-1)+1}$

(Eq. 444)

The ratios $T_5/T_2$ and ${\rho }_5/{\rho }_2$ can be obtained from Eqns.~\ref{eq:incident:tr} and \ref{eq:incident:rr} by analogy

$\frac{T_5}{T_2}=\frac{p_5}{p_2}\left(\frac{\frac{\gamma +1}{\gamma -1}+\frac{p_5}{p_2}}{1+\frac{\gamma +1}{\gamma -1}\frac{p_5}{p_2}}\right)$

(Eq. 445)

$\frac{{\rho }_5}{{\rho }_2}=\frac{1+\frac{\gamma +1}{\gamma -1}\frac{p_5}{p_2}}{\frac{\gamma +1}{\gamma -1}+\frac{p_5}{p_2}}$

(Eq. 446)

The velocity in region 2 which is the same as the induced flow velocity behind the incident shock wave can be obtained as

$u_p=W_r(\frac{{\rho }_5}{{\rho }_2}-1)=M_r{a}_2(1-\frac{{\rho }_2}{{\rho }_5})$

(Eq. 447)

Reflected Shock Relation

With the relations for the incident shock wave and reflected shock wave defined, we now have the tools to derive a relation between the incident and reflected shock waves. The induced flow velocity $u_p$ calculated using the relation obtained for the incident shock wave must of course be the same as when calculated using reflected wave properties, {\emph{i.e.}} the result of Eqn.~\ref{eq:incident:up} is identical to that of Eqn.~\ref{eq:reflected:up}

$M_r{a}_2(1-\frac{{\rho }_2}{{\rho }_5})=M_s{a}_1(1-\frac{{\rho }_1}{{\rho }_2})$

(Eq. 448)

rewriting gives

$M_r(1-\frac{{\rho }_2}{{\rho }_5})=M_s(1-\frac{{\rho }_1}{{\rho }_2})\frac{a_1}{a_2}$

(Eq. 449)

Assuming calorically perfect gas gives $a=\sqrt{\gamma RT}$ and thus

$M_r(1-\frac{{\rho }_2}{{\rho }_5})=M_s(1-\frac{{\rho }_1}{{\rho }_2})\sqrt{\frac{T_1}{T_2}}$

(Eq. 450)

Let's first look at the term on the left hand side of Eqn.~\ref{eq:relation:c}

$M_r(1-\frac{{\rho }_2}{{\rho }_5})$

(Eq. 451)

Using the ${\rho }_5/{\rho }_2$ and $p_2/p_5$ from Eqns.~\ref{eq:reflected:rr} and~\ref{eq:reflected:pr} and simplifying gives

$M_r(1-\frac{{\rho }_2}{{\rho }_5})=\left(\frac{2}{\gamma +1}\right)\left(\frac{M_r^2-1}{M_r}\right)$

(Eq. 452)

Using the same approach on the corresponding term for the incident shock wave on the right hand side of Eqn.~\ref{eq:relation:c} gives

$M_s(1-\frac{{\rho }_1}{{\rho }_2})=\left(\frac{2}{\gamma +1}\right)\left(\frac{M_s^2-1}{M_s}\right)$

(Eq. 453)

Now, inserting~\ref{eq:relation:d} and~\ref{eq:relation:e} in Eqn.~\ref{eq:relation:c} gives

$\left(\frac{2}{\gamma +1}\right)\left(\frac{M_r^2-1}{M_r}\right)=\left(\frac{2}{\gamma +1}\right)\left(\frac{M_s^2-1}{M_s}\right)\sqrt{\frac{T_1}{T_2}}$

(Eq. 454)

Simplifying and inverting gives

$\left(\frac{M_r}{M_r^2-1}\right)=\left(\frac{M_s}{M_s^2-1}\right)\sqrt{\frac{T_2}{T_1}}$

(Eq. 455)

The rightmost term in Eqn.~\ref{eq:relation:g} ($\sqrt{T_2/T_1}$) needs to be rewritten. Inserting~\ref{eq:incident:pr} in~\ref{eq:incident:tr} and expanding all terms gives

$\frac{T_2}{T_1}=\frac{2(\gamma +1)+(\gamma +1)(\gamma -1)M_s^2+4\gamma (M_s^2-1)+2\gamma (\gamma -1)M_s^2(M_s^2-1)}{(\gamma +1{)}^2{M}_s^2}=$

$=\frac{2(\gamma +1)+(\gamma +1)(\gamma -1)M_s^2+4\gamma (M_s^2-1)}{(\gamma +1{)}^2{M}_s^2}+\frac{2(\gamma -1)}{(\gamma +1{)}^2}(M_s^2-1)\gamma =$

$=\frac{2(\gamma +1)+(\gamma +1)(\gamma -1)M_s^2+4\gamma (M_s^2-1)-(2(\gamma -1)(M_s^2-1))}{(\gamma +1{)}^2{M}_s^2}$

$\frac{2(\gamma -1)}{(\gamma +1{)}^2}(M_s^2-1)(\gamma +\frac{1}{M_s^2})$

(Eq. 456)

Finally we end up with the following relation

$\frac{T_2}{T_1}=1+\frac{2(\gamma -1)}{(\gamma +1{)}^2}(M_s^2-1)(\gamma +\frac{1}{M_s^2})$

(Eq. 457)

The temperature ratio over the incident shock wave is now totally defined by the incident Mach number $M_s$ and the ratio of specific heats $\gamma$. With~\ref{eq:relation:tr} in~\ref{eq:relation:g} we get the sought relation between the reflected and incident Mach numbers.

$\left(\frac{M_r}{M_r^2-1}\right)=\left(\frac{M_s}{M_s^2-1}\right)\sqrt{1+\frac{2(\gamma -1)}{(\gamma +1{)}^2}(M_s^2-1)(\gamma +\frac{1}{M_s^2})}$

(Eq. 458)

It should be noted that Eqn.~\ref{eq:relation:final} is valid for calorically perfect gases only.

Acoustic theory

In the following we are going to derive the linear acoustic wave equation starting from the continuity and momentum equations on non-conservation differential form. The equations are repeated here for convenience.

$\frac{D\rho }{Dt}+\rho (\mathrm{\nabla }\cdot 𝐯)=0$

(Eq. 459)

$\rho \frac{D𝐯}{Dt}+\mathrm{\nabla }p=0$

(Eq. 460)

Remember that $D/Dt$ denotes the substantial derivative operator defined as follows

$\frac{D}{Dt}=\frac{\partial }{\partial t}+𝐯\cdot \mathrm{\nabla }$

(Eq. 461)

where $\partial /\partial t$ is the local temporal derivative and $𝐯\cdot \mathrm{\nabla }$ is the convective derivative.

We are going to analyze acoustic waves in one dimension, which means that the equations above reduces to

$\frac{\partial \rho }{\partial t}+u\frac{\partial \rho }{\partial x}+\rho \frac{\partial u}{\partial x}=0$

(Eq. 462)

$\rho \frac{\partial u}{\partial t}+\rho u\frac{\partial u}{\partial x}+\frac{\partial p}{\partial x}=0$

(Eq. 463)

Pressure is a thermodynamic property and thus it can be expressed as a function of two other thermodynamic properties. Let's express pressure as a function of density ($\rho$) and entropy ($s$).

$p=p(\rho ,s)\Rightarrow dp={\left(\frac{\partial p}{\partial \rho }\right)}_{s}d\rho +{\left(\frac{\partial p}{\partial s}\right)}_{\rho }ds$

(Eq. 464)

Since weak acoustic waves are considered, entropy will be constant and thus $ds=0$, which means that

$dp={\left(\frac{\partial p}{\partial \rho }\right)}_{s}d\rho =a^{2}d\rho$

(Eq. 465)

$\rho \frac{\partial u}{\partial t}+\rho u\frac{\partial u}{\partial x}+a^2\frac{\partial \rho }{\partial x}=0$

(Eq. 466)

The acoustic perturbations can be described as small deviations around a reference state

$\begin{array}{cc}& \rho ={\rho }_{\mathrm{\infty }}+\Delta \rho \\ & p=p_{\mathrm{\infty }}+\Delta p\\ & T=T_{\mathrm{\infty }}+\Delta T\\ & u=u_{\mathrm{\infty }}+\Delta u=\{u_{\mathrm{\infty }}=0\}=\Delta u\\ \end{array}$

Inserted in Eqns.~\ref{eq:unstady:acoustic:wave:cont} and \ref{eq:unstady:acoustic:wave:mom:b} and using the fact that derivatives of the constant reference state flow quantities are zero, we get

$\frac{\partial }{\partial t}(\Delta \rho )+\Delta u\frac{\partial }{\partial x}(\Delta \rho )+({\rho }_{\mathrm{\infty }}+\Delta \rho )\frac{\partial }{\partial x}(\Delta u)=0$

(Eq. 467)

$({\rho }_{\mathrm{\infty }}+\Delta \rho )\frac{\partial }{\partial t}(\Delta u)+({\rho }_{\mathrm{\infty }}+\Delta \rho )\Delta u\frac{\partial }{\partial x}(\Delta u)+a^2\frac{\partial }{\partial x}(\Delta \rho )=0$

(Eq. 468)

In the same way as pressure, being a thermodynamic variable, can be expressed as a function of two other thermodynamic variables, so can the speed of sound. Once again we will select density and entropy as the two thermodynamic variables

$a^2=a^2(\rho ,s)$

(Eq. 469)

and since entropy is constant

$a^2=a^2(\rho )$

(Eq. 470)

Taylor expansion of $a^2$ around the reference state $a_{\mathrm{\infty }}$ with $\Delta \rho =\rho -{\rho }_{\mathrm{\infty }}$ gives

$a^2=a_{\mathrm{\infty }}^2+{(\frac{\partial }{\partial \rho }(a^2))}_{\mathrm{\infty }}\Delta \rho +{(\frac{{\partial }^2}{\partial {\rho }^2}(a^2))}_{\mathrm{\infty }}(\Delta \rho {)}^2+\cdots$

(Eq. 471)

Inserted in Eqn.~\ref{eq:unstady:acoustic:wave:mom:pert}, we get

$({\rho }_{\mathrm{\infty }}+\Delta \rho )\frac{\partial }{\partial t}(\Delta u)+({\rho }_{\mathrm{\infty }}+\Delta \rho )\Delta u\frac{\partial }{\partial x}(\Delta u)+[a_{\mathrm{\infty }}^2+{(\frac{\partial }{\partial \rho }(a^2))}_{\mathrm{\infty }}\Delta \rho +\cdots ]\frac{\partial }{\partial x}(\Delta \rho )=0$

(Eq. 472)

The perturbations $\Delta u$ and $\Delta \rho$ are small, which implies that $\Delta u\ll a_{\mathrm{\infty }}$ and $\Delta \rho \ll {\rho }_{\mathrm{\infty }}$. This means that products of perturbations can be canceled and so can higher-order terms in the Taylor expansion of $a^2$. This means that the continuity and momentum equations reduces to

$\frac{\partial }{\partial t}(\Delta \rho )+{\rho }_{\mathrm{\infty }}\frac{\partial }{\partial x}(\Delta u)=0$

(Eq. 473)

${\rho }_{\mathrm{\infty }}\frac{\partial }{\partial t}(\Delta u)+a_{\mathrm{\infty }}^2\frac{\partial }{\partial x}(\Delta \rho )=0$

(Eq. 474)

Before making the assumption that the perturbations are small compared to the corresponding reference state flow quantities and thus justifying the cancelation of products of perturbations from the equations, the flow equations were still the exact fully non-linear equations. Eqns.~\ref{eq:unstady:acoustic:wave:cont:linear}. and \ref{eq:unstady:acoustic:wave:mom:linear}, however, are approximations as several terms has been removed. The equations are linear and are good approximations as long as the perturbations are small. The smaller the perturbations, the better the approximation are the linear equations. Eqns.~\ref{eq:unstady:acoustic:wave:cont:linear} and \ref{eq:unstady:acoustic:wave:mom:linear} describes the motion induced in a gas by the passage of a sound wave. By combining the temporal derivative of Eqn.~\ref{eq:unstady:acoustic:wave:cont:linear} with the divergence of Eqn.~\ref{eq:unstady:acoustic:wave:mom:linear}, it is possible to obtain a wave equation describing the propagation of acoustic waves in a quiescent surrounding.

The temporal derivative of the continuity equation:

$\frac{{\partial }^2}{\partial t^2}(\Delta \rho )+{\rho }_{\mathrm{\infty }}\frac{{\partial }^2}{\partial x\partial t}(\Delta u)=0$

(Eq. 475)

The divergence of the momentum equation:

${\rho }_{\mathrm{\infty }}\frac{{\partial }^2}{\partial x\partial t}(\Delta u)+a_{\mathrm{\infty }}^2\frac{{\partial }^2}{\partial x^2}(\Delta \rho )=0$

(Eq. 476)

The second term in the first equation is the same as the first term in the second equation. Substituting the term, the two equations reduces to one single equation

$\frac{{\partial }^2}{\partial t^2}(\Delta \rho )=a_{\mathrm{\infty }}^2\frac{{\partial }^2}{\partial x^2}(\Delta \rho )$

(Eq. 477)

which is a one-dimensional form of the classic wave equation with the general solution

$\Delta \rho =F(x-a_{\mathrm{\infty }}t)+G(x+a_{\mathrm{\infty }}t)$

(Eq. 478)

$F$ and $G$ are arbitrary functions. The function $F$ describes the shape of a wave traveling in the positive $x$-direction at the speed of sound of the ambient gas and the function $G$ describes the shape of a wave traveling in the negative $x$-direction at the same speed. In Eqn.~\ref{eq:wave} $\Delta \rho$ appears with second derivatives in space and time. Let's differentiate the proposed solution (Eqn.~\ref{eq:wave:solution}) two times in time and space, respectively, and check that it is actually a valid solution to Eqn.~\ref{eq:wave}.

$\frac{\partial }{\partial t}(\Delta \rho )=\frac{\partial F}{\partial (x-a_{\mathrm{\infty }}t)}\frac{\partial (x-a_{\mathrm{\infty }}t)}{\partial t}+\frac{\partial G}{\partial (x+a_{\mathrm{\infty }}t)}\frac{\partial (x+a_{\mathrm{\infty }}t)}{\partial t}$

(Eq. 479)

$\frac{\partial }{\partial t}(\Delta \rho )=-a_{\mathrm{\infty }}{F}^{\prime }+a_{\mathrm{\infty }}{G}^{\prime }$

(Eq. 480)

$\frac{{\partial }^2}{\partial t^2}(\Delta \rho )=-a_{\mathrm{\infty }}\frac{\partial F^{\prime }}{\partial (x-a_{\mathrm{\infty }}t)}\frac{\partial (x-a_{\mathrm{\infty }}t)}{\partial t}+a_{\mathrm{\infty }}\frac{\partial G^{\prime }}{\partial (x+a_{\mathrm{\infty }}t)}\frac{\partial (x+a_{\mathrm{\infty }}t)}{\partial t}$

(Eq. 481)

$\frac{{\partial }^2}{\partial t^2}(\Delta \rho )=a_{\mathrm{\infty }}^2{F}^{″}+a_{\mathrm{\infty }}^2{G}^{″}$

(Eq. 482)

$\frac{\partial }{\partial x}(\Delta \rho )=\frac{\partial F}{\partial (x-a_{\mathrm{\infty }}t)}\frac{\partial (x-a_{\mathrm{\infty }}t)}{\partial x}+\frac{\partial G}{\partial (x+a_{\mathrm{\infty }}t)}\frac{\partial (x+a_{\mathrm{\infty }}t)}{\partial x}$

(Eq. 483)

$\frac{\partial }{\partial x}(\Delta \rho )=F^{\prime }+G^{\prime }$

(Eq. 484)

$\frac{{\partial }^2}{\partial x^2}(\Delta \rho )=\frac{\partial F^{\prime }}{\partial (x-a_{\mathrm{\infty }}t)}\frac{\partial (x-a_{\mathrm{\infty }}t)}{\partial x}+\frac{\partial G^{\prime }}{\partial (x+a_{\mathrm{\infty }}t)}\frac{\partial (x+a_{\mathrm{\infty }}t)}{\partial x}$

(Eq. 485)

$\frac{{\partial }^2}{\partial x^2}(\Delta \rho )=F^{″}+G^{″}$

(Eq. 486)

Eqns. \ref{eq:wave:ddt} and \ref{eq:wave:ddx} inserted Eqn. \ref{eq:wave} gives

$a_{\mathrm{\infty }}^2{F}^{″}+a_{\mathrm{\infty }}^2{G}^{″}=a_{\mathrm{\infty }}^2(F^{″}+G^{″})$

(Eq. 487)

which shows that Eqn. \ref{eq:wave:solution} is a valid solution to the wave equation.

$F$ and $G$ are arbitrary functions and thus $G=0$ is a valid solution, which gives

$\Delta \rho (x,t)=F(x-a_{\mathrm{\infty }}t)$

(Eq. 488)

If $\Delta \rho$ is constant, i.e. a wave with constant amplitude, we see from Eqn.~\ref{eq:wave:solution:F} that $(x-a_{\mathrm{\infty }}t)$ is constant and thus

$x=a_{\mathrm{\infty }}t+c\Rightarrow \frac{dx}{dt}=a_{\mathrm{\infty }}$

(Eq. 489)

From Eqn.~\ref{eq:wave:solution:F}, we get

$\frac{\partial }{\partial t}(\Delta \rho )=-a_{\mathrm{\infty }}{F}^{\prime }$

(Eq. 490)

$\frac{\partial }{\partial x}(\Delta \rho )=F^{\prime }$

(Eq. 491)

and thus

$\frac{\partial }{\partial x}(\Delta \rho )=-\frac{1}{a_{\mathrm{\infty }}}\frac{\partial }{\partial t}(\Delta \rho )$

(Eq. 492)

which gives a relation between the temporal derivative of $\Delta \rho$ and the spatial derivative of $\Delta \rho$. With Eqn.~\ref{eq:wave:solution:F:b}, the linearized momentum equation Eqn.~\ref{eq:unstady:acoustic:wave:mom:linear} can be rewritten as follows

$\frac{\partial }{\partial t}(\Delta u)=-\frac{a_{\mathrm{\infty }}^2}{{\rho }_{\mathrm{\infty }}}\frac{\partial }{\partial x}(\Delta \rho )=\{\frac{\partial }{\partial x}(\Delta \rho )=-\frac{1}{a_{\mathrm{\infty }}}\frac{\partial }{\partial t}(\Delta \rho )\}=\frac{a_{\mathrm{\infty }}}{{\rho }_{\mathrm{\infty }}}\frac{\partial }{\partial t}(\Delta \rho )\Rightarrow$

(Eq. 493)

$\frac{\partial }{\partial t}(\Delta u-\frac{a_{\mathrm{\infty }}}{{\rho }_{\mathrm{\infty }}}\Delta \rho )=0\Rightarrow \Delta u-\frac{a_{\mathrm{\infty }}}{{\rho }_{\mathrm{\infty }}}\Delta \rho =const$

(Eq. 494)

In an undisturbed gas $\Delta u=\Delta \rho =0$ and thus

$\Delta u-\frac{a_{\mathrm{\infty }}}{{\rho }_{\mathrm{\infty }}}\Delta \rho =0$

(Eq. 495)

or

$\Delta u=\frac{a_{\mathrm{\infty }}}{{\rho }_{\mathrm{\infty }}}\Delta \rho$

(Eq. 496)

If instead $F$ is set to zero and $G$ is non-zero, we get

$\Delta u=-\frac{a_{\mathrm{\infty }}}{{\rho }_{\mathrm{\infty }}}\Delta \rho$

(Eq. 497)

${\left(\frac{\partial p}{\partial \rho }\right)}_s=a^2\Rightarrow \Delta p=a_{\mathrm{\infty }}^2\Delta \rho$

(Eq. 498)

Acoustic wave traveling in the positive $x$-direction:

$\Delta u=\frac{a_{\mathrm{\infty }}}{{\rho }_{\mathrm{\infty }}}\Delta \rho =\frac{1}{a_{\mathrm{\infty }}{\rho }_{\mathrm{\infty }}}\Delta p$

(Eq. 499)

Acoustic wave traveling in the negative $x$-direction:

$\Delta u=-\frac{a_{\mathrm{\infty }}}{{\rho }_{\mathrm{\infty }}}\Delta \rho =-\frac{1}{a_{\mathrm{\infty }}{\rho }_{\mathrm{\infty }}}\Delta p$

(Eq. 500)

Finite non-linear waves

Starting point: the governing flow equations on partial differential form

Continuity equation:

$\frac{\partial \rho }{\partial t}+u\frac{\partial \rho }{\partial x}+\rho \frac{\partial u}{\partial x}=0$

(Eq. 501)

Momentum equation:

$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+\frac{1}{\rho }\frac{\partial p}{\partial x}=0$

(Eq. 502)

Any thermodynamic property can be expressed as a function of two other thermodynamic properties. This means that we can get density as a function of pressure and entropy: $\rho =\rho (p,s)$ and therefore

$d\rho ={\left(\frac{\partial \rho }{\partial p}\right)}_{s}dp+{\left(\frac{\partial \rho }{\partial s}\right)}_{p}ds$

(Eq. 503)

Assuming isentropic flow $ds=0$ gives

$d\rho ={\left(\frac{\partial \rho }{\partial p}\right)}_{s}dp$

(Eq. 504)

$\begin{array}{cc}& \frac{\partial \rho }{\partial t}={\left(\frac{\partial \rho }{\partial p}\right)}_s\frac{\partial p}{\partial t}=\frac{1}{a^2}\frac{\partial p}{\partial t}\\ & \\ & \frac{\partial \rho }{\partial x}={\left(\frac{\partial \rho }{\partial p}\right)}_s\frac{\partial p}{\partial x}=\frac{1}{a^2}\frac{\partial p}{\partial x}\end{array}$

(Eq. 505)

Now, insert \ref{eq:rhotop} in \ref{eq:pde:cont} gives

$\frac{\partial p}{\partial t}+u\frac{\partial p}{\partial x}+\rho a^2\frac{\partial u}{\partial x}=0$

(Eq. 506)

Dividing \ref{eq:pde:cont:b} by $\rho a$ gives

$\frac{1}{\rho a}(\frac{\partial p}{\partial t}+u\frac{\partial p}{\partial x})+a\frac{\partial u}{\partial x}=0$

(Eq. 507)

A slightly modified form of the momentum equation is obtained by multiplying and dividing the last term by $a$

$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+\frac{1}{\rho a}\left(a\frac{\partial p}{\partial x}\right)=0$

(Eq. 508)

If the continuity equation on the form \ref{eq:pde:cont:c} is added to the momentum equation on the form \ref{eq:pde:mom:c}, we get

$[\frac{\partial u}{\partial t}+(u+a)\frac{\partial u}{\partial x}]+\frac{1}{\rho a}[\frac{\partial p}{\partial t}+(u+a)\frac{\partial p}{\partial x}]=0$

(Eq. 509)

If, instead, the continuity equation on the form \ref{eq:pde:cont:c} is subtracted from the momentum equation on the form \ref{eq:pde:mom:c}, we get

$[\frac{\partial u}{\partial t}+(u-a)\frac{\partial u}{\partial x}]+\frac{1}{\rho a}[\frac{\partial p}{\partial t}+(u-a)\frac{\partial p}{\partial x}]=0$

(Eq. 510)

Since $u=u(x,t)$, we have from the definition of a differential

$du=\frac{\partial u}{\partial t}dt+\frac{\partial u}{\partial x}dx=\frac{\partial u}{\partial t}dt+\frac{\partial u}{\partial x}\frac{dx}{dt}dt$

(Eq. 511)

Now, let $dx/dt=u+a$

$du=\frac{\partial u}{\partial t}dt+(u+a)\frac{\partial u}{\partial x}dt=[\frac{\partial u}{\partial t}+(u+a)\frac{\partial u}{\partial x}]dt$

(Eq. 512)

which is the change of $u$ in the direction $dx/dt=u+a$

In the same way

$dp=\frac{\partial p}{\partial t}dt+\frac{\partial p}{\partial x}dx=\frac{\partial p}{\partial t}dt+\frac{\partial p}{\partial x}\frac{dx}{dt}dt$

(Eq. 513)

and thus, in the direction $dx/dt=u+a$

$dp=\frac{\partial p}{\partial t}dt+(u+a)\frac{\partial p}{\partial x}dt=[\frac{\partial p}{\partial t}+(u+a)\frac{\partial p}{\partial x}]dt$

(Eq. 514)

If we go back and examine Eqn. \ref{eq:nonlin:a}, we see that Eqns. \ref{eq:du:b} and \ref{eq:dp:b} appear in the equation and thus it can now be rewritten as follows

$\frac{du}{dt}+\frac{1}{\rho a}\frac{dp}{dt}=0\Rightarrow du+\frac{dp}{\rho a}=0$

(Eq. 515)

Eqn. \ref{eq:nonlin:a:ode} applies along a $C^{+}$ characteristic, i.e., a line in the direction $dx/dt=u+a$ in $xt$-space and is called the compatibility equation along the $C^{+}$ characteristic. If we instead chose a $C^{-}$ characteristic, i.e., a line in the direction $dx/dt=u-a$ in $xt$-space, we get

$du=[\frac{\partial u}{\partial t}+(u-a)\frac{\partial u}{\partial x}]dt$

(Eq. 516)

$dp=[\frac{\partial p}{\partial t}+(u-a)\frac{\partial p}{\partial x}]dt$

(Eq. 517)

which can be identified as subsets of Eqn. \ref{eq:nonlin:b} and thus

$\frac{du}{dt}-\frac{1}{\rho a}\frac{dp}{dt}=0$

(Eq. 518)

In order to fulfil the relation above, either $du=dp=0$ or

$du-\frac{dp}{\rho a}=0$

(Eq. 519)

Eqn. \ref{eq:nonlin:b:ode} applies along a $C^{-}$ characteristic, i.e., a line in the direction $dx/dt=u-a$ in $xt$-space and is called the compatibility equation along the $C^{-}$ characteristic.

So, what we have done now is that we have have found paths through a point ($x_1$, $t_1$) along which the governing partial differential equations Eqns. \ref{eq:nonlin:a} and \ref{eq:nonlin:b} reduces to the ordinary differential equations \ref{eq:nonlin:a:ode} and \ref{eq:nonlin:b:ode}. The $C^{+}$ and $C^{-}$ characteristic lines are physically the paths of right- and left-running sound waves in the $xt$-plane.

Riemann Invariants

If the compatibility equations are integrated along respective characteristic line, i.e., integrate \ref{eq:nonlin:a:ode} along the $C^{+}$ characteristic and \ref{eq:nonlin:b:ode} along the $C^{-}$ characteristic, we get the Riemann invariants $J^{+}$ and $J^{-}$.

$J^{+}=u+\int \frac{dp}{\rho a}=const$

(Eq. 520)

$J^{-}=u-\int \frac{dp}{\rho a}=const$

(Eq. 521)

The Riemann invariants are constants along the associated characteristic line.

We have assumed isentropic flow and thus we may use the isentropic relations

$p=C_1{T}^{\gamma /(\gamma -1)}=C_2{a}^{2\gamma /(\gamma -1)}$

(Eq. 522)

where $C_1$ and $C_2$ are constants. Differentiating Eqn. \ref{eq:isentropic:a} gives

$dp=C_2\left(\frac{2\gamma }{\gamma -1}\right)a^{[2\gamma /(\gamma -1)-1]}da$

(Eq. 523)

Now, if we further assume the gas to be calorically perfect

$a^2=\gamma RT=\frac{\gamma p}{\rho }\Rightarrow \rho =\frac{\gamma p}{a^2}$

(Eq. 524)

Eqn. \ref{eq:isentropic:a} in \ref{eq:isentropic:c} gives

$\rho =C_2\gamma a^{[2\gamma /(\gamma -1)-2]}$

(Eq. 525)

and thus

$J^{+}=u+\int \frac{C_2\left(\frac{2\gamma }{\gamma -1}\right)a^{[2\gamma /(\gamma -1)-1]}}{C_2\gamma a^{[2\gamma /(\gamma -1)-2]}a}da=u+\left(\frac{2}{\gamma -1}\right)\int da$

(Eq. 526)

$J^{+}=u+\frac{2a}{\gamma -1}$

(Eq. 527)

$J^{-}=u-\frac{2a}{\gamma -1}$

(Eq. 528)

Eqns. \ref{eq:riemann:a:b} and \ref{eq:riemann:b:b} are the Riemann invariants for a calorically perfect gas. The Riemann invariants are constants along $C^{+}$ and $C^{-}$ characteristics and if the situation shown in Fig. \ref{fig:characteristics} appears, that fact can be used to calculate the flow velocity and speed of sound in the location ($x_1$, $t_1$).

$J^{+}+J^{-}=u+\frac{2a}{\gamma -1}+u-\frac{2a}{\gamma -1}=2u\Rightarrow u=\frac{1}{2}(J^{+}+J^{-})$

(Eq. 529)

$J^{+}=u+\frac{2a}{\gamma -1}=\frac{1}{2}(J^{+}+J^{-})+\frac{2a}{\gamma -1}\Rightarrow a=\frac{\gamma -1}{4}(J^{+}-J^{-})$

(Eq. 530)

Moving expansion waves

Moving Expansion Waves

The expansion wave propagation into the driver section in a shock tube can be described using characteristic lines.

The expansion is propagating into stagnant fluid in region four (the driver section), which means that the flow properties ahead of the expansion wave are constant.

$J_a^{+}=J_b^{+}$

(Eq. 531)

$J^{+}$ invariants constant along $C^{+}$ characteristics

$J_a^{+}=J_c^{+}=J_e^{+}$

(Eq. 532)

$J_b^{+}=J_d^{+}=J_f^{+}$

(Eq. 533)

Since $J_a^{+}=J_b^{+}$ this also implies $J_e^{+}=J_f^{+}$. In fact, since the flow properties ahead of the expansion are constant, all $C^{+}$ lines will have the same $J^{+}$ value.

$J^{-}$ invariants constant along $C^{-}$ characteristics

$J_c^{-}=J_d^{-}$

(Eq. 534)

$J_e^{-}=J_f^{-}$

(Eq. 535)

$\begin{array}{cc}& u_e=\frac{1}{2}(J_e^{+}+J_e^{-})\\ & u_f=\frac{1}{2}(J_f^{+}+J_f^{-})\\ & J_e^{-}=J_f^{-}\\ & J_e^{+}=J_f^{+}\end{array}\}\Rightarrow u_e=u_f\Rightarrow a_e=a_f$

(Eq. 536)

Due to the fact the $J^{+}$ is constant in the entire expansion region, $u$ and $a$ will be constant along each $C^{-}$ line.

The constant $J^{+}$ value can be used to obtain relations for the variation of flow properties through the expansion region. Evaluation of the $J^{+}$ invariant at any position within the expansion region should give the same value as in region 4.

$u+\frac{2a}{\gamma -1}=u_4+\frac{2a_4}{\gamma -1}=0+\frac{2a_4}{\gamma -1}$

(Eq. 537)

and thus

$\frac{a}{a_4}=1-\frac{\gamma -1}{2}\left(\frac{u}{a_4}\right)$

(Eq. 538)

Eqn. \ref{eq:expansion:a} and $a=\sqrt{\gamma RT}$ gives

$\frac{T}{T_4}={[1-\frac{\gamma -1}{2}\left(\frac{u}{a_4}\right)]}^2$

(Eq. 539)

Using isentropic relations, we can get pressure ratio and density ratio

$\frac{p}{p_4}={[1-\frac{\gamma -1}{2}\left(\frac{u}{a_4}\right)]}^{2\gamma /(\gamma -1)}$

(Eq. 540)

$\frac{\rho }{{\rho }_4}={[1-\frac{\gamma -1}{2}\left(\frac{u}{a_4}\right)]}^{2/(\gamma -1)}$

(Eq. 541)

Shock-tube relations

From the analysis of the incident shock, we have a relation for the induced flow behind the shock

$u_2=u_p=\frac{a_1}{{\gamma }_1}(\frac{p_2}{p_1}-1){\left(\frac{\left(\frac{2{\gamma }_1}{{\gamma }_1+1}\right)}{\left(\frac{{\gamma }_1-1}{{\gamma }_1+1}\right)+\left(\frac{p_2}{p_1}\right)}\right)}^{1/2}$

(Eq. 542)

The velocity in region 3 can be obtained from the expansion relations

$\frac{p_3}{p_4}={[1-\frac{{\gamma }_4-1}{2}\left(\frac{u_3}{a_4}\right)]}^{2{\gamma }_4/({\gamma }_4-1)}$

(Eq. 543)

Solving for $u_3$ gives

$u_3=\frac{2a_4}{{\gamma }_4-1}[1-{\left(\frac{p_3}{p_4}\right)}^{({\gamma }_4-1)/(2{\gamma }_4)}]$

(Eq. 544)

There is no change in pressure or velocity over the contact surface, which means $u_2=u_3$ and $p_2=p_3$.

$u_2=\frac{2a_4}{{\gamma }_4-1}[1-{\left(\frac{p_2}{p_4}\right)}^{({\gamma }_4-1)/(2{\gamma }_4)}]$

(Eq. 545)

Now, we have two ways of calculating $u_2$. Setting Eqn. \ref{eq:shocktube:up:a} equal to Eqn. \ref{eq:shocktube:up:d} leads to the shock tube relation

$\frac{p_4}{p_1}=\frac{p_2}{p_1}{\{1-\frac{({\gamma }_4-1)(a_1/a_4)(p_2/p_1-1)}{\sqrt{2{\gamma }_1[2{\gamma }_1+({\gamma }_1+1)(p_2/p_1-1)]}}\}}^{-2{\gamma }_4/({\gamma }_4-1)}$

(Eq. 546)