In this final post on compressible fluid flow, we will take a look at an application of the machinery developed over the previous three posts to the converging-diverging nozzle, often called a rocket nozzle or a de Laval nozzle.  The setup in straightforward even though it is probably unfamiliar.  We imagine a reservoir containing a gas at a high pressure and, not often, high temperature that serves to supply the flow through the nozzle.  By convention, the thermodynamic state variables in the reservoir are given a ‘0’ subscript and are called the stagnant conditions as there is no flow within the reservoir itself.  This reservoir feeds into a converging section of a nozzle in which the cross-sectional area decreases steadily until a minimum occurs at what is called the throat of the nozzle.  After the throat, the cross-sectional area increases again until the nozzle ends at the exit.  The nozzle and reservoir are placed into a larger pressure vessel (sometimes called the receiver) where the pressure outside of the exit, called by convention the back pressure, $P_b$.

Initially the two pressures are equal and there is no flow from the reservoir through the nozzle.  As the back pressure is dropped, flow develops from the reservoir.  de Laval’s innovation was to realize that if the flow can be brought to sonic conditions (i.e., Mach number $M=1$) at the throat then the diverging section would further speed up the flow to supersonic speeds (see the discussion at the end of the second post of this series).

The goal of this post is to discuss the flow profiles along the length of the nozzle, from inlet at the reservoir, to exit into the pressure vessel, as a function of the ratio $P_b/P_0$ and to give a flavor of how the previous equations are used to quantitatively.  Despite the relative simplicity of the system (reservoir, nozzle, and pressure vessel) a rich set of flow patterns result.  Qualitatively, these fall into two broad classes with several variations within them. 

In the first class (shown in blue traces in the figure below), the flow is everywhere isentropic from the reservoir. through the nozzle, and into the pressure vessel.  Included in this class is the flat trace where $P_b/P_0 = 1$ where no flow results.  As $P_b/P_0$ is lowered, flow commences with the speed of the flow increasing as it approaches the throat from the converging section and then subsequently slowing down in the diverging section.  The maximum speed in the nozzle (which for this class is always at the throat) continues to increase until as $P_b/P_0$ is lowered until a critical value of the back pressure, $P_{b_s}$ when the flow is sonic (i.e., $M=1$) in the throat but still subsonic everywhere else in the nozzle (including the diverging section).

At this point, lowering the $P_b/P_0$ ratio further develops the second class of flow (shown in red traces in the figure below) where some portion of the flow in the diverging section is supersonic.  A shock develops somewhere downstream of the throat in order to decelerate the flow back to subsonic conditions so that it can match boundary conditions with the environment in the pressure vessel.  The only remaining question is precisely where the shock is found.  Five possible options are available:  1) the shock is found in the nozzle between the throat and the exit, 2) the shock falls at the exit, 3) the shock is beyond the nozzle but the exit pressure $P_e < P_b$, 4) the shock is beyond the nozzle and $P_e = P_b$, and 5) the shock is beyond the nozzle but $P_e > P_b$.  Options 3), 4), and 5) are called over-expanded, perfectly-expanded, and under-expanded flows, respectively (see discussion of The Converging-Diverging Nozzle applet).

In the third post of this series, we derived expressions relating the pressure, density, and temperature anywhere in an isentropic portion of the flow to the stagnation properties in the reservoir.  The relations derived were:

\[ \frac{P_0}{P} = \left[1 + \frac{\gamma – 1}{2} M^2 \right]^{\frac{\gamma}{\gamma -1}} \; ,\]

\[ \frac{\rho_0}{\rho} = \left[1 + \frac{\gamma – 1}{2} M^2 \right]^{\frac{1}{\gamma -1}} \; ,\]

and

\[ \frac{T_0}{T} = 1 + \frac{\gamma – 1}{2} M^2 \; .\]

We also derived the area-Mach number relationship, (AMR),

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

that specified what the ratio of the local area, $A$, to the area of the throat, $A^*$, must be to achieve a target value of Mach number $M$.  A typical problem where these equations are useful comes from Fluid Mechanics Demystified by Merle Potter and reads

The solution starts with recognizing that if the flow in the diverging portion of the nozzle is to be supersonic then choked conditions occurs at the throat and, as a result, the throat area coincides with the critical area $A_{throat} = A^*$.  The area ratio of the exit to the throat is then $\frac{A}{A^*} = \left( \frac{24}{10} \right)^2=5.76$. Next we have to invert the AMR to find the value of $M$ that corresponds to this ratio.  A bisection, specifically adapted to the two possible dynamical behaviors: subsonic and supersonic, was written in python.  The resulting value was $M=3.3244$ was then plugged into the $P_0/P$ relation and the reciprocal taken to arrive at the exit pressure being $P_{exit} = 6746.659 \, Pa$, corresponding to an exit pressure only about $1.6%$ as large as the stagnant pressure.

The possibility of a shock in the possible flows requires us to add three new relations to our collection.  Our starting point for deriving the so-called jump conditions across the shock are the three basic fluid equations for mass, momentum, and energy, derived in the first post, but simplified by the tacit assumption that the shock is negligibly thick so that the cross-sectional area is the same on both sides of the shock, giving:

\[ \rho_u V_u = \rho_d V_d \; , \]

\[ (\rho_u {V_u}^2 + P_u) = (\rho_d {V_d}^2 + P_d) \; ,\]

and

\[ h_u + \frac{{V_u}^2}{2} = h_0 = h_d + \frac{{V_d}^2}{2} \; , \]

respectively.  Note that the subscripts `u’ and `d’ stand for upstream and downstream, respectively.

As before, our thermodynamics will assume the fluid to be an ideal gas with the equation of state given by $P/\rho = R T$.  The general form of enthalpy $h = e + P/\rho$ takes the simple form $h = (\gamma/\gamma-1) P/\rho$.  The state variables relate to their values at different parts of the flow via 

\[ \frac{P_0}{P} = \left( \frac{\rho_0}{\rho} \right)^{\gamma} = \left( \frac{h_0}{h} \right)^{\frac{\gamma}{\gamma-1}} \; . \]

The speed of sound in the gas, which is temperature dependent, takes on the forms

\[ c = \sqrt{\gamma R T } = \sqrt{\frac{\gamma P}{\rho}} = \sqrt{(\gamma-1) h} \; , \]

each useful in context.

Finally, the energy equation, which had been relegated to a by-stander role in incompressible fluid flow, offers two very useful relations.  The first, which will be termed the stagnant enthalpy expression (SEE), expresses the stagnant enthalpy in terms of the local sound speed via $(\gamma -1) h = c^2$ giving

\[ h_0 = h + \frac{V^2}{2} = h \left(1 + \frac{V^2}{2h} \right) = \frac{c^2}{\gamma-1} \left(1 + \frac{\gamma – 1}{2} M^2 \right) \; . \]

The second, which will be termed the local enthalpy expression (LEE), expresses the local enthalpy in terms of the local flow speed by

\[ (\gamma – 1) h = (\gamma – 1) \left( h_0 – \frac{1}{2}V^2 \right) \; . \]

One caution when using these relations.  The reservoir is operational defined as the location where the flow is isentropically brought to a stop.  Since a shock produces a great deal of entropy for any fluid element crossing from one side to another, the downstream thermodynamic conditions for a reservoir will be different from the actual reservoir attached to the nozzle and care must be taken not to confuse the notional one with the actual physical one.

To derive the normal shock relations, we follow these notes. 

Begin by dividing the momentum equation by the mass continuity equation to get

\[ V_u – V_d = \frac{1}{\gamma} \left( \frac{{c_d}^2}{V_d} – \frac{{c_u}^2}{V_u} \right) \; .  \]

Eliminate the local speeds of sound by relating them to their local enthalpies and then using the LEE to get

\[ V_u – V_d = \frac{\gamma -1}{\gamma} \left[ \frac{h_0 (V_u – V_d)}{V_u V_d} + \frac{1}{2}(V_u – V_d) \right] \; .\]

Dividing both sides by $V_u – V_d$ and rearranging yields

\[ \left( \frac{\gamma+1}{2} \right)^2 = \frac{h_0^2 (\gamma-1)^2}{{V_u}^2 {V_d}^2 } \; .\]

Next, employ a neat trick by writing the numerator on the left-hand side as ${h_0}^2 (\gamma-1)^2 = h_{0u} (\gamma-1) h_{0d} (\gamma – 1) $.  Use the SEE to express each of these expressions in terms of the local sound speed and Mach number.  Doing so yields

\[ \left( \frac{\gamma+1}{2} \right)^2 = \frac{1}{{M_u}^2} \left[1 + \frac{\gamma-1}{2} {M_u}^2 \right] \frac{1}{{M_d}^2} \left[1 + \frac{\gamma-1}{2} {M_d}^2 \right] \; . \]

Solving for $M_d$ is relatively painless with the substitutions $L = (\gamma-1)/2$ and $Q = (\gamma +1)/2$.  These auxiliary variables keep the clutter down and the observation that $L^2 – Q^2 = -\gamma$ simplifies things enormously so that we arrive at 

\[ {M_d}^2 = \frac{1 + \left(\frac{\gamma – 1}{2}\right) {M_u}^2 }{\gamma {M_u}^2 – \left( \frac{1-\gamma}{2} \right)} \; . \]

This relation allows u get the jump in all of the local thermodynamic variables.  Starting with the density, use the continuity equation to

\[ \frac{\rho_d}{\rho_u} = \frac{V_u}{V_d} = \frac{{V_u}^2}{V_u V_d} \; .\]

Using the definition of Mach number and the SEE together yields

\[ {V_u}^2 = \frac{(\gamma-1)h_0}{1 + \frac{\gamma-1}{2} {M_u}^2} {M_u}^2 \; . \]

Earlier in the derivation of the Mach jump relations we found that $1/V_u V_d = (\gamma+1)/2 h_0(\gamma-1)$.  Substituting these expressions back into the density equation yields

\[ \frac{\rho_d}{\rho_u} = \frac{(\gamma+1) {M_u}^2}{2 + (\gamma-1){M_u}^2} \; . \]

The momentum equation  gives

\[ P_d – P_u = \rho_u {V_u}^2 – \rho_d {V_d}^2 = \rho_u {V_u}^2 \left(1 – \frac{\rho_d V_d V_d}{\rho_u V_u V_u} \right) = \rho_u {V_u}^2 \left(1 – \frac{\rho_u}{\rho_d} \right)\; . \]

Substitute $\rho V^2 = \gamma P M^2$ from the definition of the speed of sound and dividing by $P_u$ gives

\[ \frac{P_d}{P_u} – 1 = \gamma {M_u}^2 \left( 1- \frac{\rho_u}{\rho_d} \right) \; . \]

Using the already obtained expression for the density ratio (and noting to flip it based on which density is in the numerator) followed with some simplifications gives

\[ \frac{P_d}{P_u} = 1 + \frac{2\gamma}{\gamma+1} \left( {M_u}^2 – 1 \right) \; . \]

Finally, when needed, the temperature ratio comes subtituting into the ideal gas equation of state the ratios for density and pressure

\[ \frac{T_d}{T_u} = \frac{P_d}{P_u} \frac{\rho_u}{\rho_d} = \left[ 1 + \frac{2\gamma}{\gamma+1} \left( {M_u}^2 – 1 \right) \right] \frac{2 + (\gamma-1){M_u}^2}{(\gamma+1) {M_u}^2} \; .\]

We are now in a position to solve another typical problem associated with converging-diverging nozzle design, again taken from Merle Potter, that reads

Since the shock is at the exit, we can assume that the upstream side of the shock is supersonic at the area of the exit.  So, the value of $A/A* = 4$.  Using the inverse AMR function, the upstream Mach number is $M_u = 2.940$.  Then using the standard isentropic relations the pressure ratio is $P_u/P_0 = 0.0298$ resulting in an upstream pressure of $P_u = 11,914.79 Pa$.  Finally, the Mach shock jump relationship yields $P_d/P_u = 9.919$ for a downstream pressure, which is equal to the receiver pressure, of $P_d = 118179.93 Pa$. 

Of course, this analysis just barely scratches the surface but these basic relations demonstrate the interplay between three basic laws of fluid mechanics: continuity, momentum, and energy and how they contribute to the practical construction of a converging-diverging nozzle.  And, so, I guess this is rocket science.

In this installment we derive a central result in compressible fluid flow showing the inner working of the converging-diverging nozzle.  The key relations are the jump conditions are the mass conservation equation

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

and the energy equation

\[ h_1 + \frac{V_1^2}{2} = h_2 + \frac{V_2^2}{2} \; .\]

As a reminder, the enthalpy is defined in terms of the specific internal energy $e$ by $h = e + P/\rho$.  Correspondingly, the first law of thermodynamics (as usually expressed in terms of the specific internal energy)

\[ de = T ds + \frac{P}{\rho^2} d\rho \]

becomes

\[ dh = T ds + \frac{dP}{\rho} \; .\]

The final set of equations come from the specific fluid that is flowing.  For simplicity, we will deal with compressible gas flow and will take the gas to be ideal.  The equation of state for an ideal gas is

\[ \frac{P}{\rho} = RT \]

with the usual relations between the heat capacity at constant pressure $c_P$ and the heat capacity at constant volume $c_V$ of

\[ c_P – c_V = R \;  \]

and

\[ c_P/c_V = \gamma \; , \]

where $\gamma$ is a constant usually taken to have the value of 1.4 for ordinary air.

For an ideal gas, the specific enthalpy (like the specific energy) is strictly a function of temperature given by

\[ dh = c_P dT \; , \]

where we will often assume that $c_P$ is constant over a wide range of conditions; an assumption that allows us to integrate to $h = c_P T$ immediately.  In other words, we assume the gas is a calorically perfect ideal gas – this is actually a reasonably good assumption given that $c_P$ for nitrogen varies only about 1.5

The speed of sound in the ideal gas is also only a function of temperature given by

\[ c = \sqrt{ \gamma R T } \; .\]

The final assumption is that the flow is isentropic, which is a good assumption since the flow happens quickly so that a negligible amount of heat flows into the system.  Under this assumption, the equation of state can be substituted into the first law to yield

\[ c_P dT = \frac{dP}{\rho} \; .\]

The next step is to eliminating the density by way of the ideal gas law to get an equation in terms of pressure and temperature

\[ c_P dT = \frac{dP}{P/RT} \; . \]

Re-arranging gives a simple set of differential equations

\[ \frac{dT}{T} = \frac{R}{c_P} \frac{dP}{P} \; ,\]

which integrates to

\[ \frac{T_1}{T_2} = \left( \frac{P_1}{P_2} \right)^{\frac{\gamma-1}{\gamma}} \; ,\]

where $R = c_P – c_V$ and $c_P/c_V = \gamma$.

When the analogous steps are followed with $P$ eliminated in favor of $\rho$, one arrives at

\[ \frac{T_1}{T_2} = \left( \frac{\rho_1}{\rho_2} \right)^{\gamma – 1} \; . \]

The value of this latter equation is that it allows the elimination of temperature (simply equation both forms of $T_1/T_2$ and the emergence of the polytropic equation

\[ \frac{P_1}{P_2} = \left( \frac{\rho_1}{\rho_2} \right)^{\gamma} \; \]

relating pressure and density for an isentropic flow.

With these assumptions and corresponding relations in hand, we now use them to relate the conditions from the inlet of the gas from a reservoir (e.g. a fuel tank) to the outlet of the nozzle.  The rocket engine will be the prototype for this situation.  We need to be able to relate the thermodynamics at each point in the flow to the thermodynamics in the reservoir.  Since the gas in the reservoir is not moving with respect to the nozzle we introduce the concept of stagnation, in which the flow is envisioned to have come to a stop and the only energy present being stored as internal energy.  The third jump condition relating enthalpy and bulk flow velocity at different points will be important.

Since there is a qualitative change between subsonic and supersonic flow (see last post’s discussion of how the area is related to Mach number), we will need to stitch together the two flows (stagnant-to-sonic with sonic-to-supersonic) at the critical point where $M = 1$.  Variables associated with this point are called critical and are decorated with a star (e.g. $A^*, P^*$ and so on).  In addition, we will decorate all variables associated with the reservoir where the flow is stagnant with a 0 subscript (e.g. $A_0, P_0$ and so on).

Let’s start on the subsonic side by relating the thermodynamics in the tank to the corresponding variables in the flow as it moves in the converging portion of the nozzle.  Since the flow is stagnant in the tank, we find that the stagnation enthalpy $h_0$ relates to the flow anywhere else in the nozzle by

\[ h_0 = h_1 + \frac{V_1^2}{2} \; . \]

Since the gas is calorically perfect, the enthalpy can be eliminated in favor of the temperatures to arrive at

\[ c_P T_0 = c_P T_1 + \frac{V_1^2}{2} \; .\]

Dividing by $c_P T_1$ and using the fact that $T_1 = {c_1}^2/\gamma R$ gives

\[ \frac{T_0}{T_1} = 1 + \frac{V_1^2}{2 c_P T_1} = 1 + \frac{V_1^2}{2 c_P ({c_1}^2/\gamma R)} \; .\]

In this last relation, the ratio of the flow velocity to the speed of sound is the Mach number.  Dropping the ‘1’ subscript and using the relations between $c_P$, $c_V$, and $R$, we arrive at what we will call the stagnation temperature relation

\[ \frac{T_0}{T} = 1 + \frac{(\gamma-1)}{2} M^2 \;. \]

Using the relationship between the temperature and pressure derived above gives the corresponding stagnation pressure equation

\[ \frac{P_0}{P} = \left[ 1 + \frac{(\gamma-1)}{2} M^2 \right]^{\frac{\gamma}{\gamma – 1}} \; \]

Return to the mass flow equation and, first eliminate the density using the equation of state and multiply by $c/c$, then express the sound speed in terms of the temperature to get

\[ {\dot m} = \rho A V = \frac{P}{R T} A \frac{V}{c} c = P A M \frac{c}{RT} = P A M \sqrt{ \frac{\gamma}{R T} } \; .\]

Next, solve the stagnation pressure equation for $P$ 

\[ P = P_0 \left( \frac{T_0}{T} \right)^{\frac{\gamma}{1 – \gamma}} \]

and then substitute into the mass flow equation

\[ {\dot m} = P_0 \left( \frac{T_0}{T} \right)^{\frac{\gamma}{1 – \gamma}} A M \sqrt{ \frac{\gamma}{R T}} \left( \frac{T_0}{T_0} \right)^{1/2} \\ = P_0 \sqrt{\frac{\gamma}{R T_0}} \left(\frac{T_0}{T}\right)^{\frac{\gamma}{1-\gamma}+\frac{1}{2}} A M \; .\]

Simplifying and using the stagnation temperature equation gives

\[ {\dot m} = P_0 \sqrt{\frac{\gamma}{RT_0}} A M \left[ 1 + \frac{\gamma-1}{2} M^2 \right]^{\frac{\gamma+1}{2(1-\gamma)}} \; .\]

This combination is constant no matter where in the flow it is evaluated.  Equate the expression evaluated at the critical point to the expression anywhere else in the flow to get

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

Eliminating the terms common to both sides and rearranging we can express the cross-sectional area at any point in the nozzle to the point where the flow goes critical as

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

This is the well-known relation that forms the basis for designing and analyzing the nozzle flow.  It will be the basis for our final post on compressible flow next month.

This month’s column, the second in a series of three on compressible flow, derives the key relationships for understanding steady, isentropic flow in a converging-diverging nozzle.  The starting point of these derivations are the three isentropic jump conditions derived in the last post:

\[ \begin{array}{ll} {\textrm{continuity:}}& {\rho}_1 V_1 A_1 = {\rho}_2 V_2 A_2 \equiv {\dot m}\\{\textrm{momentum:}}& (\rho_1 V_1^2 + P_1) A_1 = (\rho_2 V_2^2 + P_2) A_2 \\{\textrm{energy:}}& h_1 + \frac{{V_1}^2}{2} = h_2 + \frac{{V_2}^2}{2} \end{array} \; , \]

along with the first law of thermodynamics (cast in intensive variables by dividing both sides of $dE = TdS – P dV$ by mass $m = \rho V$)

\[ de = T ds + \frac{P}{\rho^2} d \rho \; ,\]

and the equation of state for an ideal gas (also in intensive variables)

\[ \frac{P}{\rho} = R T \; . \]

The $T ds$ term can be dropped since isentropic flow means that the specific entropy remains constant.  Note that the specific internal energy and specific enthalpy ($h = e + P/\rho$) are simply functions of temperature given by $de = c_{{\mathcal V}} dT$ and $dh = c_P dT$, respectively, where $c_{{\mathcal V}}$ and $c_P$ are the heat capacities at constant volume and constant pressure and are taken to be constants.  Other thermodynamic relations will be taken as given since proof will take us too far afield.

There are three ingredients in understanding the steady, isentropic flow through a converging-diverging nozzle (also known as a de Laval nozzle): 1) calculating the Mach number, 2) relating the cross-sectional area of the nozzle ($A$) to the Mach number, and 3) relating the thermodynamics variables of temperature ($T$), the pressure ($P$) and the density ($\rho$) to the Mach number.  This post will deal with items 1 and 2 with item 3 being covered in the next post.

This post and the next are a synthesis and clarification of the ideas and materials found in Chapter 9 of Merle Potter’s Fluid Mechanics Demystified and the 5-part series of YouTube lectures on compressible fluid flow from the University of Florida (part 1, part 2, part 3, part 4, and part 5)

Calculating the Mach Number

The central parameter in the theory is the Mach number, which is innocently defined as the ratio of the flow speed to the speed of sound.  What makes this definition so deceptively simple looking is we are used to thinking that the speed of sound is constant since the temperature range that we normally encounter is such that the variations in the speed of sound is negligible.  However, interesting compressible fluid flows often have very large variations in temperature over short durations or small distances and so Mach number varies depending on the location of the flow within the nozzle.

To calculate the speed of sound, which we will denote as $c$, we take the first differential of the continuity and momentum equations and since we are interested in the local speed over a column of fluid with constant cross section we can ignore variations in the area (imagine being outside).  The resulting equations are:

\[ d \rho c + \rho d c = 0 \; \]

and

\[ d\rho c^2 + 2 c \rho dc + dP = 0 \; . \]

Solving the first equation for $\rho dc = – c d \rho$, substituting this result into the second, and solving for $c^2$ gives

\[ c^2 = \frac{d P}{d \rho} \; \]

An ideal gas obeys a polytropic equation of state

\[ P = \rho^{\gamma} \; , \]

where the adiabatic index $\gamma$ is defined as

\[ \gamma = \frac{c_P}{c_{{\mathcal V}}} \; .\]

As a result,

\[ \frac{dP}{d\rho} = \gamma \rho^{\gamma – 1} = \gamma \frac{\rho^{\gamma-1}}{\rho} = \gamma \frac{P}{\rho} \; .\]

Substituting this result back into the first equation and taking the square root gives

\[ c = \sqrt{ \frac{\gamma P }{\rho} } = \sqrt{\gamma R T} \; , \]

where the ideal gas law was used in the last step.

The Mach number is the ratio of the flow speed to the speed of sound and is given by

\[ M = \frac{V}{c} \; . \]

Relating Cross-Sectional Area and the Mach Number

Most of us have had the experience with running water where we restrict the cross-sectional area of the outlet and watch the fluid flow increase (a thumb over the end of a garden hose being the usual scenario).  The reason for this speed increase is that the pressure difference between the inlet and the outlet sets up a given mass flow rate.  When the area drops, the flow speed has to increase to maintain that flow rate.

This approach works for flows in most everyday experiences but fails for compressible flow, explaining why a simple converging nozzle is insufficient to force a fluid into supersonic flow.

To see why this is we will derive a relationship between the cross-sectional area and the Mach number.

Start by taking the first differential of the third jump condition

\[ d \left( h + \frac{V^2}{2} \right) = dh + V dV = 0 \; .\]

Next, eliminate the differential of the enthalpy using the first law ($dh = dP/\rho$) to get

\[ \frac{dP}{\rho} + V dV = 0 \; .\]

Multiply the first term by $d\rho/d\rho$ to get

\[ \frac{dP d\rho}{\rho d\rho} + VdV = 0 \; .\]

Recognizing that $c^2 = dP/d\rho$ and dividing both sides by $V^2$ gives

\[ \frac{c^2}{V^2} \frac{d\rho}{\rho} + \frac{dV}{V} = 0 \;. \]

Rearranging this equation gives

\[ \frac{d\rho}{\rho} = -M^2 \frac{dV}{V} \; . \]

Next take the differential of the first jump condition and divide by $\rho V A$ to get

\[ \frac{d\rho}{\rho} + \frac{dV}{V} + \frac{dA}{A} = 0 \;. \]

Finally substitute into this equation the previous one and re-arrange to arrive at

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

This equation contains a wealth of information:

When $M > 1$ the flow slows down (speeds up) when the area decreases (increases)
When $M < 1$ the flow speeds up (slows down) when the area decreases (increases)
When $M = 1$ the flow remains at the same speed even as the area changes

These observations explain why a simple converging nozzle can’t force the flow to supersonic speeds. Once the pressure difference between the inlet and outlet is large enough that the flow speed equals the speed of sound ($M=1$), the flow is choked and can get no faster.  Practitioners usually speak of this situation is that choked flow can’t communicate a further pressure drop at the outlet upstream so that the flow can respond.

It is also worth pointing out something that is easy to overlook in the above manipulation.  While the steps may seem haphazard and uninspired there is a method to the madness. There is a hierarchy of thermodynamic relations that are combined: the first law of thermodynamics is a universal law; the third jump condition represents energy conservation for a general fluid flow; and the equation of state for an ideal gas is a specific law for a very special kind of fluid.  The first law links a fluid’s energy conservation with the speed of sound for an ideal fluid, linking flow velocity to Mach number and the compressibility of the fluid (as represented by $d\rho$).  The final step uses another universal law – mass continuity – to arrive at a relationship between a fractional change in the flow speed to a fractional change in volume (as represented by $dA$).

Next post puts these pieces together with a set of relations between the fluid’s thermodynamic variables and the Mach number to explain the function of a converging-diverging nozzle and how it can push fluid flows past the speed of sound.

This month’s column briefly examines compressible fluid flow.  This journey will be a short excursion into a rather complicated discipline with only two aims.  The first is to give a flavor of how and why the energy equation, which has been essentially unused in previous analyses, plays a role and how additional thermodynamic information (spatial and temporal variations in density and temperature) comes into play.  The second is to setup all the equations needed to understand the basic physics associated with that object that ushered in the space age: the rocket nozzle, which will be the subject of next month’s entry.

Following the usual introductory approaches (e.g. Eric R. Pardyjak’s online notes and Chapter 9 of Merle C. Potter’s Fluid Mechanics Demystified), we will look at only simplified geometries, such as the converging nozzle shown here

with a flow from left-to-right as indicated.

This analysis assumes that all of the flow variables are constant over each cross-section but can differ from cross-section to cross-section.  The first condition means that the flow is effectively inviscid, since the boundary layer is ignorable for the vast majority of the flow, and it is steady.  All of these simplifying assumptions result in algebraic flow equations that come from the integral forms of the fluid equations rather than the differential equations.  In addition, we will ignore body forces (i.e. gravity) and only allow pressure-driven flows.

The conservation (i.e. Eulerian) form of the equations is slightly easier to work with and these will be our starting point for the various simplifications that we impose (although with one extra wrinkle in the case of energy).  All of the derivations will involve some vector quantity (say the velocity $${\vec V}$$) being integrated over a surface after being projected onto the local surface normal.  This process results in an average value of the normal component of the vector over the entire face.  This quantity will be denoted simply as a scalar (e.g. $$V_1$$) with appropriate labeling indicating which segment of the bounding surface.

Continuity Equation

Our starting equation is

\[ \int_{{\mathcal Vol}} d {\mathcal Vol} \, \frac{\partial \rho}{\partial t} + \int_{{\mathcal S}} d {\mathcal S} \, \rho {\vec V} \cdot {\hat n} = 0 \; .\]

Since the flow is steady, the first term vanishes.  If we imagine the surface as shown in the nozzle above, then only the flow through the endcaps ($${\mathcal S}_1$$ and $${\mathcal S}_2$$) is non-zero and

\[ \int_{{\mathcal Vol}} d {\mathcal Vol} \, {\vec V} = {\rho}_1 {\hat n}_1 \cdot {\vec V}_1 A_1 + {\rho}_2 {\hat n}_2 \cdot {\vec V}_2 A_2 = 0 \; .\]

Accounting for the direction of the flow relative to the respective normals we arrive at

\[ {\rho}_1 V_1 A_1 = {\rho}_2 V_2 A_2 \equiv {\dot m} \; . \]

Momentum Equation

Next, we consider the momentum equation expressed as

\[ \int_{{\mathcal Vol}} d {\mathcal Vol} \, \frac{\partial \left(\rho {\vec V}\right)}{\partial t} + \int_{{\mathcal S}} d {\mathcal S} \, (\rho {\vec V}) {\vec V} \cdot {\hat n} = \int_{{\mathcal S}} d {\mathcal S} \, {\mathbf T} \cdot {\hat n} \; ,\]

where $${\mathbf T}$$ is the Cauchy stress tensor.  Following the same approach as above, and recalling that the body forces are zero (i.e. $${\mathbf T} = – P {\mathbf I}$$), the Cauchy momentum equation takes the form:

\[ {\rho}_1 {\vec V}_1 ( {\vec V}_1 \cdot {\hat n}_1) A_1 + {\rho}_2 {\vec V}_2 ( {\vec V}_2 \cdot {\hat n}_2) A_2 = -P_1 A_1 {\hat n}_1 – P_2 A_2 {\hat n}_2 \; .\]

Again taking into account the direction of the flow relative to the normals, we arrive at

\[ {\rho}_2 V_2 A_2 {\vec V}_2 – {\rho}_1 V_1 A_1 {\vec V}_1 = ( P_1 A_1 – P_2 A_2 ) {\hat n}_2 \; ,\]

where we’ve used $${\hat n}_1 = – {\hat n}_2$$.  Two things are worth noting about this equation.  First, the direction of the flow must be along (either parallel or anti-parallel) to the end-caps, nicely in keeping with the assumption that the flow varies only from cross-section to cross-section and not within a cross-section.  Second, the continuity equation can be used to eliminate many of the terms in favor of the constant mass flow rate $${\dot m}$$.  Keeping all this in mind brings us to the scalar equation

\[ {\dot m} (V_2 – V_1) = (P_1 A_1 – P_2 A_2) \; . \]

Energy Equation

Up to this point, we’ve not had to occasion to use the energy equation for the simple fact that in incompressible flows without explicit heat flows and generation, the internal energy of each fluid element remains constant; only the kinetic energy changed.  Indeed, the energy equation for such a flow follows immediately by taking the dot product of the Cauchy momentum equation with the velocity.  In these simple cases, the density is a given and only pressure and velocity are unknown.  The four equations of continuity and 3 components of the Cauchy momentum equation sufficed to solve for the four unknowns.

Once we allow the flow to be compressible, the density becomes an unknown and the also the internal energy enters into the picture.  The energy equation provides one more equation.  The final equation comes from a thermodynamic equation of state.  For this latter equation, we will assume the ideal gas law.  Often, it is more convenient to switch from internal energy to temperature by using the ideal gas relation

\[ e = c_v T \]

and

Fourier’s law relating heat flow to temperature via

\[ {\vec q} = – k \nabla T \; .\]

There is an additional thermodynamic consideration (the one wrinkle mentioned at the beginning of this post).  It is the introduction of the enthalpy $$h = e + P/\rho$$ that combines the internal energy with the pressure work.  Its introduction is really a matter of convenience rather than an important insight but, as the majority of the discipline’s practioners uses it, it would be a mistake to omit.

To see how the energy equation admits its inclusion, start with the energy equation in its full glory

\[ \int_{{\mathcal Vol}} d {\mathcal Vol} \, \frac{\partial}{\partial t} \left[ \rho \left(e + \frac{V^2}{2} \right) \right] + \int_{{\mathcal S}} d {\mathcal S} \, \rho \left(e + \frac{V^2}{2} \right) {\vec V} \cdot {\hat n} = \\ \int_{{\mathcal Vol}} d {\mathcal Vol} \, \rho \left( {\vec a}_{body} \cdot {\vec V} + {\dot Q} \right) + \int_{{\mathcal S}} d {\mathcal S} \, \left( {\mathbf T} \cdot {\vec V} + {\vec q} \right) \cdot {\hat n} \; .\]

By the above assumptions, the volumetric heating and the heat flow are zero, the body forces are zero, the stress tensor simply accounts for the pressure, and the flow is steady.  Under these conditions, the equation reduces to

\[ \int_{{\mathcal S}} d {\mathcal S} \, \rho \left(e + \frac{V^2}{2} \right) {\vec V} \cdot {\hat n} = – \int_{{\mathcal S}} d {\mathcal S} \, P {\vec V} \cdot {\hat n} \; .\]

Collecting terms gives

\[ \int_{{\mathcal S}} d {\mathcal S} \, \rho \left( e + \frac{P}{\rho} + \frac{V^2}{2} \right) {\vec V} \cdot {\hat n} = 0 \; , \]

or, when integrated, accounting for the normal on the endcaps, and using the definition of enthalpy

\[ \rho_2 V_2 A_2 \left( h_2 + \frac{{V_2}^2}{2} \right) – \rho_1 V_1 A_1 \left( h_1 + \frac{{V_1}^2}{2} \right)= 0 \; .\]

Summary

Finally, we can summarize all the relations

\[ \begin{array}{ll} {\textrm{continuity:}}& {\rho}_1 V_1 A_1 = {\rho}_2 V_2 A_2 \equiv {\dot m}\\{\textrm{momentum:}}& (\rho_1 V_1^2 + P_1) A_1 = (\rho_2 V_2^2 + P_2) A_2 \\{\textrm{energy:}}& h_1 + \frac{{V_1}^2}{2} = h_2 + \frac{{V_2}^2}{2} \end{array} \; , \]

which are basically the Rankine-Hugonoit jump conditions used in shock wave propagation analysis (assuming $$A_1 = A_2$$).

Next month we will apply these equations to elementary compressible fluid and the production of shocks.

This month’s column finishes up the study of elementary viscous flow and presents the final two examples, both dealing with shear-induced flow.  This type of flow (called Couette flow) typically results from the no-slip condition in which a moving boundary drags the fluid that touches it into motion.  That layer of the fluid then drags the fluid layer next to it and so on.

For both of these examples, we are going to ignore body forces and external pressure differences and focus on the boundary motion alone.  These assumptions eliminate several terms from the incompressible Navier-Stokes equations leaving us with

\[ \frac{\partial {\vec V}}{\partial t} + ( {\vec V} \cdot \nabla) {\vec V} = \nu \nabla^2 {\vec V} \; , \]

augmented with the usual incompressibility assumption $$\nabla \cdot {\vec V} = 0 $$.

Cylindrical Couette Flow

In this example, we imagine a fluid trapped between two cylindrical shells, an inner and an outer shell, of radii $$R_{inner} \equiv R_{1}$$ and $$R_{outer} \equiv R_{2}$$, respectively.

For sake of generality, both shells will be allowed to rotate about the vertical axis (out of the page) with angular velocities $$\Omega_{inner} \equiv \Omega_{1}$$ and $$\Omega_{outer} \equiv \Omega_{2}$$.  We seek to describe the velocity profile of the fluid in steady motion.  With the radii and angular velocities assumed to be known there will be a fairly simple relationship between the velocity profile and the kinematic viscosity $$\nu$$.  A viscometer uses this experimental setup to measure kinematic viscosity by tracking the torque required to keep the inner shell moving at a constant rate (the out shell is usually fixed).

The Navier-Stokes equations are best presented in cylindrical coordinates since these coordinates match the geometry involved.  The added benefit is that expressing the derivatives in cylindrical coordinates is instructive in its own right.

The key point to focus on is that in cylindrical coordinates, the basis vectors in directions perpendicular to the $$z$$-axis change as a function of azimuthal angle $$\theta$$.  The radial unit vector changes as

\[ \partial_{\theta} {\hat e}_r = {\hat e}_{\theta} \; \]

and the azimuthal unit vector as

\[ \partial_{\theta} {\hat e}_{\theta} = – {\hat e}_r \; .\]

In steady flow, by symmetry, the velocity field is a function of radius alone

\[ {\vec V} = V_{\theta}(r) {\hat e}_{\theta} \; \].

The nabla operator takes the well-known form

\[ \nabla = \frac{1}{r} \partial_{r} ( r ) + \frac{1}{r} \partial_{\theta} + \partial_z \; .\]

Under the flow geometry assumptions, the convective portion of the material derivative becomes

\[ ({\vec V} \cdot \nabla) = \frac{V_{\theta}}{r} \partial_{\theta} \; .\]

Applying this operator to the velocity field yields

\[ ({\vec V} \cdot \nabla) {\vec V} = \left( \frac{V_{\theta}}{r} \partial_{\theta}\right) (V_{\theta} {\hat e}_{\theta}) = – \frac{V_{\theta}^2}{r} {\hat e}_{r} \; \]

The Laplacian of the vector field becomes

\[ \nabla^2 {\vec V} = \left( \partial_{r}^2 + \frac{1}{r} \partial_{r} + \frac{1}{r^2} \partial_{\theta}^2 \right) ( V_{\theta} {\hat e}_{\theta} ) = \left( \partial_{r}^2 V_{\theta}+ \frac{1}{r} \partial_{r} V_{\theta} – \frac{V_{\theta}}{r^2} \right) {\hat e}_{\theta} \; . \]

The assumption of steady flow then leaves the incompressible Navier-Stokes equations, in component form, as

\[ – \frac{V_{\theta}^2}{r} = \frac{1}{\rho} \partial_{r} P \; , \]

\[ 0 = – \frac{1}{\rho r} \partial_{\theta} P + \nu \left( \partial_{r}^2 V_{\theta} + \frac{1}{r} \partial_{r} V_{\theta} – \frac{V_{\theta}}{r} \right) \; , \]

\[ 0 = – \frac{1}{\rho} \partial_{z} P \; . \]

Note that the incompressibility condition is automatically satisfied by the form of the velocity field.

From the last equation the pressure cannot be a function of $$z$$ and since everything else in the second equation is only a function of $$r$$ so must $$\partial_{\theta} P = f(r)$$.  Integrating gives

\[ P(r) = f(r) \theta + g(r) \; , \]

the form of which requires $$f(r) = 0$$ otherwise the pressure would be a multi-valued function.  The second equation now becomes

\[ 0 = \nu \left( \partial_{r}^2 V_{\theta} + \frac{1}{r} \partial_{r} V_{\theta} – \frac{V_{\theta}}{r} \right) \; , \]

The general solution is

\[ V_{\theta} = A r + \frac{B}{r} \; .\]

In the interest of brevity, the solutions for the constants are

\[ A = \frac{ \Omega_{2} r_{2}^2 – \Omega_{1} r_{1}^2 }{r_{2}^2 – r_{1}^2} \; \]

and

\[ B = \frac{ (\Omega_{2} – \Omega{1}) r_{1}^2 r_{2}^2 }{ r_{2}^2 – r_{1}^2} \; , \]

which are easily obtained by applying the boundary conditions to the general relation $$V_{\theta} = r \Omega$$.

Once $$V_{\theta}$$ is known, the first equation gives the radial pressure gradient due to the rotation.

Linear Couette Flow

In this last example, we will solve the time varying Couette flow setup up between two plates separated by a distance $$w$$ and extending far in the $$x$$- and $$z$$-directions so that we can ignore the edge effects and model the flow as two-dimensional.

We imagine that the fluid is at rest prior to $$t=0$$ at which time the upper plate is set into motion with velocity $${\vec V} = U {\hat e}_x$$.  There is a well-known steady solution in which the velocity linearly varies from zero on the lower plate to $$U$$ on the upper plate but it is interesting to see how the transient flow settles into this steady state.

By our earlier, assumptions, the flow is of the form (automatically satisfying the incompressibility condition)

\[ {\vec V} = V_x(t,y) {\hat e}_x \; , \]

which has to satisfy the following form of the Navier-Stokes equations

\[ \rho \partial_{t} V_x = \mu \partial_{y}^2 V_x \; . \]

The boundary conditions and initial condition are

\[ \textrm{BCs:} \;\;\;\;\;\;\;\; \left\{ \begin{array}{l} V_x(t,0) = 0 \\ V_x (t,w) = U \end{array} \right. \; \]

and

\[ \textrm{IC:} \;\;\;\;\;\;\;\; V_x(0,y) = 0 \;\;\;\;\; 0 \leq y \leq w \; , \]

respectively.

The solution strategy is to employ separation of variables.  To that end, following Acheson, write the velocity as a sum of the steady solution and a time-varying piece where the latter has homogeneous boundary conditions

\[ V_x(t,y) = \frac{U}{w} y + {\tilde V}_x(t,y) = \frac{U}{w} y + f(t)g(y) \; . \]

After performing the separation of variables, we arrive at

\[ \frac{1}{\nu} \frac{1}{f} \partial_{t} f = – s_{n}^2 = \frac{1}{g} \partial_{y}^2 g \; , \]

where $$s_n$$ is a separation constant.

The spatial equation, $$g^{\prime\prime} + s_{n}^2 g = 0$$ has general solutions $$g(y) = A \sin (s_n y) + B \cos (s_n y) $$.  Applying the boundary conditions requires $$B=0$$ and $$s_n = n \pi / w \;\; n = 1, 2, 3, \ldots$$.  The time equation is easily integrated to $$f(t) = \exp( -s_n\nu t )$$.  The general solution is

\[ V_x(t,y) = \frac{U}{w} y + \sum_{n=1}^{\infty} A_n \exp \left( – \frac{ n^2 \pi^2}{w^2} \nu t \right) \sin \left( \frac{n \pi}{w} y \right) \; . \]

The initial conditions are satisfied by picking the $$A_n$$ appropriately.  The key to doing this is to exploit the orthogonality

\[ \int_0^w \sin \left( \frac{m \pi}{w} y \right) \sin \left( \frac{n \pi}{w} y \right) = \frac{w}{2} \delta_{mn} \; \]

to get

\[ A_m = -\frac{2}{w} \int_0^w \frac{Uy}{w} \sin \left( \frac{m \pi}{w} y \right) dy = \frac{ 2 U (-1)^n }{n \pi} \; . \]

The following plot shows this solution for $$w = 5mm$$ and $$\nu = 0.8926 mm^2/s$$.

Note the majority of the velocity is confined near the top plate at short times (this is essentially the boundary layer) but how quickly the solution tends to the steady state.

As a way of seeing this boundary layer as a feature of the transient behavior consider what happens if the top plate moves in an oscillatory fashion with velocity $$U \cos (\omega t) $$.  Assuming that the fluid velocity has the form $$ V_x(t,y) = f(y) e^{i \omega t} $$, the Navier-Stokes equations yield

\[ f” – i \omega/\nu f = 0 \; .\]

Following the usual strategy $$f = exp(\lambda y)$$ gives values

\[ \lambda_{\pm} = \pm e^{i \pi/4} \sqrt{\frac{\omega}{\nu}} \; .\]

Some spirited complex algebra leads to

\[ V_x(t,y) = U e^{-ky} e^{i (\omega t – k y) } \; ,\]

where $$k = \sqrt{\omega/2 \nu}$$.  As the following animation shows, the velocity at the upper plate follows the motion of the plate but the main action is confined to a thin layer.  Essentially, the motion is always transient and the influence of the fluid motion near the moving plate isn’t able to diffuse into the bulk of the fluid.

This month’s column deals with viscous flow in a variety of situations.  Three main physical effects can drive elementary fluid flow: the pull of gravity coaxing a fluid downhill; a pressure differential, such as in a pipe, causing the fluid to flow from inlet to outlet; or shearing motion, which, in the usual scenario, involves a moving boundary dragging the fluid.  We’ll be looking at four examples covering different aspects of these three effects.  The first two will illustrate the effect of gravitational body forces and pressure differentials to drive the fluid and will be covered in this column.  The next two examples deal with steady and unsteady shear-induced flows and will be featured next month.

But, before we deal with specific examples, we will need to settle on a model of the fluid.  Usually, when seeking analytic solutions describing the motion of non-ideal Newtonian fluids (already a simplification compared to the amazingly rich and diverse types of fluids that nature offers), one either chooses to assume the flow is viscous and incompressible or that it is inviscid and compressible.  For these examples, we will choose the former.

The starting point will be the incompressible Navier-Stokes equations, of which the flux-conservative form is

\[ \frac{\partial \rho {\vec V}}{\partial t} + \nabla \cdot (\rho {\vec V} {\vec V}) = – \nabla P + \mu \nabla^2 {\vec V} + \frac{1}{3} \mu \nabla \left( \nabla \cdot {\vec V} \right) + \rho {\vec g} \; .\]

Applying the incompressibility assumption ($$\nabla \cdot {\vec V} = 0$$) is easy to do using the Lagrangian form of the equations, but it is more instructive to play with the flux-conservative form.  Most of the action takes place on the left-hand side so let’s focus on that by expanding the terms to get

\[ \rho \frac{\partial {\vec V}}{\partial t} + {\vec V} \frac{\partial \rho}{\partial t} + {\vec V} ( {\vec V} \cdot \nabla \rho ) + \rho ({\vec V} \cdot \nabla) {\vec V} + {\vec V} \rho (\nabla \cdot {\vec V}) \; .\]

The last term vanishes under the incompressibility assumption.  The second and third terms, combined, are

\[ {\vec V} \left( \frac{\partial \rho}{\partial t} + ( {\vec V} \cdot \nabla ) \rho \right) \; .\]

The quantity in the parentheses vanishes when applying incompressibility to the continuity equation, as follows:

\[ \frac{\partial \rho}{\partial t} + ({\vec V} \cdot \nabla ) \rho + \rho \nabla \cdot {\vec V} = \frac{\partial \rho}{\partial t} + ({\vec V} \cdot \nabla ) \rho = 0 \; ,\]

where the first expression is the full continuity equation and the second is what is left after the incompressibility assumption is applied.  The left-hand side is now simply

\[\rho \left( \frac{\partial {\vec V}}{\partial t} + ({\vec V} \cdot \nabla) {\vec V} \right) \;. \]

The third term on the right-hand side vanishes due to the incompressibility assumption.  Dividing by the density gives the general form of the Navier-Stokes equations for incompressible fluid flow:

\[ \frac{\partial {\vec V}}{\partial t} + ({\vec V} \cdot \nabla) {\vec V} = – \frac{1}{\rho} \nabla P + \nu \nabla^2 {\vec V} + {\vec g} \; , \]

where $$\nu = \mu/\rho$$ is the kinematic viscosity.

This equation is what we will use in the following examples.

Gravity-driven Flow

This example derives from the very nice YouTube lecture by Rick Sellens.  Here we imagine flow of a fluid, such as oil, down an inclined plane (flow in the $$x$$ direction), as shown in the figure below.  The height of the oil is larger at the top of the ramp but it quickly thins as gravity accelerates the oil.  Eventually, the viscous force equals the component of the gravitational force parallel to the plane and, thereafter, the height of the oil is constant (the flow is said to be fully developed) and the velocity profile depends only on $$y$$.  Once the flow establishes itself, this height profile remains constant in time and the flow is steady, allowing us to set all time derivatives to zero.  We will assume that the incline is wide enough in the $$z$$-direction (not pictured) that we can approximate the flow as being two-dimensional.  Next, we will confine our description to the fully developed portion of the flow ($$x > 0$$).  Finally, we will assume a no-slip boundary condition at the incline plane and, since the flow is open to the air on top and air is not particularly viscous, we will assume the shear stress on the upper surface to be zero.

Applying the steady, fully-developed, and 2-D flow assumptions leads to the continuity equation

\[ \partial_x V_x = 0 \; \]

and the $$x$$- and $$y$$-component Navier-Stokes equations of

\[ 0 = g_x + \nu \frac{\partial^2}{\partial_y^2} V_X \; \]

and

\[ 0 = -\partial_y P + \rho g_y \; ,\]

where $$g_x = g \sin(\theta) $$ and $$g_y = g \cos(\theta)$$, are the components of the gravitational acceleration parallel and perpendicular to the surface of the incline, respectively.

The continuity equation merely confirms what we already assumed, $$V_x = V_x(y)$$.  The $$y$$ component of Navier-Stokes says that the pressure at the surface of the incline is higher than at the top due to the weight of the oil above.  All the action is in the $$x$$-component equation.

Integrating the equation once yields

\[ \frac{d}{dy} V_x = – \frac{g_x}{\nu} y + C_1 \; . \]

The constant of integration is determined by the shear-stress-free boundary condition at the free surface and gives $$C_1 = g_x h/\nu$$.

Integrating the equation a second time yields

\[ V_x = \frac{g_x}{\nu} \left(h y – \frac{1}{2} y^2 \right) + C_2 \; .\]

The constant of integration is identically zero after the application of the no-slip boundary condition at the incline ($$V_x(y=0)=0$$).

Pressure-driven Flow

This example draws from the Faculty of Kahn video on Poiseuille’s law.  In this scenario, we imagine a horizontal pipe with an inlet pressure, $$P_{in}$$, and an outlet pressure, $$P_{out}$$, with $$P_{in} > P_{out}$$ and that the flow has had time to become steady.  Cylindrical geometry clearly suggests itself and we expect that the only flow is in the $$z$$-direction ($$V_r = V_{\theta} = 0$$) and its only variation is radial ($$V_z = V_z(r)$$).

As discussed in an earlier post, usually one must take care to differentiate the basis vectors as well as the components when setting up fluid equations in curvilinear coordinates.  However, the above assumptions are sufficiently limiting that only the $$z$$-component equation survives in the form

\[ \frac{d P}{d z} = \nu \frac{1}{r} \frac{d}{dr} \left( r \frac{d V_z}{dr} \right) \; . \]

The general solution to this is

\[ V_z(r) = \frac{P_{in} – P_{out}}{4 \mu L} r^2 + C_1 \ln r + C_2 \; .\]

The constants of integration follow from a boundary condition that requires $$V_z(r=0)$$ to be finite and the usual no-slip boundary condition at the walls.  The final form is then

\[ V_z(r) = \frac{P_{in} – P_{out}}{4 \mu L}(r^2 – R^2) \; .\]

Integrating this form over the surface area of the pipe gives the famous Poisueille flow rate relationship that says, for a given flow rate, the pressure differential is inversely proportional to the radius to the fourth power.  This relationship has profound consequences in the field of medicine.

This month’s column begins a new chapter in the ongoing study of fluid dynamics.  Up to this point, the general equations of fluid motion have been derived from Newton’s laws using the transformation between a fluid-following viewpoint (Lagrangian point-of-view) and a field concept (Eulerian point-of-view).  General considerations of mass conservation and the principles of momentum and energy led to a system of non-linear partial differential equations.  Ignoring viscosity allowed the equations to simplify to Euler’s equations for an ideal fluid and several specific examples drove home some of the kinematics, dynamics, and methods of solution.

Unfortunately, the situation changes quite dramatically when viscosity is brought back into consideration and the full Navier-Stokes equations (which are themselves something of a simplification) begin to stare back from the page.  What is the general mode of attack, and how does one think about the solutions that result?  The enormous difficulty of finding solutions becomes even more acute when two facts are driven home.  First, it is extremely difficult to find solutions for the simpler set of equations for an ideal fluid, and the solutions shown to date (e.g. the Rankine Vortex), complicated as they may seem, are special cases passed down from one generation to the next precisely because they had a solution.  Second, general analytic solutions remain elusive nearly 200 years after Navier first derived these equations in 1822 despite millions of dollars offered to those who could derive such solutions.

To set the stage for an exploration that will consume many columns to come, this installment will present some of the observations and follow many of the arguments gleaned from the article Ludwig Prandtl’s Boundary Layer by John D. Anderson Jr published in the December 2005 issue of Physics Today.

According to Anderson, Ludwig Prandtl

revolutionized fluid dynamics with a very simple to understand and yet highly non-intuitive concept of the boundary layer that exists in a fluid near a rigid surface.

To understand the importance of Prandtl’s result one has to start the story in the latter half of the 18^th century, around the time Euler derived the ideal fluid equations bearing his name.  The French mathematician Jean le Rond d’Alembert derived a paradox associated with ideal fluid flow.  He proved that, for an incompressible fluid with no viscosity (i.e., an ideal fluid), the drag force on an object placed within this flow, moving with constant relative velocity with respect to the fluid, is identically zero, an obviously incorrect result.

Acheson presents a mathematical derivation of d’Alembert’s paradox In Section 4.13 of his book Elementary Fluid Dynamics.  This derivation examines the force and momentum balance in the fluid and, using Bernoulli’s streamline theorem, concludes that, for the fluid to be moving at a uniform speed far upstream and far downstream from the body, the drag force on the body must vanish.  This result contradicts everyday observations about the motion of solid objects in air or water.

Clearly something was wrong with the theory of ideal fluid, but the specific root of the problem escaped identification for roughly 150 years until Prandtl introduced his boundary layer at a conference in 1904.  The modern picture of the boundary layer is best understood in terms of the motion of air over a cambered wing as shown in the figure below (based on Figure 2 of Anderson).

Away from the wing, the ideal fluid model works well in describing the streamlines.  However, close to the wing’s surface, viscous effects become dominant in a thin layer whose base attaches to the wing’s surface.  This no-slip condition demands that the fluid be at rest at the surface and that the fluid velocity rapidly rise from zero (where the layer attaches to the wing) to the ‘free stream’ velocity of the bulk inviscid flow outside the layer (shown in the inset).  In this way, d’Alembert’s paradox is conceptually resolved.  The drag force is non-zero, as observed, and the bulk flow is inviscid or ideal (depending on whether one allows the viscosity-free flow to be compressible or not), as observed, but the mathematical conditions by which d’Alembert derived the paradox are not true everywhere – the thin boundary layer makes all the difference.

Prandtl also predicted that, in certain cases, this layer can separate from the surface, which can cause additional drag.  To quote Prandtl

Indeed, modern theory classifies the boundary layer behavior as falling into two broad classes.  Laminar motion describes the flow when the layer stays fixed to the surface.  Turbulent motion describes the flow when the layer separates and produces vortices in wake.  Laminar flow produces less drag than turbulent flow but it is unstable and thus harder to produce or control.  Turbulence changes the characteristics of the flow far from the surface causing additional drag.

In addition to these two critical physical observations, Anderson points out that the concept of the boundary layer also enabled the first quantitative calculations of aerodynamic drag.  The reason for this is that the boundary layer equations exhibit an entirely different type of mathematical behaviors than the equations describing the bulk flow outside.  Although both sets of equations are coupled, nonlinear, partial differential equations, the Navier-Stokes equations used to describe the bulk flow external to the layer are elliptic (at least when the flow is subsonic) while the simplification of these equations in the boundary layer are parabolic.  This distinction allows the solution methods used to describe the bulk motion to ‘decoupled’ from the methods used to describe the drag and since the methods for solving parabolic equations are much easier to deal with than those for elliptic equations real progress can be made.

What is particularly amazing about Prandtl’s contribution is that the original presentation in 1904 lasted no longer than 10 minutes.  I would say that the return on investment resulting from that rather short talk have been nothing short of phenomenal.

One of the more difficult things for any student of fluid dynamics to understand is the distinction between streamlines and streaklines and pathlines.  It ranks up there with mastering the difference between the Eulerian and Lagrangian viewpoints and, in fact, the same difficulty underpins both – the ability to tell the difference between an attribute of the field (Eulerian point-of-view) and an attribute of the particle or fluid element (Lagrangian point-of-view).  And, I’ll confess, that even though I understand the definitions, it isn’t always obvious how to apply them in different contexts.

Wikipedia has a somewhat useful article on the distinction between streamlines, streaklines, and pathlines but the math presented obscures the underlying physics and the foundation (either directly or by reference) is sorely lacking.

The first thing to note is that if the flow is steady, then all three lines result in the same trace in the fluid.  This is because each fluid element, flowing through a given point, follows exactly the same trajectory as the fluid element that precedes it or follows it.  The distinction between these lines (and possibly a few more differently defined traces) comes when the flow is unsteady (varies in time).

The formal definitions of these lines is as follows:

• Streamline – a trace formed so that it is everywhere tangent to the velocity field of some local region of the fluid at a given time
• Streakline – a trace of the front of all fluid elements, at a given time, that originated from the same physical place (for example a nozzle that flows smoke into a wind tunnel)
• Pathline – a trace of the motion through time of a particular fluid element

These formal definitions don’t do a lot to promote understanding the situation.  So, we turn to an introductory diagram to explain the difference between streakline and pathline (streamlines will be discussed below) and then we will look at a simple example of unsteady flow to see the mathematical difference of all three.

Consider fluid starting to come out of a garden hose as the faucet is turned on.  The flow may start weakly and then climb until it levels off at the steady state value.  Below is a diagram of such a situation (inspired by the lecture from the Cowan Academy).  The individual pathlines are quite different and the time the fluid elements came through the outlet of the hose are different.  But their motion is stopped at a common time and the line joining them is the streakline.  It is particularly useful in that it represents an experiment that can be performed and measured.

Streamlines are best shown using a quiver plot.  To be concrete consider the velocity field in 2 dimensions to be given by

\[ V_x= \alpha x \; \]

and

\[ V_y = -\alpha y \; ,\]

where $$\alpha$$ is some real number.  A plot of the flow velocity at a grid of points is called a quiver plot and the results for the flow above are (with $$\alpha = 1$$)

Note that the color of the vectors match their length with red corresponding the largest flow rate and blue the smallest, that they suggest that hyperbolic pattern with flow coming in along the y-axis and flowing along the x-axis, and the is a stagnation point at the origin where the fluid elements there are at rest.

Acheson, in Chapter 1 of Elementary Fluid Dynamics, has a nice problem associated with this flow that accentuates the difference between Lagrangian and Eulerian viewpoints.  In the Eulerian viewpoint, the concentration of a pollutant is given by

\[ c(x,y,t) = \beta x^2 y e^{-\alpha t} \; ,\]

which clearly varies in space and time and fades in strength at a given point as time progresses.  The question Acheson asks is does the concentration in a fluid element change.  To answer, we employ the material derivative, since it follows the flow of a given material element.  Taking the material derivative of the concentration gives

\[ \frac{D}{Dt} c(x,y,z) = \frac{\partial c(x,y,t)}{\partial t} + \left( \vec V \cdot \nabla \right) c(x,y,t) \; .\]

The right-hand side of the equation is a set of partial derivatives that, when expanded, is

\[ \partial_t c + \alpha x \partial_x c – \alpha y \partial_y c = -\alpha \beta x^2 y e^{-\alpha t} + 2 \alpha \beta x^2 y e^{-\alpha t} – \alpha \beta x^2 y e^{-\alpha t} = 0 \; . \]

So, the concentration of a fluid element is constant.  The reason the concentration decays in strength at a given point is simply that the flow shunts the more polluted elements away from any point and towards infinity.  Alternatively, one can use Lagrangian coordinates to describe the flow.  The position of a given fluid element is given by

\[ x = X e^{\alpha t} \; \; y = Y e^{-\alpha t} \; , \]

where $$X$$ and $$Y$$ are the x- and y-position of the fluid element at $$t=0$$.  In terms of these coordinates the concentration of a fluid element is

\[ c = \beta X^2 e^{-2 \alpha t} Y e^{\alpha t} e^{\alpha t} = \beta X^2 Y \; ,\]

which is manifestly time independent.

With this warm up exercise under our belt, let’s look at an unsteady flow where the difference between streamlines and pathlines becomes obvious.

Consider the flow

\[ V_x = V_{x_0} \; \; V_y = k t \; \]

The quiver plots for two different times are shown here (with $$V_{x_0} = 1$$ and $$k=1$$).

In both cases, the streamlines are straight lines but with a time varying slope, as evidenced by the different slopes at different times.

The pathlines come from integrating

\[ \frac{\partial x}{\partial t} = V_{x_0} \;\; \textrm{and} \;\; \frac{\partial y}{\partial t} = k t \; .\]

Performing the integration and and setting the integration constants ($$x(t=0) = X$$ and $$y(t=0) = Y $$) to zero gives

\[ x(t) = V_{x_0} t \;\; \textrm{and} \;\; y(t) = \frac{1}{2} k t^2 \; .\]

Solving the first equation for $$t$$ and substituting into the second equation gives the geometric path of

\[ y = \frac{1}{2} \frac{k}{V_{x_0}^2} x^2 \; , \]

which shows that the particles follow a parabolic pathline even though the streamlines are always linear.

This month’s column closes out the study of ideal fluids. The focus will be on two particular results that have only been touched upon briefly in other posts: 1) why the energy equation plays no role in solving flows for a given configuration and 2) what can happen when the density is not assumed to be constant. The latter point will only be touched upon in terms of how the relaxation of this constraint modifies the vorticity. The aim of this column is to provide a better understanding of how the physics of vortex generation and flow changes as we move from ideal fluids to viscous fluids. As a result, the analysis in this post and two preceding ones should serve as a baseline against the results for similar situations that arise from the Navier-Stokes equations. This post is strongly influenced by David Acheson’s Elementary Fluid Mechanics.

To start, let’s look at the momentum equation. In order to make this post self-contained, details will be repeated from other posts. In non-conservation form, the momentum equation is

\[ \rho \frac{D {\vec u}}{D t} = -\nabla P – \rho \vec g \; , \]

which is easy to remember based since it looks like Newton’s law. However, to really understand the content, we need to expand the material derivative to get

\[ \rho \left( \partial_t \vec u + (\vec u \cdot \nabla) \vec u \right) = -\nabla P – \rho \vec g \; .\]

The next step is to use the vector identity

\[ ( \vec u \cdot \nabla ) \vec u = \nabla \times (\nabla \times \vec u) – \frac{1}{2} \nabla( u^2 ) ;, \]

and to identify the vorticity as $$\vec \omega = \nabla \times \vec u$$.

Substituting these results back into the expanded momentum equation, gives

\[ \rho \left( \partial_t \vec u + \nabla \times \vec \omega \right) = -\nabla P – \rho \vec g – \frac{1}{2} \nabla (u^2) \; .\]

Now we take the dot product of the above equation with $$\vec u$$ to get

\[ \rho \partial_t \left( \frac{d}{dt} \frac{1}{2} u^2 \right) = – \vec u \cdot \nabla (P + gz) – \frac{1}{2} \vec u \cdot \nabla (u^2) \; .\]

Now we employ the vector identity

\[ \vec u \cdot \nabla H = H \nabla \cdot \vec u – \nabla \cdot (H \vec u) \; , \]

where $$H$$ is shorthand defined by

\[ H \equiv P + \rho z g + \frac{1}{2} u^2 \; .\]
Substituting this result in gives, mindful of the incompressibility of the fluid, $$\nabla \cdot \vec u$$, gives

\[ \partial_t \left( \frac{1}{2} \rho u^2 \right) = -\nabla \cdot \left[ H \vec u \right] \; . \]

In some sense, we are done, since the left-hand side has the power (time derivative of the kinetic energy) and the right-hand side has the flow of energy out of the fluid, since $$H$$ contains the energy density in the form of pressure (units of energy per unit volume), the gravitational potential energy density, and the kinetic energy density.

But a clear picture emerges if we employ the divergence theorem, which bring us to

\[ \frac{d}{dt} \int d {\mathcal Vol} \frac{1}{2} \rho u^2 = – \int d {\mathcal A} H \vec u \cdot {\hat n} \; , \]

which shows that the time rate of change of the kinetic energy of a fluid parcel is equal to the flow of the energy out from that parcel. Note that the energy here doesn’t involve any local thermodynamics and relies strictly on the bulk effects. That is why the solution to ideal fluid flow can ignore the energy equation and simply use momentum equation (as was done in the last two posts).

As we move to the Navier-Stokes equations, the presence of viscosity means that energy and momentum can flow between the smallest of fluid elements because viscosity implies that two fluid elements can interact in a ‘shearing’ fashion. Said differently, viscosity is a type of friction that excites internal degrees of freedom (as does temperature – but more on that in future posts).

Now, we look at the relaxing the assumption that an ideal fluid has constant density. Our starting point is Euler’s equation above, repeated here for convenience

\[ \partial_t \vec u + \vec \omega \times \vec u = -\frac{1}{\rho} \nabla \left( p + E_{mech} \right) \; ,\]

where $$E_{mech} = \Xi + 1/2 u^2$$ and $$\Xi$$ is the gravitational potential and the mass conservation equation

\[ \frac{D\rho}{Dt} + \rho \nabla \cdot {\vec u} = 0 \; . \]

Now take the curl of Euler’s equation to get
\[ \partial_t \vec \omega + \nabla \times \left( \vec \omega \times \vec u \right) = – \frac{1}{\rho} \nabla \times \vec F \; , \]

where $$\vec F = \nabla \left(p + E_{mech} \right)$$, for convenience.

Expand $$\nabla \left( \vec \omega \times \vec u \right)$$ as

\[ \nabla \times \left( \vec \omega \times \vec u \right) = (\vec u \cdot \nabla) \vec \omega – (\vec \omega \cdot \nabla ) \vec u + \omega (\nabla \cdot u) – \vec u (\nabla \cdot \vec \omega) \; .\]

If the fluid were incompressible $$\nabla \cdot \vec u = 0$$, but under these assumptions we need to keep this term, using the mass conservation to eliminate the divergence of the velocity in terms of the material derivative of the density. The one term that is identically zero is $$\vec u (\nabla \cdot \vec \omega)$$ since the gradient of a curl is zero.

The next ingredient is to use the identity

\[ \nabla \times \left( \frac{1}{\rho} \vec F \right) = \vec F \times \nabla \left( \frac{1}{\rho} \right) + \frac{1}{\rho} \nabla \times \vec F \]

Putting all these pieces together

\[ \frac{D}{Dt} \left( \frac{\vec \omega}{\rho} \right) = \left( \frac{\omega}{\rho} \cdot \nabla \right) \vec u – \frac{1}{\rho} \nabla \left(\frac{1}{\rho} \right) \times \nabla P \; \]

This equation generalizes the vorticity equation to compressible fluid flow. Note that if the pressure is a function of the density only, the gradient of the pressure is parallel to the gradient of the density and the second term on the right-hand vanishes.

This month’s column is going to look at a simple model of a vortex called the Rankine vortex.  The model is extraordinarily simple and it is somewhat astonishing how well it matches results from real world observations.  Giaiotti and Stel provide vivid comparisons between the Doppler observations of tornado’s radial velocity distribution and the agreement is strong.   This, despite the fact, that the Rankine vortex is really an ideal fluid model that doesn’t carry over into the realm of the Navier-Stokes equations where viscosity is important.

The presentation here follows David Acheson’s textbook Elementary Fluid Dynamics with some inspiration and direction from a lecture by Luca Brandt on the Rankine Vortex available here.

The Rankine vortex consists of two parts: an inner core that moves as a rigid body and a vortex free exterior.  Mathematically, the Rankine vortex is specified as

\[ R = \left\{ \begin{array}{c} \Omega r \;\; r < a \\ \frac{\Omega a^2}{r} \; \; r >a \end{array} \right. \, \]

with $$u_r = u_z = 0$$.

By its geometry, the Rankine vortex demands us to work in cylindrical coordinates.  There are two steps that one must pay attention to in order to setup the problem.  First, by standard vector calculus results, the convective derivative is given by

\[ {\vec u} \cdot \nabla = u_r \partial_r + \frac{u_{\theta}}{r} \partial_\theta + u_z \partial_z \; .\]

Second, the fluid velocity, given in terms of the radial, azimuthal, and vertical unit vectors, is

\[ {\vec u} = u_r {\hat e}_r + u_\theta {\hat e}_\theta + u_z {\hat e}_z \; ,\]

with the realization that the radial and azimuthal unit vectors change as a function of space according to

\[ \partial_\theta {\hat e_r} = {\hat e}_\theta \]

and

\[ \partial_\theta {\hat e_\theta} = -{\hat e}_r \; .\].

Now applying the convective derivative to the velocity,

\[ (\vec u \cdot \nabla) \left[ u_r {\hat e}_r + u_\theta {\hat e}_\theta + u_z {\hat e}_z \right] \; ,\]

formally, gives a lot of terms

\[ (\vec u \cdot \nabla) {\vec u} = \left[ u_r \partial_r u_r + \frac{u_\theta}{r} \partial_\theta u_r + u_z \partial_z u_r – \frac{u_\theta^2}{r} \right] {\hat e_r} \\ + \left[ u_r \partial_r u_\theta + \frac{u_r u_\theta}{r} + \frac{u_\theta}{r} \partial_\theta u_\theta + u_z \partial_z u_\theta \right] {\hat e}_\theta \\ + \left[ u_r \partial_r u_z + \frac{u_\theta}{r} \partial_\theta u_z + u_z \partial_z u_z \right] {\hat e}_z \; ,\]

but, by the assumptions in the Rankine vortex, most are zero.  Putting all these pieces together results in the convective derivative taking the form

\[ \left( \vec u \cdot \nabla \right) \vec u = -\frac{u_\theta^2}{r} {\hat e}_r \; . \]

Given that the flow is steady, the resulting model is the following partial, coupled, differential equations

\[ -\rho \frac{u_\theta^2}{r} {\hat e}_r = -\frac{\partial P}{\partial r} {\hat e}_r \; , \]

\[ 0 = – \frac{\partial P}{\partial \theta} {\hat e}_\theta \; , \]

and

\[ 0 = \left( -\frac{\partial P}{\partial z} – \rho g \right) {\hat e}_z \; . \]

The $$\theta$$ equation merely tells us that $$P = P(r,z)$$.  The $$z$$ equation is easily solved for both the inner and outer regions of the vortex because there is no discontinuity.  This gives

\[ P(r,z) = – \rho g z + f(r) \; .\]

The radial profile requires a bit more work as the inner and outer regions have to have separate integrations.  For the inner

\[ \frac{1}{\rho} \frac{\partial P}{\partial r} = \frac{ \Omega^2 r^2}{r} = \Omega^2 r \;\; r < a \; . \] In similar fashion, the outer integration yields \[ \frac{1}{\rho} \frac{\partial P}{\partial r} = \frac{\frac{\Omega^2 a^r}{r^2}}{r} = \frac{\Omega^2 a^4}{r^3} \;\; r > a \; . \]

Solving both of these equations gives

\[ P(r,z) = \left\{ \begin{array}{c} \frac{1}{2} \rho \Omega^2 r^2 – \rho g z + C_1 \;\; r < a \\ – \frac{1}{2} \rho \frac{\Omega^2 a^4}{r^2} – \rho g z + C_2 \;\; r > a \end{array} \right. \]

Continuity at $$r=a$$ gives

\[ -\frac{1}{2} \rho \Omega^2 a^2 – \rho g z + C_2 = \frac{1}{2} \rho \Omega^2 a^2 – \rho g z + C_1 \; , \]

which leads to

\[ C_2 – C_1 = \rho \Omega^2 a^2 \; .\]

We can determine the pressure difference between the core and far from the core by comparing the interrogating the above result.  The pressure at the core is

\[ P(0,z) = -\rho g z + C_1 \; \]

and the pressure at infinity is

\[ P(\infty,z) = – \rho g z + C_2 \; .\]

Their difference is

\[ P(\infty,z) – P(0,z) = C_2 – C_1 = \rho \Omega^2 a^2 \; , \]

which explains why the pressure at the center of a tornado (vortex) is lower than the pressure far from the core.

The free surface (the shape) of the vortex as a function of radius is found by setting $$P(r,z) = P_0$$ with the height being assumed to be $$z_0$$ at $$r=0$$.