The primary text is P.K. Kundu and I.M. Cohen, Fluid Mechanics, 3^rd edition (Academic Press, 2004), which covers all topics in the undergraduate portion. G.K. Batchelor, An Introduction to Fluid Dynamics (Cambridge University Press, 1967), remains the authoritative mathematical reference and is cited throughout for depth and rigour. For boundary layer theory and the physical intuition behind high-Reynolds-number flows, D.J. Acheson, Elementary Fluid Dynamics (Oxford, 1990), provides an accessible companion. The chapters on geophysical fluid dynamics draw on J. Pedlosky, Geophysical Fluid Dynamics, 2^nd ed. (Springer, 1987). The two graduate extension chapters at the end — on nonlinear waves (AMATH 867) and hydrodynamic stability (AMATH 863) — draw on G.B. Whitham, Linear and Nonlinear Waves (Wiley, 1974), and P.G. Drazin and W.H. Reid, Hydrodynamic Stability, 2^nd ed. (Cambridge, 2004), respectively. Open educational resources consulted include MIT OpenCourseWare 2.20 (Marine Hydrodynamics) and 12.800 (Fluid Dynamics of the Atmosphere and Ocean).

---

Part I — Core Theory of Fluid Mechanics

Chapter 1: Equations of Motion

1.1 The Continuum Hypothesis and the Fluid Element

Before writing down any equations, we must justify the very concept of a fluid described by smooth fields. Matter is discrete at the molecular scale — a cubic centimetre of water contains roughly \(10^{22}\) molecules. Yet we speak of a smooth velocity field \(\mathbf{u}(\mathbf{x},t)\) and a smooth pressure \(p(\mathbf{x},t)\) as if fluid properties varied continuously in space. The justification is the continuum hypothesis: we consider length scales much larger than the mean free path \(\ell_{mfp}\) but much smaller than the macroscopic scale \(L\) of the flow. For water at sea level, \(\ell_{mfp} \sim 3 \times 10^{-10}\) m — so the continuum approximation is valid down to nanometre scales for liquids and to roughly a micrometre for gases at sea level.

The Knudsen number \(Kn = \ell_{mfp}/L\) quantifies the validity of the continuum approximation: fluid mechanics holds when \(Kn \ll 1\). Satellite re-entry vehicles, microelectromechanical systems (MEMS), and high-altitude flight all push against this limit — but for the overwhelming majority of fluid flows encountered in engineering and geophysics, \(Kn\) is negligibly small and the continuum description is exact to any practical precision.

1.2 The Reynolds Transport Theorem

The fundamental conservation laws are most naturally stated for a material volume \(\mathcal{V}(t)\) — a volume that moves with the fluid, always enclosing the same fluid particles. To derive Eulerian (fixed-frame) equations from these integral laws, we need to differentiate an integral over a moving domain.

The first term is the rate of change of the integral at fixed volume; the second accounts for the volume sweeping up additional fluid as its boundary moves with velocity \(\mathbf{u}\). Applying the divergence theorem to the surface integral converts this to the volume integral of \(\nabla \cdot (f\mathbf{u})\). Since the result must hold for every material volume, the integrand itself must vanish pointwise (in smooth flows), yielding the Eulerian form of the conservation law.

1.3 Conservation of Mass

Mass is conserved: the total mass in a material volume does not change. The Reynolds transport theorem with \(f = \rho\) gives:

\[ \frac{\partial\rho}{\partial t} + \nabla \cdot (\rho\mathbf{u}) = 0 \]

Expanding the divergence using the product rule:

\[ \frac{D\rho}{Dt} + \rho(\nabla \cdot \mathbf{u}) = 0 \]

For an incompressible flow, the density of each fluid parcel is constant along its trajectory: \(D\rho/Dt = 0\). The continuity equation then reduces to the incompressibility constraint:

\[ \nabla \cdot \mathbf{u} = 0 \]

1.4 Conservation of Momentum — The Navier–Stokes Equations

Newton’s second law for a material volume, combined with Cauchy’s theorem for the stress tensor \(\tau_{ij} = -p\delta_{ij} + 2\mu\dot{e}_{ij}\) (for an incompressible Newtonian fluid), yields after the Reynolds transport theorem:

In component form with kinematic viscosity \(\nu = \mu/\rho\):

\[ \frac{\partial u_i}{\partial t} + u_j\frac{\partial u_i}{\partial x_j} = -\frac{1}{\rho}\frac{\partial p}{\partial x_i} + \nu\frac{\partial^2 u_i}{\partial x_j \partial x_j} + g_i \]

The Navier–Stokes equations are the complete description of an incompressible Newtonian fluid. Every topic in this course — potential flow, water waves, boundary layers, turbulence, geophysical dynamics — is a different limit or specialisation of these equations. The richness of fluid mechanics arises entirely from the nonlinear advection term \((u_j \partial u_i / \partial x_j)\) interacting with the linear viscous term and the pressure gradient.

1.5 Boundary Conditions

Boundary conditions are as important as the equations themselves. At a rigid solid wall:

• No-slip: the fluid velocity equals the wall velocity, \(\mathbf{u}_{fluid} = \mathbf{u}_{wall}\) (from molecular sticking)
• For inviscid flow (Euler equations), only the normal component is specified: \(\mathbf{u} \cdot \hat{n} = \mathbf{u}_{wall} \cdot \hat{n}\)

At a free surface (air–water interface):

• The normal velocity is continuous: no fluid passes through the surface
• The tangential stress is continuous (the air exerts a very small shear on the water surface)
• The pressure jump across the surface equals the surface tension times the curvature: \(\Delta p = \sigma \kappa\)

At infinity for an exterior flow problem:

• The velocity approaches the undisturbed far-field: \(\mathbf{u} \to \mathbf{U}_\infty\) as \(|\mathbf{x}| \to \infty\)

---

---

Chapter 2: Vorticity Dynamics

Vorticity is the fundamental measure of local rotation in a fluid. Defined as the curl of the velocity field, \(\boldsymbol{\omega} = \nabla \times \mathbf{u}\), it encodes the spin that each infinitesimal fluid element carries as it is swept through the flow. Understanding how vorticity is generated, transported, stretched, and destroyed is arguably the deepest single thread running through all of classical fluid mechanics: potential flow lives in the absence of vorticity, turbulence is dominated by its generation and cascade, geophysical flows are controlled by its planetary-scale budget, and stability theory asks when ordered vorticity fields break into disorder.

This chapter develops the dynamics of vorticity from first principles. We begin by deriving the vorticity transport equation directly from the Navier–Stokes equations, identifying the key physical mechanisms. We then prove Kelvin’s circulation theorem and Helmholtz’s vortex theorems — the cornerstones of inviscid vortex dynamics — and develop the Biot–Savart law for the velocity field induced by a given vorticity distribution. The chapter closes by extending the framework to rotating reference frames and introducing potential vorticity, the conserved quantity that anchors large-scale geophysical dynamics (taken up in detail in Chapter 8).

2.1 The Vorticity Equation

Derivation from the Navier–Stokes Equations

The vorticity equation is obtained by taking the curl of the incompressible Navier–Stokes equations. Starting from

\[ \frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla)\mathbf{u} = -\frac{1}{\rho}\nabla p + \nu\nabla^2\mathbf{u} + \mathbf{g}, \]

we use the vector identity \((\mathbf{u}\cdot\nabla)\mathbf{u} = \nabla(\tfrac{1}{2}|\mathbf{u}|^2) - \mathbf{u}\times\boldsymbol{\omega}\) to rewrite the advection term, then apply \(\nabla\times\) to the entire equation. Since \(\nabla\times(\nabla\phi) = \mathbf{0}\) for any scalar \(\phi\), the pressure gradient and the kinetic energy gradient both vanish under the curl. The gravitational term \(\mathbf{g} = -\nabla(gz)\) also vanishes. For a fluid of variable density, however, the pressure term contributes a baroclinic torque: the curl of \(-\nabla p / \rho\) does not vanish when \(\nabla \rho\) and \(\nabla p\) are not parallel. The result is the vorticity equation:

Each term on the right-hand side represents a distinct physical mechanism:

Vortex stretching and tilting — the term \((\boldsymbol{\omega}\cdot\nabla)\mathbf{u}\) is a \(3\times 3\) tensor acting on the vorticity vector. It has two geometric interpretations. When the component of \(\nabla\mathbf{u}\) along \(\boldsymbol{\omega}\) is extensional (the fluid element is being stretched along the vortex line), the vorticity magnitude increases by angular momentum conservation — this is vortex stretching. When the off-diagonal components are nonzero, the vortex line is tilted into a new orientation — this is vortex tilting. Vortex stretching is the primary mechanism by which turbulent flows generate intense small-scale vorticity from weaker large-scale vorticity, and it is the mechanism responsible for the tornado-like intensification of a vortex as it is drawn into a convergent region.

Baroclinic generation — the term \(\nabla\rho \times \nabla p / \rho^2\) (often called the baroclinic torque or solenoidal term) generates vorticity wherever surfaces of constant density are not parallel to surfaces of constant pressure. A classic example is the sea breeze: over a heated coastline, the density of air decreases faster over land than over sea, but the pressure gradient is nearly horizontal, so \(\nabla\rho \times \nabla p \neq \mathbf{0}\) and a coastal circulation is driven. For a barotropic fluid — one in which \(p = p(\rho)\), so that \(\nabla p\) is always parallel to \(\nabla\rho\) — the baroclinic term vanishes identically.

Viscous diffusion — the term \(\nu\nabla^2\boldsymbol{\omega}\) diffuses vorticity in exactly the same way that heat diffusion spreads a thermal anomaly. In the absence of the other terms, vorticity would simply spread out over a diffusion length \(\ell \sim \sqrt{\nu t}\). The dimensionless ratio of vortex stretching to viscous diffusion is the Reynolds number \(Re = UL/\nu\): at high \(Re\), diffusion is negligible and the inviscid dynamics dominate.

Two-Dimensional Flows

In two-dimensional flow in the \((x,y)\) plane, the velocity field is \(\mathbf{u} = (u,v,0)\) with no \(z\)-dependence, so the vorticity is purely in the \(z\)-direction: \(\boldsymbol{\omega} = \omega\hat{\mathbf{z}}\) with \(\omega = \partial v/\partial x - \partial u/\partial y\). The stretching term vanishes because \((\boldsymbol{\omega}\cdot\nabla)\mathbf{u} = \omega\, \partial\mathbf{u}/\partial z = \mathbf{0}\). The vorticity equation reduces to the 2D vorticity transport equation:

\[ \frac{D\omega}{Dt} = \nu\nabla^2\omega, \]

which is a scalar advection–diffusion equation. In the inviscid 2D limit, \(D\omega/Dt = 0\): vorticity is simply carried along by the flow without change. This is why 2D turbulence behaves so differently from 3D turbulence — the absence of vortex stretching prevents the forward cascade of enstrophy to small scales and instead permits an inverse energy cascade from small to large scales.

2.2 Kelvin’s Circulation Theorem

The Circulation Integral

The circulation \(\Gamma\) around a closed loop \(\mathcal{C}\) is defined as the line integral of velocity:

\[ \Gamma = \oint_{\mathcal{C}} \mathbf{u}\cdot d\boldsymbol{\ell}. \]

By Stokes’ theorem, \(\Gamma = \iint_{\mathcal{S}} \boldsymbol{\omega}\cdot d\mathbf{S}\), where \(\mathcal{S}\) is any surface bounded by \(\mathcal{C}\). Thus circulation measures the net vorticity flux threading through the loop — it is the “amount of spin” enclosed by the curve.

The key question is how \(\Gamma\) evolves when the loop \(\mathcal{C}(t)\) is a material loop — one whose points move with the fluid. The rate of change has two contributions: the change in the integrand \(\mathbf{u}\) as fluid accelerates, and the change in the path element \(d\boldsymbol{\ell}\) as the loop deforms.

Proof. The material derivative of the circulation is

\[ \frac{D\Gamma}{Dt} = \oint_{\mathcal{C}(t)} \frac{D\mathbf{u}}{Dt}\cdot d\boldsymbol{\ell} + \oint_{\mathcal{C}(t)} \mathbf{u}\cdot \frac{D(d\boldsymbol{\ell})}{Dt}. \]

For the second integral, the material derivative of the line element satisfies \(D(d\boldsymbol{\ell})/Dt = (d\boldsymbol{\ell}\cdot\nabla)\mathbf{u}\), so \(\mathbf{u}\cdot D(d\boldsymbol{\ell})/Dt = \mathbf{u}\cdot(d\boldsymbol{\ell}\cdot\nabla)\mathbf{u} = d(\tfrac{1}{2}|\mathbf{u}|^2)\), which integrates to zero around a closed loop. For the first integral, the Euler equation gives \(D\mathbf{u}/Dt = -\nabla p/\rho + \mathbf{g}\). The body force integrates to zero (conservative). The pressure term gives

\[ \oint_{\mathcal{C}} \frac{\nabla p}{\rho}\cdot d\boldsymbol{\ell} = \oint_{\mathcal{C}} \frac{dp}{\rho}. \]

For a barotropic fluid, \(p = p(\rho)\), so \(dp/\rho = dP\) is an exact differential of the function \(P(\rho) = \int dp/\rho(\cdot)\). Its integral around a closed loop vanishes. Hence \(D\Gamma/Dt = 0\). \(\square\)

Physical Implications

Kelvin’s theorem has profound consequences. If a flow is initially irrotational everywhere — \(\boldsymbol{\omega} = \mathbf{0}\) at \(t=0\) — then every material loop has zero circulation initially. By Kelvin’s theorem, every material loop maintains zero circulation for all time. Since this holds for every loop (and loops can be shrunk arbitrarily small), the vorticity remains zero everywhere for all time. An inviscid, barotropic flow that starts irrotational remains irrotational — the justification for potential flow theory (Chapter 3).

Conversely, if vorticity is created somewhere (by viscosity at a wall, or by the baroclinic term in a density-stratified flow), Kelvin’s theorem tells us that the vorticity cannot simply appear in the interior of an inviscid barotropic flow without violating the circulation constraint. Vorticity must enter through boundaries or through baroclinic generation.

2.3 Helmholtz’s Vortex Theorems

The vorticity field \(\boldsymbol{\omega}\) has zero divergence — \(\nabla\cdot\boldsymbol{\omega} = \nabla\cdot(\nabla\times\mathbf{u}) = 0\) always — and is therefore solenoidal. This constrains the geometry of vortex lines (curves tangent to \(\boldsymbol{\omega}\)) in the same way that \(\nabla\cdot\mathbf{B} = 0\) constrains magnetic field lines. Helmholtz’s vortex theorems, formulated in 1858, describe the kinematics and dynamics of vortex tubes in inviscid, barotropic flow.

This follows immediately from \(\nabla\cdot\boldsymbol{\omega} = 0\): by the divergence theorem applied to the volume between the two cross-sections, the net flux of \(\boldsymbol{\omega}\) through the closed surface is zero. Since the lateral surface of the tube contributes nothing (vortex lines are tangent to it), the two end fluxes must be equal and opposite. As a corollary, a vortex tube cannot end in the interior of a fluid — it must either close on itself (a vortex ring), end on a boundary, or extend to infinity.

This is simply Kelvin’s circulation theorem applied to the material loop that forms the boundary of a cross-section of the vortex tube. Combined with the first theorem, it says that the entire vortex tube has a fixed, time-independent circulation — vortex tubes are indestructible in ideal flow.

This is the statement that vorticity is frozen into the fluid in ideal flow. The mathematical proof shows that if \(\mathbf{x}(a,t)\) is the trajectory of a fluid particle initially at \(\mathbf{a}\), then the vorticity \(\boldsymbol{\omega}(\mathbf{x},t)\) satisfies the same equation as the material line element \(\partial\mathbf{x}/\partial\mathbf{a}\), so that vortex lines deform exactly as infinitesimal material line elements do.

Vortex Stretching as an Intensification Mechanism

Helmholtz’s theorems have a remarkable consequence for vortex dynamics. Because the strength \(\Gamma\) of a vortex tube is constant in time, while its cross-sectional area \(A\) can change (as the fluid stretches the tube), the average vorticity magnitude must change: \(\Gamma = \omega_{\text{avg}} \cdot A\), so if \(A\) decreases (the tube is stretched and thinned), then \(\omega_{\text{avg}}\) must increase proportionally. This is vortex stretching: a vortex that is stretched by the strain field of the surrounding flow intensifies. The phenomenon is visually familiar in the bathtub vortex, where water draining into a narrow outlet accelerates the angular velocity dramatically, and in tornadoes, where air converging toward the funnel stretches vertical vorticity into an intense concentrated core.

2.4 Vortex Dynamics: The Biot–Savart Law

Velocity Induced by a Vortex Distribution

The solenoidal condition \(\nabla\cdot\boldsymbol{\omega} = 0\) allows us to write \(\boldsymbol{\omega} = \nabla\times\mathbf{A}\) for a vector potential \(\mathbf{A}\). Combined with \(\boldsymbol{\omega} = \nabla\times\mathbf{u}\) and the Helmholtz decomposition, one can show that in an unbounded domain the velocity field induced by a given vorticity distribution satisfies the vector Poisson equation \(\nabla^2\mathbf{u} = -\nabla\times\boldsymbol{\omega}\). The solution is given by the Biot–Savart law, in exact analogy with magnetostatics:

The analogy with the Biot–Savart law of magnetostatics (where current \(\mathbf{J}\) plays the role of \(\boldsymbol{\omega}\)) is exact: vorticity acts as the “source” of the velocity field in the same way that electric current acts as the source of the magnetic field.

The Infinite Straight Vortex

For an infinite straight vortex filament of circulation \(\Gamma\) aligned with the \(z\)-axis, integrating the Biot–Savart law gives an azimuthal velocity

\[ u_\theta = \frac{\Gamma}{2\pi r}, \]

where \(r\) is the radial distance from the filament. This is the potential vortex of Chapter 3, recovered here from the fundamental integral law. The velocity field is irrotational everywhere except at \(r = 0\), where the vorticity is a delta function.

Self-Induced Motion of a Vortex Ring

A vortex ring of radius \(R\) and core radius \(a \ll R\) (and circulation \(\Gamma\)) translates along its symmetry axis at speed

\[ U_{\text{ring}} = \frac{\Gamma}{4\pi R}\left(\ln\frac{8R}{a} - \frac{1}{4}\right), \]

a result first obtained by Kelvin and derived by applying the Biot–Savart law to the toroidal geometry. The logarithm reflects the divergence of the self-induced velocity as the core size shrinks; real cores have finite size and internal structure that cut off this divergence. Smoke rings, which are vortex rings in air, translate at speeds of order \(\Gamma/(4\pi R)\) and expand slowly due to viscous diffusion.

Interaction of Two Parallel Vortices

Consider two parallel line vortices of circulations \(\Gamma_1\) and \(\Gamma_2\) separated by distance \(d\). Each vortex induces a velocity field at the location of the other, causing each to move.

Co-rotating vortices (\(\Gamma_1 = \Gamma_2 = \Gamma > 0\)): Each vortex moves in the velocity field of the other. The induced speed at distance \(d\) is \(U = \Gamma/(2\pi d)\), directed perpendicular to the line joining the two vortices. Both vortices therefore orbit their common centroid at angular velocity

\[ \Omega = \frac{\Gamma}{\pi d^2}, \]

tracing circular paths of radius \(d/2\). The combined system is a rotating doublet. In geophysical flows, merging of co-rotating vortices is a dominant mechanism for the inverse energy cascade in two-dimensional turbulence.

Counter-rotating vortices (\(\Gamma_1 = \Gamma > 0\), \(\Gamma_2 = -\Gamma\)): the two vortices translate together in the direction perpendicular to the line joining them at speed \(U = \Gamma/(2\pi d)\). This is the mechanism underlying the trailing vortex pair behind an aircraft wing — the two tip vortices of opposite sign propagate downward and away from the aircraft. A swimmer’s hands pushing backward form a counter-rotating pair that propels them forward.

2.5 The Vorticity Equation in Rotating Frames

The Coriolis Modification

In a frame rotating with angular velocity \(\boldsymbol{\Omega}\) (for Earth, \(|\boldsymbol{\Omega}| = 7.27\times 10^{-5}\) rad s\(^{-1}\)), the equations of motion gain Coriolis and centrifugal terms. The absolute vorticity — the vorticity measured in an inertial frame — is

\[ \boldsymbol{\omega}_{\text{abs}} = \boldsymbol{\omega} + 2\boldsymbol{\Omega}, \]

where \(\boldsymbol{\omega} = \nabla\times\mathbf{u}\) is the relative vorticity computed in the rotating frame and \(2\boldsymbol{\Omega}\) is the planetary vorticity. The factor of 2 arises from the relation \(\nabla\times(\boldsymbol{\Omega}\times\mathbf{x}) = 2\boldsymbol{\Omega}\).

The vorticity equation in a rotating frame, for an inviscid barotropic fluid, takes the form

\[ \frac{D\boldsymbol{\omega}_{\text{abs}}}{Dt} = (\boldsymbol{\omega}_{\text{abs}}\cdot\nabla)\mathbf{u}, \]

which is identical to the non-rotating form but with absolute vorticity replacing relative vorticity. This means that all of Kelvin’s and Helmholtz’s theorems carry over: the absolute circulation is conserved on material loops, and absolute vortex lines move with the fluid.

The Coriolis Parameter and \(f\)-Plane Approximation

On a sphere of radius \(a\) at latitude \(\phi\), the component of planetary vorticity in the local vertical direction is \(f = 2\Omega\sin\phi\), the Coriolis parameter. On the \(f\)-plane approximation (treating a small patch of Earth as flat with constant \(f\)), the relevant vorticity component is the vertical one:

\[ \zeta = \frac{\partial v}{\partial x} - \frac{\partial u}{\partial y}, \]

and the absolute vorticity is \(\zeta + f\).

Potential Vorticity Conservation

For shallow-water flow of depth \(H(\mathbf{x},t)\) in a rotating frame, combining conservation of absolute vorticity with conservation of mass leads to the most important invariant in geophysical fluid dynamics:

The physical interpretation is elegant: as a fluid column is stretched vertically (increasing \(H\)), it must spin up (increasing \(\zeta + f\)) to conserve potential vorticity — the fluid mechanical analogue of a figure skater pulling in their arms. Conversely, column squashing reduces vorticity. This mechanism is responsible for the generation of Rossby waves (§8.4), which propagate because a displaced fluid column creates a potential vorticity anomaly that drives the restoring motion.

Potential vorticity conservation is the organizing principle behind most of geophysical fluid dynamics. Its three-dimensional generalisation (Ertel’s potential vorticity theorem) plays the same role in the stratified atmosphere and ocean that it plays here in shallow water. Chapter 8 develops these ideas in detail; the present section is the conceptual foundation for that discussion.

Chapter 3: Potential Flow and Aerofoil Theory

3.1 Irrotational Flow and the Velocity Potential

Recall from AMATH 361: a flow is irrotational if \(\boldsymbol{\omega} = \nabla \times \mathbf{u} = \mathbf{0}\). By Kelvin’s circulation theorem, inviscid flows that start irrotational remain irrotational. For a simply connected domain, irrotationality implies the existence of a velocity potential \(\phi\) with \(\mathbf{u} = \nabla\phi\). Combined with incompressibility:

\[ \nabla^2\phi = 0 \qquad \text{(Laplace's equation)} \]

The theory of potential flow is therefore the theory of harmonic functions applied to fluid mechanics. Solutions are smooth, satisfy the superposition principle, and can be constructed by adding elementary solutions.

3.2 Elementary Potential Flows

A powerful feature of Laplace’s equation is linearity: superpositions of solutions are solutions. The elementary 2D complex potential solutions are:

Uniform stream: \(w(z) = Uz\), giving \(\phi = Ux\), \(\psi = Uy\), velocity \((U, 0)\).

Line source/sink of strength \(m\) at the origin: \(w(z) = (m/2\pi)\ln z\). Fluid radiates outward from (or converges toward) the origin like a point source in 2D.

Irrotational vortex of circulation \(\Gamma\): \(w(z) = -i(\Gamma/2\pi)\ln z\). Streamlines are circles; the speed \(|\mathbf{u}| = \Gamma/(2\pi r)\) is singular at the origin. This is the 2D analogue of a tornado vortex column.

Doublet (dipole) of strength \(\mu_d\): \(w(z) = \mu_d/(2\pi z)\). Obtained as the limit of a source-sink pair approaching each other with fixed product of strength and separation.

Flow past a circular cylinder: superposing a uniform stream of speed \(U\) and a doublet of strength \(\mu_d = Ua^2\) gives:

\[ w(z) = U\!\left(z + \frac{a^2}{z}\right) \]

The streamline \(\psi = 0\) passes through the circle \(|z| = a\), so this is the potential flow past a cylinder of radius \(a\). The velocity on the cylinder surface is:

\[ u_\theta = -2U\sin\theta \]

3.3 The Blasius Theorem and Lift

For a body in a potential flow, the net force per unit length can be computed from the complex potential without knowing the full flow field:

Applied to a cylinder with a superimposed irrotational vortex of circulation \(\Gamma\):

\[ X = 0, \qquad Y = \rho U \Gamma \]

This is the Kutta–Joukowski theorem: the lift per unit span is \(L' = \rho U \Gamma\). Zero drag (d’Alembert’s paradox — a consequence of neglecting viscosity and the resulting flow symmetry). The generation of lift requires circulation, and the Kutta condition (see below) determines the appropriate amount of circulation for a real aerofoil.

3.4 Conformal Mapping and Joukowski Aerofoils

A conformal map \(\zeta = f(z)\) is a complex analytic function that maps one flow domain to another while preserving angles. This is extraordinary: if we can solve Laplace’s equation in a simple geometry (a circle), a conformal map immediately gives us the solution in a complicated geometry (an aerofoil cross-section).

The Joukowski transformation:

\[ \zeta = z + \frac{c^2}{z} \]

maps a circle \(|z - z_0| = a\) (with \(z_0\) chosen appropriately) to a smooth curve in the \(\zeta\)-plane that resembles a wing cross-section. The trailing edge of the Joukowski aerofoil corresponds to the critical point of the map, where the derivative \(d\zeta/dz\) vanishes and the map is no longer conformal — the angle doubles, creating the characteristic cusped trailing edge.

The Kutta condition selects the physically correct solution from the family of potential flows past an aerofoil. Mathematically, it requires the velocity to remain finite at the trailing edge; this uniquely determines the circulation \(\Gamma\), and hence the lift:

\[ L' = \rho U \Gamma = \rho U \cdot 4\pi U a \sin(\alpha + \beta) \]

where \(\alpha\) is the angle of attack and \(\beta\) accounts for the aerofoil geometry. The physical mechanism: viscosity causes the rear stagnation point (which would otherwise sit on the upper surface) to move to the trailing edge, generating the correct circulation. Once the circulation is established, potential theory gives an excellent prediction of the lift even though viscosity has been neglected everywhere else.

---

Chapter 4: Surface Gravity Waves

4.1 The Linear Theory of Water Waves

Surface gravity waves are generated by any disturbance of a free water surface — a wind gust, a boat, an earthquake-driven seafloor motion. The restoring force is gravity, which acts to flatten disturbances of the free surface. We derive the linear (small-amplitude) wave theory for an inviscid, incompressible fluid with a free surface.

Consider a fluid of depth \(H\) (or infinite depth) occupying \(-H \leq z \leq 0\) in the rest state, with the free surface at \(z = \eta(x,t)\) in the disturbed state. The linearised free-surface boundary conditions (applied at \(z = 0\) rather than \(z = \eta\), since \(\eta\) is small) are:

• Kinematic: \(\partial\eta/\partial t = \partial\phi/\partial z\) — the surface moves with the normal velocity of the fluid beneath
• Dynamic: \(\partial\phi/\partial t + g\eta = 0\) — the pressure equals atmospheric pressure at the free surface (from Bernoulli’s equation)

Seeking wave solutions \(\phi \propto e^{ikx - i\omega t} Z(z)\), the Laplace equation \(\nabla^2\phi = 0\) gives \(Z'' - k^2 Z = 0\). The bottom boundary condition \(\partial\phi/\partial z = 0\) at \(z = -H\) gives \(Z(z) = \cosh k(z+H)\).

Combining the two free-surface conditions:

Two important limits follow immediately:

Deep water (\(kH \gg 1\), i.e., wavelength \(\ll\) depth): \(\tanh(kH) \to 1\), giving \(\omega^2 = gk\). Long waves travel faster than short waves; the group velocity \(c_g = \frac{1}{2}c_p\).

Shallow water (\(kH \ll 1\), i.e., wavelength \(\gg\) depth): \(\tanh(kH) \approx kH\), giving \(\omega^2 = gH k^2\), or \(\omega = \sqrt{gH}\,k\). All wavelengths travel at the same speed \(c = \sqrt{gH}\) — shallow-water waves are non-dispersive. Tsunamis are shallow-water waves (ocean depth \(\sim\) 4 km, wavelength \(\sim\) 200 km): they travel at \(\sqrt{gH} \approx 200\) m/s \(\approx\) 700 km/h.

Shallow (kH≪1)Deep (kH≫1)

kH=1

4.2 Group Velocity and Energy Transport

The concept of group velocity is one of the most profound in wave physics. For a dispersive medium with dispersion relation \(\omega = \omega(k)\), a localised wave packet — a superposition of plane waves with wavenumbers centred near \(k_0\) — propagates as a whole at the group velocity:

\[ c_g = \frac{d\omega}{dk}\bigg|_{k_0} \]

while individual crests move at the phase velocity \(c_p = \omega/k\). The group velocity is the velocity of energy transport. The split \(c_g = \frac{1}{2}c_p\) for deep-water gravity waves means that the energy packet travels at half the speed of the crests — a crest is born at the trailing edge of the packet, travels through to the leading edge, and disappears.

For shallow-water waves (\(c_g = c_p\)), there is no dispersive spreading: the wave packet propagates unchanged in shape. This is why a tsunami is devastating — it does not disperse during the thousands of kilometres of open-ocean travel.

4.3 Stokes Waves — A Glimpse of Nonlinearity

The linear theory of water waves assumes small amplitude \(a\) compared to wavelength \(\lambda\). The relevant parameter is the wave steepness \(\epsilon = ak\). For finite steepness, the free-surface boundary conditions must be applied at the actual free surface \(z = \eta\) rather than at \(z = 0\), introducing nonlinear corrections.

G.G. Stokes (1847) showed that the dispersion relation for finite-amplitude waves is modified:

\[ \omega^2 = gk\!\left(1 + \epsilon^2 + \cdots\right) \]

Higher-amplitude waves travel slightly faster. Moreover, the wave profile is no longer a pure sinusoid but has sharper crests and flatter troughs — the Stokes wave. As \(\epsilon \to \epsilon_{max} \approx 0.443\), the crest angle reaches 120° (the Stokes limiting wave) and the wave breaks. This nonlinear wave theory connects directly to the subject of AMATH 867.

---

Chapter 5: Laminar Flow — Exact Solutions

5.1 The Meaning of Exact Solutions

The Navier–Stokes equations are nonlinear PDEs and admit very few exact analytical solutions. When they do exist, the nonlinear advection term vanishes — usually because the geometry forces the flow to be unidirectional or because the streamlines are circles or spirals. These exact solutions are not mere mathematical exercises: they describe real flows that arise in engineering and nature, they provide benchmarks for numerical solvers, and they build physical intuition about how viscosity acts.

5.2 Plane Couette Flow

The simplest exact solution: fluid between two infinite parallel plates at \(y = 0\) and \(y = h\), with the lower plate at rest and the upper plate moving at speed \(U\). No imposed pressure gradient. The flow is \(\mathbf{u} = u(y)\hat{x}\), incompressibility is satisfied, and the Navier–Stokes equations reduce to:

\[ \nu\frac{d^2 u}{dy^2} = 0 \implies u(y) = U\frac{y}{h} \]

A perfect linear profile. The shear stress on each plate is \(\tau = \mu U/h\). Couette flow is the defining experiment for measuring viscosity: by measuring the torque on a rotating cylinder (the Couette viscometer), one infers \(\mu\) from the known geometry.

5.3 Hagen–Poiseuille Flow

Steady, pressure-driven flow in a straight circular pipe of radius \(R\) along the \(z\)-axis. Symmetry forces \(\mathbf{u} = u(r)\hat{z}\). In cylindrical coordinates, the Navier–Stokes equations reduce to:

\[ \frac{1}{r}\frac{d}{dr}\!\left(r\frac{du}{dr}\right) = \frac{1}{\mu}\frac{dp}{dz} \]

with no-slip \(u(R) = 0\) and regularity at the axis \(u'(0) = 0\). Integrating:

\[ u(r) = \frac{1}{4\mu}\!\left(-\frac{dp}{dz}\right)(R^2 - r^2) \]

The flow rate is:

\[ Q = \int_0^R u(r)\, 2\pi r\, dr = \frac{\pi R^4}{8\mu}\!\left(-\frac{dp}{dz}\right) \]

This \(R^4\) scaling — Hagen–Poiseuille’s law — has enormous practical consequences. A 10% reduction in pipe radius reduces flow rate by \((0.9)^4 \approx 66\%\) at the same pressure gradient. It governs blood flow in the circulatory system, water supply in municipal networks, and flow through porous media.

5.4 The Stokes Equations and Creeping Flow

At very low Reynolds number (\(Re \ll 1\)), the nonlinear advection term \((\mathbf{u}\cdot\nabla)\mathbf{u}\) is negligible compared to the viscous term. The Navier–Stokes equations linearise to:

\[ \nabla p = \mu\nabla^2\mathbf{u}, \qquad \nabla \cdot \mathbf{u} = 0 \]

These are the Stokes equations. They are linear, elliptic, and time-reversible (reversing the boundary velocities reverses the flow exactly). The famous consequence is that a swimming microorganism (bacterium, spermatozoon) cannot propel itself by a reciprocal stroke — it must break time-reversal symmetry, which is why microorganisms use rotating helical flagella rather than the back-and-forth paddling motion that works at high Reynolds number.

G.G. Stokes derived the drag on a sphere of radius \(a\) moving at speed \(U\) in Stokes flow:

\[ F_{drag} = 6\pi\mu a U \]

This formula underlies Millikan’s oil-drop experiment (1909), Stokes settling of particles in centrifugation, and the motion of aerosol particles in the atmosphere.

5.5 Oscillatory Stokes Layer

Consider the Stokes problem of a flat plate oscillating in its own plane at frequency \(\omega\): \(u(0,t) = U\cos\omega t\). The governing equation is:

\[ \frac{\partial u}{\partial t} = \nu\frac{\partial^2 u}{\partial y^2} \]

Seeking solutions of the form \(u = \text{Re}[\hat{u}(y) e^{-i\omega t}]\):

\[ u(y,t) = U e^{-\delta y}\cos(\omega t - \delta y), \qquad \delta = \sqrt{\frac{\omega}{2\nu}} \]

The velocity oscillation decays exponentially with a characteristic Stokes layer thickness \(\delta^{-1} = \sqrt{2\nu/\omega}\). For water at 1 Hz, \(\delta^{-1} \approx 0.6\) mm — the oscillatory boundary layer is very thin. This result is crucial for acoustic damping, oscillating bodies in water, and the dynamics of blood flow in the aorta.

---

Chapter 6: Boundary Layer Theory

6.1 The Boundary Layer Concept

One of the great discoveries of twentieth-century fluid mechanics is due to Ludwig Prandtl (1904): at high Reynolds number, viscosity is important only in a thin layer adjacent to solid boundaries. Away from this boundary layer, the flow is approximately inviscid and can be described by the Euler equations. Within the boundary layer, viscosity dominates and enforces the no-slip condition.

The thickness of the boundary layer \(\delta\) can be estimated by balancing the inertial and viscous terms in the Navier–Stokes equations:

\[ U\frac{\partial u}{\partial x} \sim U\frac{U}{L}, \qquad \nu\frac{\partial^2 u}{\partial y^2} \sim \nu\frac{U}{\delta^2} \]

Setting these equal: \(\delta \sim \sqrt{\nu L / U} = L/\sqrt{Re}\). The boundary layer becomes thinner as \(Re\) increases — which is why viscous effects appear to diminish at high Reynolds number even though \(\nu\) is fixed.

6.2 Prandtl’s Boundary Layer Equations

Prandtl derived simplified equations valid inside the boundary layer by scaling arguments. Let \(x\) be along the wall, \(y\) normal to it, \(U(x)\) the outer (inviscid) flow speed at the wall, \(u\) the streamwise velocity, and \(v\) the normal velocity. The full Navier–Stokes equations, after rescaling \(y \to y/\delta\) with \(\delta \sim Re^{-1/2}\), reduce in the limit \(Re \to \infty\) to:

The key simplification: the pressure gradient across the boundary layer is negligible (\(\partial p / \partial y \approx 0\)), so the pressure at the wall equals the inviscid outer pressure, given by Bernoulli’s equation along the outer streamline: \(dP/dx = -\rho U\, dU/dx\).

6.3 The Blasius Flat-Plate Solution

For a flat plate (\(dU/dx = 0\), so the outer flow speed \(U\) is constant), the boundary layer equations admit a similarity solution. The key observation is that the equations have no built-in length scale when \(U = \text{const}\). This suggests that the velocity profile at different \(x\)-locations is self-similar: rescaling \(y\) by \(\delta(x) \sim \sqrt{\nu x/U}\) collapses all profiles onto a single curve.

Introducing the similarity variable \(\eta = y\sqrt{U/(\nu x)}\) and stream function \(\psi = \sqrt{\nu U x}\, f(\eta)\), the boundary layer equations reduce to the Blasius ordinary differential equation:

\[ f''' + \frac{1}{2}ff'' = 0, \qquad f(0) = f'(0) = 0, \quad f'(\infty) = 1 \]

This third-order nonlinear ODE must be solved numerically (a shooting method works well). The boundary layer thickness (defined as the 99% thickness where \(u = 0.99U\)) is:

\[ \delta_{99}(x) \approx \frac{5.0\, x}{\sqrt{Re_x}}, \qquad Re_x = \frac{Ux}{\nu} \]

The skin friction coefficient (drag per unit area normalised by dynamic pressure):

\[ c_f = \frac{\tau_w}{\frac{1}{2}\rho U^2} = \frac{0.664}{\sqrt{Re_x}} \]

These are exact asymptotic results, valid for large \(Re_x\). They agree very well with experiments up to the transition to turbulent boundary layers (typically around \(Re_x \approx 5\times 10^5\)).

6.4 Boundary Layer Separation

When the outer flow decelerates (\(dU/dx < 0\)), the adverse pressure gradient (\(dP/dx > 0\)) decelerates the boundary layer fluid. Near the wall, where the fluid has been slowed by friction, the momentum may be insufficient to overcome the adverse pressure gradient, and the flow reverses direction. At the point where the wall shear stress vanishes, the boundary layer is said to separate.

After separation, the thin boundary layer approximation breaks down, the flow leaves the surface, and a large recirculating wake forms. This catastrophic phenomenon — boundary layer separation — is responsible for:

• The dramatic drag increase on a bluff body (a sphere, a car at high speed)
• The stall of an aircraft wing at high angle of attack
• Flow separation at pipe bends and over hills

The position of separation depends delicately on the surface geometry and the upstream pressure distribution. For a circular cylinder in potential flow, the outer flow decelerates after the midpoint, and separation occurs around 80° from the leading stagnation point. The dramatic drag reduction of a dimpled golf ball (compared to a smooth one) works by triggering an early turbulent transition in the boundary layer, which can sustain higher adverse pressure gradients before separating — the separation point moves backward and the wake narrows.

---

Chapter 7: Turbulence

7.1 The Nature of Turbulent Flow

At sufficiently high Reynolds number, laminar flows become unstable and transition to turbulence: an apparently chaotic, three-dimensional, time-dependent motion characterised by a broad range of active length and time scales. Turbulence is not random — it is deterministic at every instant (the Navier–Stokes equations have a unique solution for smooth initial data) — but it is extremely sensitive to initial conditions, making long-time prediction effectively impossible.

The Reynolds number at which transition occurs depends on geometry and the level of background disturbances:

• Pipe flow: \(Re \approx 2300\) (but can be delayed to \(10^5\) in very quiet conditions)
• Flat-plate boundary layer: \(Re_x \approx 5 \times 10^5\)
• Channel flow: \(Re \approx 2000\)

Turbulence is characterised by vortex stretching (absent in 2D) transferring energy to progressively smaller scales until viscous dissipation converts kinetic energy to heat. This energy cascade (Kolmogorov 1941) determines the statistics of turbulent flows.

7.2 The Reynolds Decomposition

Osborne Reynolds (1895) introduced the decomposition that still bears his name: split each quantity into a mean (time-averaged) and a fluctuation:

\[ u_i = \overline{u}_i + u'_i, \qquad p = \bar{p} + p' \]

where \(\overline{u'_i} = 0\) by definition of the mean. Substituting into the Navier–Stokes equations and averaging:

The term \(-\rho\overline{u'_i u'_j}\) is the Reynolds stress tensor: the mean momentum flux due to turbulent fluctuations. It plays the same formal role as the viscous stress, but it is orders of magnitude larger in fully turbulent flows. The fundamental difficulty of turbulence theory — the closure problem — is that the RANS equations for the means involve unknowns (the Reynolds stresses), and equations for the Reynolds stresses involve higher-order correlations, and so on without end.

7.3 Turbulent Kinetic Energy and the Closure Problem

The turbulent kinetic energy (TKE) is:

\[ k = \frac{1}{2}\overline{u'_i u'_i} \]

Its transport equation (derived by multiplying the fluctuation equation by \(u'_i\) and averaging) takes the form:

\[ \frac{Dk}{Dt} = \mathcal{P} - \varepsilon + \nabla \cdot (\text{diffusive flux}) \]

where \(\mathcal{P} = -\overline{u'_i u'_j}\,\partial\overline{u}_i/\partial x_j\) is the production of TKE by interaction of Reynolds stresses with the mean flow shear, and \(\varepsilon = \nu\overline{(\partial u'_i/\partial x_j)^2}\) is the viscous dissipation rate — the rate at which turbulent kinetic energy is converted to heat.

The closure problem is acute for the dissipation rate. Kolmogorov’s 1941 theory provides a statistical framework: in the inertial subrange (scales between the energy-containing large eddies and the dissipation-scale Kolmogorov eddies), the energy spectrum follows:

\[ E(k) = C\varepsilon^{2/3} k^{-5/3} \]

where \(C \approx 1.5\) is the Kolmogorov constant. This \(k^{-5/3}\) power law is one of the most universally observed results in fluid turbulence, confirmed over many decades of wavenumber in atmospheric, oceanic, and laboratory experiments.

EnergyinjectionInertialsubrangeDissipation

k_Lk_η

E(k) = Cε²/³k⁻⁵/³

7.4 Wall Turbulence and the Logarithmic Layer

Near a solid wall, turbulence has a characteristic multi-layer structure. Define the friction velocity \(u_* = \sqrt{\tau_w/\rho}\) (where \(\tau_w\) is the wall shear stress) and the viscous length \(\ell_\nu = \nu/u_*\). The dimensionless wall distance is \(y^+ = y/\ell_\nu = u_* y/\nu\).

Three distinct regions exist:

Viscous sublayer (\(y^+ \lesssim 5\)): turbulent fluctuations are suppressed by the wall, viscous stress dominates, and the velocity profile is linear: \(u^+ = y^+\) where \(u^+ = \overline{u}/u_*\).

Buffer layer (\(5 \lesssim y^+ \lesssim 30\)): transition region where both viscous and turbulent stresses are comparable.

Logarithmic (inertial) layer (\(y^+ \gtrsim 30\), \(y \ll \delta\)): both viscosity and outer boundary conditions are irrelevant. By dimensional analysis alone, the only possible form of the velocity gradient is:

\[ \frac{d\overline{u}}{dy} = \frac{u_*}{\kappa y} \]

Integrating:

\[ u^+ = \frac{1}{\kappa}\ln y^+ + B \]

with von Kármán constant \(\kappa \approx 0.41\) and intercept \(B \approx 5.2\). This logarithmic law of the wall is remarkably universal: it holds for smooth and rough walls, pipes and boundary layers, gases and liquids, from laboratory experiments to atmospheric surface layers.

---

Chapter 8: Geophysical Fluid Dynamics

8.1 The Rotating Reference Frame

The dynamics of Earth’s atmosphere and oceans cannot be understood without accounting for the rotation of the Earth. In a reference frame rotating at angular velocity \(\boldsymbol{\Omega}\), additional apparent forces arise:

The Coriolis force acts perpendicular to the velocity: in the Northern Hemisphere it deflects moving air to the right; in the Southern Hemisphere, to the left. This deflection is responsible for the counter-clockwise rotation of Northern Hemisphere cyclones (low-pressure systems) and the clockwise rotation of anticyclones. It is also responsible for trade winds, ocean gyres, and the jet streams.

On a local \(f\)-plane approximation at latitude \(\phi\), the local vertical component of the Coriolis parameter is:

\[ f = 2\Omega\sin\phi \]

At the North Pole \(f = 2\Omega \approx 1.45 \times 10^{-4}\) s\(^{-1}\); at the equator \(f = 0\).

8.2 Geostrophic Balance

At large scales in the atmosphere and ocean (\(L \sim 10^3\) km, \(U \sim 10\) m/s), the Rossby number:

\[ Ro = \frac{U}{fL} \sim \frac{10}{10^{-4} \times 10^6} = 0.1 \]

is small. To leading order in \(Ro\), the Coriolis force balances the pressure gradient: this is geostrophic balance.

In the atmosphere, low-pressure systems (cyclones) have counter-clockwise geostrophic circulation in the Northern Hemisphere: wind flows around the low, not into it. This is completely opposite to what you might naively expect (wind blowing from high to low pressure) — and is entirely due to the Coriolis force.

8.3 The Taylor–Proudman Theorem

A remarkable consequence of rotation concerns three-dimensional flow at small Rossby number.

The experimental demonstration by G.I. Taylor (1921) is beautiful: obstacles placed in a rotating tank of water generate columns of fluid that extend vertically far above the obstacle, deflecting the flow horizontally as if the column were a solid cylinder. The fluid “sees” the obstacle even where it cannot see it directly. This column rigidity has profound consequences for the dynamics of the Earth’s liquid outer core and for the structure of large-scale atmospheric and oceanic eddies.

8.4 Rossby Waves

The most important wave in large-scale geophysical fluid dynamics is the Rossby wave (or planetary wave), driven by the variation of the Coriolis parameter with latitude. On a \(\beta\)-plane approximation, \(f \approx f_0 + \beta y\) where \(\beta = df/dy = 2\Omega\cos\phi/R_E\) (with \(R_E\) the Earth’s radius).

Considering linearised perturbations about a state of rest with barotropic (depth-independent) flow, the vorticity equation gives the Rossby wave dispersion relation:

Two features are immediately striking. First, \(\omega < 0\) always (since \(\beta > 0\) in both hemispheres and the phase speed in the \(k\)-direction is \(c_p = \omega/k < 0\)): Rossby waves always propagate westward relative to the mean flow. Second, the group velocity in the \(x\)-direction is:

\[ c_{gx} = \frac{\partial\omega}{\partial k} = \frac{\beta(k^2 - l^2)}{(k^2 + l^2 + f_0^2/c^2)^2} \]

which can be positive or negative depending on the wavenumber. Long Rossby waves carry energy westward; short Rossby waves carry energy eastward.

8.5 Stratification and Internal Waves

In addition to rotation, geophysical fluids are typically stratified: their density varies with depth due to temperature and salinity gradients. The buoyancy force on a vertically displaced fluid parcel provides a restoring mechanism that supports internal gravity waves.

Define the buoyancy frequency (Brunt–Väisälä frequency):

\[ N^2 = -\frac{g}{\rho_0}\frac{d\bar{\rho}}{dz} \]

where \(\bar{\rho}(z)\) is the background density profile. When \(N^2 > 0\) (density decreasing upward — stably stratified), a vertically displaced parcel oscillates at frequency \(N\); when \(N^2 < 0\) (density increasing upward), the stratification is unstable and convection occurs.

The dispersion relation for internal gravity waves in a uniformly stratified fluid is:

\[ \omega^2 = N^2 \frac{k^2 + l^2}{k^2 + l^2 + m^2} = N^2\sin^2\theta \]

where \(m\) is the vertical wavenumber and \(\theta\) is the angle of the wave vector from the vertical. This has a remarkable geometric interpretation: the frequency depends only on the angle of propagation, not the magnitude of the wavenumber. Equally remarkable: the group velocity is perpendicular to the phase velocity — energy travels horizontally while wave crests tilt. Internal waves are responsible for the stirring of the deep ocean interior, the generation of microstructure at density interfaces, and the mixing that drives the global thermohaline circulation.

---

---

Chapter 9: Compressible Flow and Sound Waves

All the flow physics developed in Chapters 1–8 assumed that the fluid density is either constant (incompressible liquids) or at most weakly varying (the Boussinesq approximation for stratified flows). This assumption breaks down whenever the fluid velocity becomes comparable to the speed of sound, or whenever thermal effects couple significantly to the dynamics. The resulting field — compressible flow — is the domain of aeroacoustics, high-speed aerodynamics, shock physics, and atmospheric acoustics.

This chapter follows Kundu Chapter 16. We begin with the thermodynamic framework needed to close the equations of motion for a compressible gas, then derive the compressible Euler equations and linearise them to recover the acoustic wave equation. We examine the energy and intensity of sound waves, discuss the Mach number as the key parameter governing compressibility effects, and close with the physics of shock waves — the most dramatic consequence of nonlinear compressible dynamics.

9.1 Thermodynamics of Gases

The Ideal Gas Law

A perfect (ideal) gas satisfies the equation of state

\[ p = \rho R T, \]

where \(p\) is the absolute pressure, \(\rho\) the density, \(T\) the absolute temperature (Kelvin), and \(R\) the specific gas constant (\(R = R_{\text{univ}}/M\), with \(R_{\text{univ}} = 8.314\) J mol\(^{-1}\) K\(^{-1}\) and \(M\) the molar mass). For dry air, \(R \approx 287\) J kg\(^{-1}\) K\(^{-1}\). The ideal gas law is an excellent approximation for gases at moderate pressures and temperatures — conditions met throughout atmospheric dynamics and most engineering flows.

First and Second Laws of Thermodynamics

The first law of thermodynamics for a fluid parcel states that the change in specific internal energy \(e\) equals the heat added minus the work done by the parcel:

\[ De = \delta q - p\, D\!\left(\frac{1}{\rho}\right) = \delta q + \frac{p}{\rho^2}D\rho, \]

where \(\delta q\) is the specific heat added and \(1/\rho\) is the specific volume. For an ideal gas, \(e = c_v T\) where \(c_v\) is the specific heat at constant volume. The specific enthalpy is \(h = e + p/\rho = c_p T\), where \(c_p = c_v + R\) is the specific heat at constant pressure.

The second law introduces the specific entropy \(s\): for a reversible process, \(Ds = \delta q / T \geq 0\) (with equality only for reversible processes). For an ideal gas:

\[ Ds = c_v \frac{DT}{T} - R\frac{D\rho}{\rho} = c_p\frac{DT}{T} - \frac{R}{p}Dp. \]

Isentropic Relations

An isentropic (adiabatic and reversible) process has \(Ds = 0\). Setting the entropy expression to zero and integrating:

\[ \frac{p}{\rho^\gamma} = \text{const}, \quad \text{or equivalently} \quad p \propto \rho^\gamma, \]

where \(\gamma = c_p/c_v\) is the ratio of specific heats (\(\gamma = 7/5 = 1.4\) for diatomic gases such as air at moderate temperatures). These isentropic (Poisson) relations connect the thermodynamic variables along a fluid parcel trajectory in the absence of heat transfer.

The Speed of Sound

Consider a small pressure perturbation propagating through a gas at rest. The propagation speed is determined by how the gas responds to compression. For an isentropic process (sound waves are rapid enough that heat conduction is negligible):

\[ c^2 = \left.\frac{\partial p}{\partial \rho}\right|_s = \frac{\gamma p}{\rho} = \gamma R T. \]

The square-root dependence on temperature means that the speed of sound increases with temperature: sound travels faster in warm air near the ground than in the cold upper atmosphere, leading to refraction of sound toward the ground in temperature inversions. The isentropic assumption is critical: the isothermal speed \(\sqrt{RT}\) (Newton’s original estimate) is lower by a factor of \(\sqrt{\gamma}\), a discrepancy that puzzled scientists until Laplace correctly identified the adiabatic character of sound propagation in 1816.

9.2 The Compressible Euler Equations

The Full System

For an inviscid compressible gas, conservation of mass, momentum, and energy gives the compressible Euler equations:

The first equation is the mass continuity equation — now in its fully compressible form, with no simplification to \(\nabla\cdot\mathbf{u} = 0\). The second is the momentum equation in conservative form; expanding the divergence and using mass conservation recovers \(\rho D\mathbf{u}/Dt = -\nabla p\). The third is the energy equation: the quantity in parentheses on the left is the total energy density (internal plus kinetic), and the flux on the right is the total energy flux plus work done by pressure forces.

Isentropic Flow

For smooth (shock-free) adiabatic flows with no heat sources, the entropy of each fluid parcel is conserved: \(Ds/Dt = 0\). The energy equation is then equivalent to the isentropic relation \(p\propto\rho^\gamma\), and the system reduces to four equations (mass, three momentum) for the four unknowns \((\rho, \mathbf{u})\) with \(p = p(\rho) = p_0(\rho/\rho_0)^\gamma\). This simplification is valid for most acoustic phenomena and for smooth supersonic flows away from shocks.

9.3 Small-Amplitude Sound Waves

Linearisation about a Uniform State

We linearise the compressible Euler equations about a uniform rest state \((\rho_0, p_0, \mathbf{u}_0 = \mathbf{0})\) by writing \(\rho = \rho_0 + \rho'\), \(p = p_0 + p'\), \(\mathbf{u} = \mathbf{u}'\), where primed quantities are small. Substituting and retaining only first-order terms:

\[ \frac{\partial \rho'}{\partial t} + \rho_0\nabla\cdot\mathbf{u}' = 0, \]\[ \rho_0\frac{\partial \mathbf{u}'}{\partial t} = -\nabla p'. \]

For isentropic perturbations, \(p' = c_0^2\rho'\) where \(c_0^2 = \gamma p_0/\rho_0\). Substituting into the continuity equation and differentiating with respect to time:

\[ \frac{\partial^2 \rho'}{\partial t^2} = -\rho_0\frac{\partial(\nabla\cdot\mathbf{u}')}{\partial t} = \nabla\cdot(\nabla p') = c_0^2\nabla^2\rho'. \]

Plane Wave Solutions and Dispersion Relation

A plane wave solution \(p' = P\exp(i\mathbf{k}\cdot\mathbf{x} - i\omega t)\) satisfies the wave equation if and only if

\[ \omega^2 = c_0^2|\mathbf{k}|^2, \quad \text{i.e.,} \quad \omega = c_0 k, \]

where \(k = |\mathbf{k}|\). This is a non-dispersive dispersion relation: the phase velocity \(\omega/k = c_0\) and the group velocity \(\partial\omega/\partial k = c_0\) are both equal to the sound speed and independent of wavenumber. All frequencies travel at the same speed — sound pulses do not spread out as they propagate (unlike water waves). This is why the pitch of a distant instrument is unaffected by distance, and why sonar can achieve sharp time resolution.

The velocity perturbation associated with the plane wave is

\[ \mathbf{u}' = \frac{p'}{\rho_0 c_0}\hat{\mathbf{k}}, \]

where \(\hat{\mathbf{k}}\) is the unit vector in the propagation direction. Sound waves are longitudinal — the fluid oscillates parallel to the direction of wave propagation, in contrast to transverse waves (e.g., electromagnetic waves or shear waves in solids).

9.4 Energy and Intensity of Acoustic Waves

Acoustic Energy Density

For a plane sound wave, the time-averaged acoustic energy density (energy per unit volume) is

\[ \langle E \rangle = \frac{\langle p'^2\rangle}{\rho_0 c_0^2} = \rho_0\langle|\mathbf{u}'|^2\rangle, \]

consisting of equal contributions from potential energy (compression) and kinetic energy (particle motion), in analogy with the equipartition of energy in a simple harmonic oscillator.

Acoustic Intensity and Energy Flux

The acoustic intensity \(\mathbf{I}\) is the energy flux — the rate at which energy is transported per unit area:

\[ \mathbf{I} = p'\mathbf{u}' = \frac{\langle p'^2\rangle}{\rho_0 c_0}\hat{\mathbf{k}}. \]

This is the acoustic analogue of the Poynting vector in electromagnetism. The time-averaged intensity for a sinusoidal wave of amplitude \(P\) is

\[ \langle I\rangle = \frac{P^2}{2\rho_0 c_0}. \]

The quantity \(Z = \rho_0 c_0\) is the acoustic impedance of the medium; it controls the reflection and transmission of sound at interfaces between different materials (the origin of echo, sonar returns, and the muffled sound when one’s ears are submerged).

The Inverse-Square Law

For a point source of acoustic power \(\mathcal{P}\) in a free field, the wave energy spreads uniformly over spherical surfaces. At radius \(r\) from the source, the surface area is \(4\pi r^2\), so by energy conservation:

\[ \langle I \rangle = \frac{\mathcal{P}}{4\pi r^2}. \]

This inverse-square law is the basis of the decibel scale in acoustics: each doubling of distance reduces intensity by 6 dB. The inverse-square law holds for free-field propagation; real environments involve reflections, absorption, and diffraction, but the \(r^{-2}\) decay remains the dominant trend.

9.5 The Mach Number and Compressibility Effects

The Mach Number

The Mach number is the ratio of a characteristic flow speed \(U\) to the local speed of sound:

\[ Ma = \frac{U}{c}. \]

It is the fundamental parameter governing compressibility effects.

The Breakdown of the Incompressible Approximation

For flow speeds much smaller than the sound speed, compressibility effects are negligible and the incompressible approximation is excellent. The fractional density change in a flow of speed \(U\) over an obstacle is of order \(Ma^2\):

\[ \frac{\Delta\rho}{\rho_0} \sim Ma^2. \]

For \(Ma < 0.3\), this is less than 9% — negligible for engineering purposes. Above \(Ma \approx 0.3\), compressibility corrections become measurable, and above \(Ma \approx 0.8\), they become dominant.

The Prandtl–Glauert Compressibility Correction

For subsonic flow past an aerofoil at Mach number \(Ma\), the lift coefficient is enhanced relative to its incompressible value by the Prandtl–Glauert rule:

\[ C_L = \frac{C_{L,0}}{\sqrt{1 - Ma^2}}, \]

where \(C_{L,0}\) is the incompressible lift coefficient. This correction follows from a linearised potential flow analysis: compressibility effectively increases the camber of the aerofoil as seen by the flow. The rule diverges as \(Ma \to 1\), signalling the breakdown of linear theory and the onset of transonic wave drag — the “sound barrier” encountered by early jet aircraft.

9.6 Shock Waves in Compressible Flow

Steepening of Finite-Amplitude Waves

In linear acoustics, all Fourier components travel at the same speed \(c_0\), and a waveform propagates without distortion. For finite-amplitude waves, however, the local speed of sound depends on the local density: \(c = c_0 + (\gamma+1)u'/2 + O(u'^2)\). Compressed regions (higher \(c\)) overtake rarefied regions (lower \(c\)), causing the waveform to progressively steepen. Eventually, the gradient \(\partial u'/\partial x\) becomes infinite in finite time — the nonlinear wave breaks, forming a shock wave. This is exactly the wave-breaking mechanism analysed in §10.3 for the inviscid Burgers equation, here in a compressible gas context.

A shock is an extremely thin (of order mean free path) layer across which the flow variables change discontinuously on macroscopic scales. The internal structure is controlled by viscosity and heat conduction, but the jump conditions are determined purely by conservation laws.

The Rankine–Hugoniot Relations

Consider a normal shock — a shock perpendicular to the incoming flow — with the shock stationary and the upstream state \((\rho_1, u_1, p_1, T_1)\) entering from the left and the downstream state \((\rho_2, u_2, p_2, T_2)\) exiting to the right. Conservation of mass, momentum, and energy across the shock (with no work done or heat added in the shock itself) gives:

These three equations plus the equation of state determine the four downstream variables given the upstream state and the shock speed. The key result is the normal shock relations expressed in terms of the upstream Mach number \(Ma_1 > 1\):

\[ \frac{\rho_2}{\rho_1} = \frac{(\gamma+1)Ma_1^2}{(\gamma-1)Ma_1^2 + 2}, \qquad \frac{p_2}{p_1} = \frac{2\gamma Ma_1^2 - (\gamma-1)}{\gamma+1}, \qquad \frac{T_2}{T_1} = \frac{p_2}{p_1}\cdot\frac{\rho_1}{\rho_2}. \]

The downstream Mach number satisfies \(Ma_2 < 1\): a normal shock always converts supersonic flow to subsonic flow. For strong shocks (\(Ma_1 \to \infty\)):

\[ \frac{\rho_2}{\rho_1} \to \frac{\gamma+1}{\gamma-1} = 6 \quad (\text{for air}), \qquad \frac{p_2}{p_1} \to \frac{2\gamma}{\gamma+1}Ma_1^2. \]

The density ratio is bounded (by at most a factor of 6 for a diatomic gas), while the pressure ratio grows without limit.

Entropy Production at a Shock

The Rankine–Hugoniot relations are consistent with the second law of thermodynamics only for shocks moving into supersonic flow (\(Ma_1 > 1\)). The entropy jump across the shock is

\[ s_2 - s_1 = c_v\ln\frac{p_2/p_1}{(p_2/p_1)^{\gamma}} = c_v\ln\frac{T_2^\gamma}{p_2^{\gamma-1}} - c_v\ln\frac{T_1^\gamma}{p_1^{\gamma-1}} > 0, \]

which is positive for \(Ma_1 > 1\) and negative (physically forbidden) for \(Ma_1 < 1\). Thus the second law selects only compressive shocks, ruling out expansion shocks. The entropy production represents irreversible dissipation within the shock layer — even though the Rankine–Hugoniot conditions do not explicitly involve viscosity, the shock owes its existence to molecular dissipation on the sub-mean-free-path scale.

Oblique Shocks

When a supersonic flow encounters a wedge-shaped surface at an angle, the resulting shock is oblique — inclined at an angle to the incoming flow direction. The normal Rankine–Hugoniot relations still apply to the velocity component normal to the shock; the tangential component is unchanged. Oblique shocks produce weaker jumps than normal shocks at the same upstream Mach number, and they are the predominant shock type in practical supersonic flows (over wing leading edges, in engine intakes, and in supersonic jets). The shock-wave/boundary-layer interaction at oblique shocks is an important and still-active research problem in high-speed aerodynamics.

Chapter 10: Nonlinear Waves — AMATH 867

This chapter extends the course into graduate territory, covering the topics of AMATH 867 (Nonlinear Waves). The central question: what happens when the small-amplitude assumption of linear wave theory breaks down?

10.1 Beyond Linearity — The Need for Nonlinear Theory

Linear wave theory, which dominated Chapters 3 and 7, rests on the assumption that wave amplitudes are small enough that nonlinear terms in the governing equations can be neglected. This approximation captures the dispersion relation, the group velocity, and the spatial structure of waves — but it fundamentally misses phenomena that arise from wave-wave interactions:

• A wave steepens into a shock (as in a sonic boom or a breaking ocean wave)
• Two pulses pass through each other and emerge unchanged (solitons)
• A broad spectrum of waves redistributes energy across wavenumbers (wave turbulence)
• Modulational instability breaks a uniform wavetrain into localised packets (the Benjamin–Feir instability of ocean waves)

The mathematics of nonlinear waves draws on the method of characteristics, weakly nonlinear asymptotics, and the theory of integrable systems. The goal is not merely to find solutions but to understand the qualitative character of wave dynamics when amplitude matters.

10.2 Method of Characteristics and First-Order Conservation Laws

The simplest nonlinear wave equation is the inviscid Burgers equation (or the scalar conservation law):

\[ \frac{\partial u}{\partial t} + u\frac{\partial u}{\partial x} = 0 \]

This is the nonlinear advection equation: each value of \(u\) advects at its own speed \(u\). Unlike the linear advection equation \(\partial_t u + c\, \partial_x u = 0\) (where all parts of the signal move at the same speed \(c\)), faster-moving parts of the wave profile overtake slower-moving parts, causing the wave to steepen.

The method of characteristics converts this PDE into a family of ODEs. Define characteristic curves \(x = x(t)\) along which:

\[ \frac{dx}{dt} = u, \qquad \frac{du}{dt} = 0 \]

Along each characteristic, \(u\) is constant. So characteristic \(C\) starting at \(x_0\) carries the initial value \(u_0 = u(x_0, 0)\) and has equation \(x = x_0 + u_0 t\). The solution is implicitly:

\[ u(x,t) = u_0(x - ut) \]

Characteristics are straight lines in the \((x,t)\) plane. When characteristics from different initial points converge — when \(\partial u_0/\partial x_0 < 0\) somewhere — they eventually cross at a finite time:

\[ t_{break} = \frac{-1}{\min_x\bigl(\partial u_0/\partial x\bigr)} \]

At \(t = t_{break}\), the solution becomes multi-valued: the wave has broken. Physically, the steep wave front has steepened to a vertical slope and then to an overhang — which is unphysical. The solution beyond breaking must be interpreted as a shock: a discontinuity satisfying a jump condition (the Rankine–Hugoniot relation) derived from the integral conservation law.

10.3 Shock Formation and the Rankine–Hugoniot Condition

A shock at position \(x = s(t)\) satisfies the integral conservation law for the original PDE. For the scalar conservation law \(\partial_t u + \partial_x F(u) = 0\):

\[ \frac{ds}{dt} = \frac{F(u_R) - F(u_L)}{u_R - u_L} \]

where \(u_L\) and \(u_R\) are the values of \(u\) immediately to the left and right of the shock. For Burgers’ equation \(F(u) = u^2/2\), this gives shock speed \(s'(t) = (u_L + u_R)/2\) — the average of the speeds on either side.

An additional condition — the entropy condition — selects the physically correct (stable) shock from among all weak solutions. For Burgers’ equation: \(u_L > u_R\) (faster fluid behind, slower ahead — this is the steepening condition). Shocks satisfying the entropy condition are stable and those violating it are unstable (they split into rarefactions).

10.4 Weakly Nonlinear Waves — The Stokes Expansion

For wave problems, the characteristic approach gives shocks, but many physical waves do not shock: they are dispersive, meaning different wavenumbers travel at different speeds, and dispersion can balance nonlinear steepening. The regime of weakly nonlinear waves is where both effects are comparable.

The strategy is to expand the wave amplitude in a small parameter \(\epsilon\):

\[ u = \epsilon u_1 + \epsilon^2 u_2 + \epsilon^3 u_3 + \cdots \]

At order \(\epsilon\), we recover the linear solution. At order \(\epsilon^2\), the nonlinear terms in the original equation drive a correction. Secular terms (terms that grow without bound in time) arise at successive orders unless the linear solution’s amplitude and phase are allowed to evolve slowly.

This is the method of multiple scales or modulation theory. We introduce a slow time \(T = \epsilon t\) and slow space \(X = \epsilon x\) alongside the fast oscillatory variables. The linear solution has the form:

\[ u_1 = A(X,T)\, e^{i(kx - \omega t)} + \text{c.c.} \]

where the complex amplitude \(A(X,T)\) varies slowly. Requiring the \(\epsilon^2\) and \(\epsilon^3\) equations to be non-secular imposes equations on \(A\). At leading order: \(A\) propagates at the group velocity, \(\partial A/\partial T + c_g\, \partial A/\partial X = 0\). At next order, the interplay of dispersion and nonlinearity yields the nonlinear Schrödinger equation (NLS):

\[ i\frac{\partial A}{\partial T} + \frac{1}{2}\omega''\frac{\partial^2 A}{\partial X^2} + \gamma|A|^2 A = 0 \]

where \(\omega'' = d^2\omega/dk^2\) is the dispersion coefficient and \(\gamma\) is a nonlinear coefficient depending on the specific physical problem. The NLS governs the slow evolution of the wave envelope for a broad class of dispersive nonlinear wave problems, from water waves to optical pulses in fibres.

10.5 Modulational Instability (Benjamin–Feir)

The NLS has a remarkable consequence: a uniform wavetrain \(A = A_0 e^{i\gamma A_0^2 T}\) (a wave of constant amplitude) is unstable to long-wavelength perturbations when \(\omega'' \gamma > 0\). This is the Benjamin–Feir instability (1967), discovered simultaneously by Lighthill.

The instability works as follows: a slight amplitude modulation grows exponentially because the nonlinearity shifts the local frequency, creating phase focusing that concentrates energy. The initially uniform wavetrain disintegrates into a series of wave groups — localized packets of high amplitude separated by nearly calm water. In oceanography, this is the mechanism behind rogue waves: apparently spontaneous extreme waves that emerge from a nearly uniform ocean swell.

The growth rate of the Benjamin–Feir instability for perturbation wavenumber \(K\) is:

\[ \sigma = \frac{\omega''}{2}\sqrt{K^2\!\left(\frac{4\gamma A_0^2}{\omega''} - K^2\right)} \]

Maximum growth occurs at \(K = K_{max} = \sqrt{2\gamma A_0^2/\omega''}\). Perturbations with \(K > K_{max}\) are stable; the instability band is \(0 < K < 2K_{max}\).

10.6 The Korteweg–de Vries Equation and Solitons

A different balance — between nonlinear steepening and weak dispersion — produces the Korteweg–de Vries (KdV) equation. Consider long, weakly nonlinear shallow-water waves. The full nonlinear shallow-water equations, expanded in the wave amplitude \(\epsilon\) and the longwave parameter \(\mu = (kH)^2\), give at leading order (the Boussinesq regime \(\epsilon \sim \mu \ll 1\)):

\[ \frac{\partial\eta}{\partial t} + c_0\frac{\partial\eta}{\partial x} + \frac{3c_0}{2H}\eta\frac{\partial\eta}{\partial x} + \frac{c_0 H^2}{6}\frac{\partial^3\eta}{\partial x^3} = 0 \]

In a frame moving at \(c_0 = \sqrt{gH}\), this simplifies to the KdV equation:

\[ \frac{\partial\eta}{\partial t} + \frac{3c_0}{2H}\eta\frac{\partial\eta}{\partial x} + \frac{c_0 H^2}{6}\frac{\partial^3\eta}{\partial x^3} = 0 \]

The second term is nonlinear steepening; the third is dispersive spreading from the third-derivative term.

The KdV equation admits an exact solitary wave (soliton) solution:

\[ \eta(x,t) = a\, \text{sech}^2\!\left(\sqrt{\frac{3a}{4H^3}}\,(x - ct)\right), \qquad c = c_0\!\left(1 + \frac{a}{2H}\right) \]

This is a localised wave of elevation, propagating without change of shape at speed \(c > c_0\). Crucially:

• Taller solitons are narrower and faster
• Solitons interact elastically: two solitons collide, interpenetrate, and emerge with their original shapes — only a phase shift remains

This elastic interaction was the first observation (by John Scott Russell in 1834, theorised by Boussinesq and finally by Korteweg and de Vries in 1895) that suggested waves could behave like particles. The discovery in the 1960s that the KdV equation is completely integrable — solvable by the inverse scattering transform — sparked the modern theory of integrable systems and revealed an infinite hierarchy of conserved quantities.

10.7 Inverse Scattering and Complete Integrability

The inverse scattering transform (IST, developed by Gardner, Greene, Kruskal, and Miura in 1967) is to nonlinear wave equations what the Fourier transform is to linear PDEs. The idea: map the initial condition \(\eta(x,0)\) to the scattering data of an associated linear Schrödinger equation (the direct problem), evolve the scattering data forward in time using simple linear equations, then reconstruct \(\eta(x,t)\) from the evolved scattering data (the inverse problem).

The scattering data consists of:

• Bound states (discrete eigenvalues): each corresponds to a soliton. There are exactly \(N\) solitons in the solution if the initial condition admits \(N\) bound states.
• Radiation (continuous spectrum): dispersive waves that spread and decay.

For any initial condition, the KdV equation evolves toward a superposition of solitons plus a decaying dispersive tail — the solitons emerge, separate by speed, and propagate to infinity. This gives a complete, nonlinear analogue of normal-mode decomposition.

10.8 WKB Theory: Dispersive Waves in Inhomogeneous Media

When the medium through which a wave propagates varies slowly in space — stratification changing with depth, water shoaling toward a beach, a slowly varying background current — the wave cannot be a pure plane wave with fixed wavenumber and frequency. Yet if the scale of variation \(L\) of the medium is much larger than the wavelength \(\lambda\), something close to a locally plane wave holds at each point.

The WKB method (after Wentzel, Kramers, and Brillouin, though the technique was known to Liouville and Green a century earlier) provides a systematic expansion in the small parameter \(\epsilon = \lambda/L \ll 1\). We seek solutions of the form:

\[ u(\mathbf{x}, t) = A(\mathbf{x}, t)\, e^{i\theta(\mathbf{x}, t)/\epsilon} + \text{c.c.} \]

where \(A\) is a slowly varying complex amplitude and \(\theta\) a rapidly varying phase. The local wavenumber and frequency are \(\mathbf{k} = \nabla\theta\) and \(\omega = -\partial\theta/\partial t\). These satisfy the compatibility condition \(\partial\mathbf{k}/\partial t + \nabla\omega = 0\), and together with the local dispersion relation \(\omega = \Omega(\mathbf{k}, \mathbf{x}, t)\) (now dependent on the slowly varying medium), they yield the ray equations:

\[ \frac{d\mathbf{x}}{dt} = \frac{\partial\Omega}{\partial\mathbf{k}} = \mathbf{c}_g, \qquad \frac{d\mathbf{k}}{dt} = -\frac{\partial\Omega}{\partial\mathbf{x}} \]

These are Hamilton’s equations in disguise: wavenumber plays the role of momentum, position the role of coordinates. Wave packets travel along rays — the characteristics of the kinematic equations — at the local group velocity.

For shoaling ocean waves approaching a beach, the depth \(H(x)\) decreases shoreward, so \(c_g = \sqrt{gH}\) decreases, the wavenumber \(k\) increases (wavelengths shorten), and amplitude grows. Whitham’s (1965) wave action conservation gives \(\mathcal{A} = E/\omega = A^2/(2\omega) = \text{const}\) along a ray:

\[ \frac{\partial\mathcal{A}}{\partial t} + \nabla\cdot(\mathbf{c}_g \mathcal{A}) = 0 \]

Wave action (not wave energy) is the adiabatic invariant conserved in a slowly varying medium — the direct analogue of the classical mechanical adiabatic invariant \(J = \oint p\, dq\). Energy is not conserved when the medium varies in time; action is.

For internal waves in a stratified ocean, rays refract as \(N(z)\) varies with depth, channelling wave energy into narrow beams and driving diapycnal mixing. At a caustic — a surface where rays converge and the WKB amplitude formally diverges — an Airy-function analysis replaces the leading-order WKB, capturing the wave tunnelling and the structure of the focal region. Caustics are the generic focusing surfaces of ray theory, ubiquitous in acoustics, optics, and geophysical wave propagation.

10.9 Nonlinear Resonant Interactions

When two or more waves coexist in a dispersive medium, they interact through the nonlinear terms in the governing equations. Most interactions are non-resonant: they produce oscillations at combination frequencies that quickly average to zero. For special combinations satisfying resonance conditions, however, the interaction is sustained and leads to significant energy exchange.

Three-wave resonance requires simultaneously:

\[ \mathbf{k}_1 + \mathbf{k}_2 = \mathbf{k}_3, \qquad \omega(\mathbf{k}_1) + \omega(\mathbf{k}_2) = \omega(\mathbf{k}_3) \]

When both conditions hold, a third wave grows secularly from the quadratic interaction of the first two, or a primary wave decays into two daughter waves (parametric decay instability). The triad amplitude equations take the universal form:

\[ \dot{A}_1 = \gamma A_2^* A_3, \qquad \dot{A}_2 = \gamma A_1^* A_3, \qquad \dot{A}_3 = -\gamma A_1 A_2 \]

with conserved Manley–Rowe relations \(|A_1|^2 + |A_3|^2 = \text{const}\) and \(|A_2|^2 + |A_3|^2 = \text{const}\). Whether three-wave resonance is possible depends on the geometry of the dispersion surface. For deep-water gravity waves (\(\omega = \sqrt{gk}\)), the dispersion curve is convex and three-wave resonance among waves of the same type is impossible — four-wave resonance dominates instead.

Four-wave resonance requires:

\[ \mathbf{k}_1 + \mathbf{k}_2 = \mathbf{k}_3 + \mathbf{k}_4, \qquad \omega_1 + \omega_2 = \omega_3 + \omega_4 \]

For deep-water gravity waves, these resonant quartets drive the slow evolution of the ocean wave spectrum (Hasselmann 1962): energy is transferred from the spectral peak toward both lower and higher frequencies. This is the dominant nonlinear process in operational wave forecasting models.

Wave turbulence (weak turbulence theory) treats the statistics of many resonantly interacting waves with random phases. Zakharov (1965–1968) found that the kinetic equations for the wave action spectrum \(n(\mathbf{k},t)\) admit exact stationary power-law solutions — the Kolmogorov–Zakharov spectra — analogous to Kolmogorov’s \(k^{-5/3}\) spectrum in strong turbulence but derived from first principles. These spectra describe the steady state in which wave action cascades through wavenumber space via resonant interactions, from the forcing scale to the dissipation scale.

10.10 Weakly Nonlocal Solitary Waves and Beyond-All-Orders Asymptotics

Classical solitons — the \(\text{sech}^2\) solutions of the KdV equation — are spatially localised and propagate without radiation. This perfection is special to completely integrable equations. For nearby but non-integrable equations — the fifth-order KdV, capillary-gravity water waves near the Bond-number resonance \(Bo = 1/3\) — true solitons do not exist. Instead, the system supports weakly nonlocal solitary waves: structures with a soliton-like core but small-amplitude oscillatory tails extending to infinity.

The canonical example is the fifth-order KdV equation:

\[ u_t + uu_x + u_{xxx} + \epsilon^2 u_{xxxxx} = 0 \]

The classical KdV soliton has phase velocity within the linear wave spectrum — linear waves of the same speed exist — so the soliton continuously radiates into these waves. The radiating tail has amplitude:

\[ A_{\text{tail}} \sim C\, e^{-\pi/(2\epsilon)} \]

This contains \(e^{-1/\epsilon}\), which vanishes faster than any power of \(\epsilon\). Such terms are beyond all orders of the power series in \(\epsilon\): they are invisible at every finite order of perturbation theory. This is why decades of perturbation calculations missed them entirely.

Beyond-all-orders asymptotics (exponential asymptotics) is the framework for computing these effects. The formal power series in \(\epsilon\) is divergent: the \(n\)-th coefficient grows like \(n!\). Optimal truncation — stopping at the smallest term — gives an exponentially accurate approximation. The remainder at optimal truncation is precisely the exponentially small tail, and can be computed by tracking the Stokes lines in the complex \(x\)-plane — curves across which the exponentially small contribution switches on (the classical Stokes phenomenon). Evaluating this switching gives the prefactor \(C\) explicitly (Yang and Akylas 1996; Boyd 1990).

Beyond solitary waves, exponential asymptotics appears in: quantum tunnelling rates (imaginary-time paths in Feynman integrals); Stokes phenomenon in the Airy equation and its generalisations; resurgence in quantum field theory, where non-perturbative effects (instantons) are related to the large-order growth of perturbation series; and the theory of water waves near parametric forcing resonances.

---

Chapter 11: Hydrodynamic Stability — AMATH 863

This chapter extends the course into graduate territory, covering the topics of AMATH 863 (Hydrodynamic Stability and Turbulence). The central question: given a known laminar flow, is it stable to small perturbations?

11.1 The Stability Problem

Every exact solution of the Navier–Stokes equations represents a possible flow — but not every possible flow is stable and therefore observable in practice. A flow is stable if small perturbations applied to it decay in time; it is unstable if some perturbation grows. The transition from laminar to turbulent flow is the most practically important manifestation of hydrodynamic instability.

The general approach is linear stability analysis:

1. Take a known base flow \(\overline{\mathbf{u}}(\mathbf{x})\)
2. Perturb it: \(\mathbf{u} = \overline{\mathbf{u}} + \varepsilon\mathbf{u}'(\mathbf{x},t)\) with \(|\varepsilon| \ll 1\)
3. Linearise the Navier–Stokes equations in \(\varepsilon\)
4. Seek modal solutions \(\mathbf{u}' \propto e^{ikx - i\omega t}\) (normal modes)
5. Determine whether \(\text{Im}(\omega) > 0\) (growing) or \(\text{Im}(\omega) < 0\) (decaying)

A flow is linearly stable if all normal modes decay (\(\text{Im}(\omega) < 0\) for all \(k\)). It is linearly unstable if any mode grows. Linear theory describes the onset and initial growth of instability; nonlinear theory is needed to describe the saturated state.

11.2 Normal Mode Analysis — The Orr–Sommerfeld Equation

Consider a 2D parallel shear flow \(\overline{\mathbf{u}} = (U(y), 0)\) — plane Poiseuille flow, plane Couette flow, a boundary layer profile. Introducing the stream function perturbation \(\psi' = \phi(y)e^{i(kx - \omega t)}\), the linearised Navier–Stokes equations reduce to a single ODE for the complex amplitude \(\phi(y)\):

This fourth-order ODE with complex coefficients requires four boundary conditions: \(\phi = \phi' = 0\) at each wall (no-slip and no-penetration). For a given \(k\) and \(Re\), the eigenvalue \(c\) is found. The flow is unstable if \(\text{Im}(c) > 0\) for any \(k\).

The Orr–Sommerfeld equation is the central object of viscous stability theory. In the inviscid limit \(Re \to \infty\), it reduces to the Rayleigh equation:

\[ (U-c)(\phi'' - k^2\phi) - U''\phi = 0 \]

11.3 Rayleigh’s Inflection Point Criterion

The proof proceeds by multiplying the Rayleigh equation by \(\phi^* / (U-c)^*\), integrating, and taking the imaginary part. The existence of an inflection point is not sufficient — a sharper condition is due to Fjørtoft (1950) — but its absence is sufficient for stability.

Physical interpretation: inflection points are locations of maximum shear, where vorticity gradients can drive an exponential Kelvin–Helmholtz-like instability. Profiles without inflection points — parabolic Poiseuille flow, Couette flow — are inviscidly stable. Yet they transition to turbulence in practice! Viscosity, counterintuitively, can be destabilising for these flows, as it modifies the phase relationship between the perturbation velocity and the vorticity, allowing the Orr–Sommerfeld modes to extract energy from the mean flow.

11.4 Kelvin–Helmholtz Instability

The Kelvin–Helmholtz (KH) instability arises at the interface between two fluids moving at different velocities — a shear layer. The idealised problem: fluid of density \(\rho_1\) moving at velocity \(U_1\) above fluid of density \(\rho_2\) moving at velocity \(U_2\), separated by a flat interface. Linearising the inviscid equations and matching boundary conditions at the interface:

The flow is unstable (\(\text{Im}(\omega) \neq 0\)) when the expression under the square root is negative:

\[ k > k_c = \frac{g(\rho_1^2 - \rho_2^2)}{\rho_1\rho_2(U_1-U_2)^2} \]

For \(\rho_1 = \rho_2\) (no density difference): always unstable for any non-zero velocity difference, at all wavenumbers. Surface tension modifies this by stabilising short waves.

The KH instability is ubiquitous: it is responsible for the formation of ocean surface waves by wind, for billows at the tropopause (clear-air turbulence), for the generation of vortices in mixing layers, and for instabilities in astrophysical jets. The beautiful billowing cat’s-eye patterns visible in clouds and in laboratory stratified shear layers are the nonlinear manifestation of this linear instability.

11.5 Rayleigh–Taylor Instability

The Rayleigh–Taylor (RT) instability arises when a denser fluid overlies a lighter fluid under gravity — the classic “heavy fluid over light fluid” problem. Even in the absence of mean flow, the configuration is unstable.

For two semi-infinite, inviscid, irrotational layers with densities \(\rho_2 > \rho_1\) (heavy over light), with gravity \(g\) acting downward:

\[ \omega^2 = gk\frac{\rho_1 - \rho_2}{\rho_1 + \rho_2} < 0 \quad \text{(since } \rho_2 > \rho_1\text{)} \]

So \(\omega\) is purely imaginary, with growth rate:

\[ \sigma = \text{Im}(\omega) = \sqrt{gk\frac{\rho_2 - \rho_1}{\rho_1 + \rho_2}} \]

Short waves grow faster: the RT instability is ultraviolet (high-\(k\)) unstable. Surface tension \(\sigma_T\) stabilises short waves with \(k > k_c = \sqrt{g(\rho_2-\rho_1)/\sigma_T}\), but long waves (\(k < k_c\)) remain unstable. The nonlinear evolution produces the spectacular mushroom clouds and fingers characteristic of this instability, seen in inertial confinement fusion implosions, supernova remnants, and ink dropping into water.

11.6 Rayleigh–Bénard Convection

Rayleigh–Bénard convection is the instability of a fluid layer heated uniformly from below. A horizontal layer of fluid of depth \(d\) is maintained at temperature \(T_0 + \Delta T\) at the bottom and \(T_0\) at the top. The density is lower at the bottom (since hot fluid expands), so this is a potentially unstable density stratification — but viscosity and thermal diffusion work against it.

The competition is characterised by the Rayleigh number:

\[ Ra = \frac{g\alpha\Delta T d^3}{\nu\kappa} \]

where \(\alpha\) is the thermal expansion coefficient, \(\kappa\) the thermal diffusivity, \(\nu\) the kinematic viscosity. The Rayleigh number is the ratio of buoyancy driving to diffusive damping.

The onset occurs at a critical wavenumber \(k_c\), corresponding to convection cells with a specific aspect ratio. This is a beautiful example of a bifurcation: a qualitative change in the nature of the steady state as a parameter (\(Ra\)) crosses a threshold.

Beyond \(Ra_c\), convection rolls are steady; at higher \(Ra\), they undergo further bifurcations — oscillatory convection, then chaotic convection, and ultimately turbulent convection. Rayleigh–Bénard turbulence at very high \(Ra\) is the model for convection in the Earth’s mantle, the Sun’s convection zone, and the atmospheres of giant planets.

11.7 Energy Methods and Nonlinear Stability

Linear stability theory determines whether infinitesimal perturbations grow. But it says nothing about the stability of flows to finite-amplitude perturbations. A flow can be linearly stable yet jump to turbulence when perturbed sufficiently strongly — this is precisely what happens in pipe flow (Hagen–Poiseuille), which is linearly stable for all \(Re\) yet turbulent in practice above \(Re \approx 2300\).

The energy method (Serrin 1959, Joseph 1976) provides rigorous nonlinear stability results. Define the perturbation energy:

\[ E(t) = \frac{1}{2}\int_\mathcal{V} |\mathbf{u}'|^2\, dV \]

Taking the time derivative and using the Navier–Stokes equations for the perturbation:

\[ \frac{dE}{dt} = -\int_\mathcal{V} u'_i u'_j\frac{\partial\overline{u}_i}{\partial x_j}\, dV - \nu\int_\mathcal{V} |\nabla\mathbf{u}'|^2\, dV \]

The first integral is energy production from the mean flow shear; the second is viscous dissipation. If the ratio of production to dissipation is bounded above by 1, then \(dE/dt \leq 0\) and the perturbation must decay — regardless of its amplitude.

This yields a global stability threshold \(Re_E\) (the energy Reynolds number): for \(Re < Re_E\), all perturbations, however large, decay monotonically. For pipe flow, \(Re_E \approx 81.5\). The gap between \(Re_E = 81.5\) and the experimental transition threshold \(Re \approx 2300\) is a region where small perturbations decay but large ones can sustain themselves — the realm of subcritical transition and turbulent puffs, which remains an active research area.

11.8 Centrifugal Instability — Taylor–Couette Flow

Centrifugal instability arises in rotating flows when the centrifugal force on a radially displaced fluid ring exceeds the restoring pressure gradient, so the ring continues outward rather than returning. It is the rotational analogue of Rayleigh–Taylor instability, and it produces the most visually stunning sequence of bifurcations in fluid mechanics.

The canonical example is Taylor–Couette flow: viscous fluid between two coaxial cylinders of radii \(R_1 < R_2\), rotating at angular velocities \(\Omega_1\) and \(\Omega_2\). The base flow is an azimuthal profile:

\[ V(r) = Ar + \frac{B}{r}, \qquad A = \frac{\Omega_2 R_2^2 - \Omega_1 R_1^2}{R_2^2 - R_1^2}, \quad B = \frac{(\Omega_1 - \Omega_2)R_1^2 R_2^2}{R_2^2 - R_1^2} \]

The inviscid stability criterion was derived by Lord Rayleigh (1916) using an elegant physical argument. Consider a fluid ring at radius \(r_0\) with angular momentum \(L_0 = r_0 V(r_0)\). Displaced to \(r_0 + \delta r\) with angular momentum conserved (inviscid), the ring’s velocity becomes \(V' = L_0/(r_0 + \delta r)\). The centrifugal acceleration of the displaced ring exceeds that of the ambient fluid if and only if \(L(r)^2 = (rV)^2\) decreases outward — in which case the ring feels a net outward force and the flow is unstable.

With viscosity, the critical threshold is quantified by the Taylor number:

\[ Ta = \frac{\Omega_1^2 R_1^2 (R_2 - R_1)^3}{\nu^2 R_2} \]

G.I. Taylor (1923) solved the viscous stability problem exactly and found instability at \(Ta_c \approx 1708\), predicting that the flow should develop Taylor vortices — a stack of toroidal counter-rotating cells stacked axially. Taylor famously photographed these vortices, providing one of the first quantitative agreements between linear stability theory and experiment. The match was precise to within a few percent: a landmark moment in the history of fluid mechanics.

Beyond the first bifurcation, Taylor–Couette flow undergoes a remarkably rich sequence as \(Ta\) increases: steady Taylor vortices → wavy vortices (the azimuthal symmetry breaks) → modulated wavy vortices → turbulent Taylor vortices → featureless turbulence. This cascade has been exhaustively studied as a paradigm for routes to turbulence, and each transition has been identified with a specific bifurcation type (pitchfork, Hopf, torus bifurcation). The system remains an important experimental and theoretical benchmark.

11.9 Barotropic Instability

Barotropic instability is the two-dimensional analogue of Kelvin–Helmholtz instability in a rotating geophysical fluid. It occurs in flows with no vertical structure (barotropic flows — pressure surfaces are horizontal and coincide with density surfaces) when perturbations can extract energy from the horizontal shear of the mean flow. It is responsible for the break-up of atmospheric jets and oceanic currents into synoptic-scale eddies.

On the \(\beta\)-plane, the linearised quasi-geostrophic potential vorticity equation for stream function perturbations \(\psi'\) to a zonal mean flow \(\overline{U}(y)\) is the Rayleigh–Kuo equation:

\[ (\overline{U} - c)\!\left(\frac{\partial^2\psi'}{\partial y^2} - k^2\psi'\right) + \left(\beta - \frac{\partial^2\overline{U}}{\partial y^2}\right)\psi' = 0 \]

This generalises the Rayleigh equation of §11.3 by replacing the vorticity gradient \(-U''\) with the absolute vorticity gradient \(\beta - U''\). The necessary condition for instability generalises accordingly:

When \(q_y = 0\) somewhere, Rossby waves propagating on the vorticity gradient can be over-reflected, extracting energy from the mean flow. The condition \(q_y = 0\) requires \(U_{yy} = \beta > 0\) — a region of strong positive curvature in the jet. Atmospheric jets (particularly the polar jet stream) can satisfy this condition at their flanks, generating barotropically unstable modes that grow into synoptic-scale eddies. In the ocean, the Gulf Stream and Kuroshio extension are subject to barotropic instability, producing the energetic meanders and rings observed by satellite.

11.10 Baroclinic Instability — The Eady Model

Baroclinic instability is the dominant instability mechanism in mid-latitude weather systems and one of the most important dynamical processes in geophysical fluid dynamics. Unlike barotropic instability, which involves purely horizontal shear of a vertically uniform flow, baroclinic instability operates in stratified, rotating fluids where the density surfaces (isopycnals) are tilted — as they are in the atmosphere, where temperature decreases poleward, and in the ocean, where salinity and temperature both create lateral density gradients.

The tilted isopycnal structure stores available potential energy — the potential energy that would be released if the density surfaces were flattened to horizontal. Baroclinic instability converts this potential energy into kinetic energy of growing eddies. The physical process is best visualised through the concept of slantwise convection: a fluid parcel can move along a slope that cuts across isopycnals, converting potential energy to kinetic energy, if the slope of the parcel path is shallower than the slope of the isopycnals. This is a weaker condition than simple vertical convection and is satisfied over vast regions of the mid-latitude atmosphere and ocean.

The Eady Model

The simplest and most influential model of baroclinic instability is due to E.T. Eady (1949). The setup: a Boussinesq fluid between rigid horizontal lids at \(z = 0\) and \(z = H\), with uniform Coriolis parameter \(f\) (no \(\beta\)) and uniform buoyancy frequency \(N\). The mean flow is a vertical shear \(\overline{U}(z) = \Lambda z\) — a “thermal wind” balanced by a horizontal temperature gradient via the thermal wind relation \(f\partial U/\partial z = -(g/\rho_0)\partial\rho/\partial y\). The perturbation quasi-geostrophic potential vorticity equation in the interior is satisfied trivially (the interior QGPV gradient vanishes), so the dynamics are entirely controlled by temperature anomalies at the two boundaries.

Seeking normal mode solutions \(\psi' \propto e^{ik(x-ct)}\), the dispersion relation is:

The maximum growth rate — achieved at \(\mu \approx 1.61\) — is:

\[ \sigma_{max} = \text{Im}(kc)_{max} \approx 0.31\frac{f\Lambda}{N} \]

The corresponding most unstable wavelength is \(L^* \sim NH/f = L_R\) — the Rossby deformation radius. For the atmosphere, \(L_R \approx 1000\) km, consistent with the observed scale of mid-latitude cyclones. The e-folding growth time \(\sigma_{max}^{-1} \sim 1\) day matches the observed rapid intensification of extratropical storms.

Mechanism: Coupled Edge Waves

The physical mechanism of Eady instability is beautifully transparent in terms of boundary Rossby waves. Each horizontal boundary supports a temperature (Rossby) edge wave: a wave whose amplitude decays away from the boundary, propagating along the boundary at a speed set by the mean temperature gradient there. In isolation, the lower-boundary wave propagates eastward and the upper-boundary wave propagates westward relative to the mean flow at each level.

For a specific vertical tilt between the two edge waves, they phase-lock: their relative propagation speeds cancel and they remain stationary relative to each other. When phase-locked with the correct tilt (upper wave lagging the lower by a quarter wavelength), the temperature anomaly of the upper wave is maximally correlated with the vertical velocity induced by the lower wave — the upper wave amplifies the lower, and vice versa. The result is mutual amplification: an exponentially growing mode. The Eady instability is thus a constructive interference between two vertically separated Rossby waves, held together against their natural tendency to drift apart by the shear of the mean flow.

The Phillips Model and the Charney–Stern Condition

The Phillips (1954) two-layer model retains the \(\beta\)-effect, representing the atmosphere as two fluid layers each of uniform density. The stability analysis yields a generalised Rayleigh–Kuo condition, and the most important theoretical result is the:

This result unifies barotropic instability (sign change of \(q_y = \beta - U_{yy}\) in the interior) and baroclinic instability (sign change at the boundaries — Eady model) into a single framework. It is the deepest and most general necessary condition for instability in quasi-geostrophic flows.

11.11 Homogeneous Isotropic Turbulence

The AMATH 863 curriculum treats the full spectrum of turbulence, from onset to the statistical theory of fully developed turbulence. The idealised case of homogeneous isotropic turbulence (HIT) — statistically uniform in all positions and directions — provides the theoretical foundation, even though it is never exactly realised in geophysical flows.

Kolmogorov’s 1941 Theory

Kolmogorov argued that at sufficiently high Reynolds number, there exists an inertial subrange of scales — larger than the dissipation scale but smaller than the energy-injection scale — in which neither forcing nor viscosity plays a direct role. In this range, the only relevant dimensional quantity is the energy dissipation rate \(\varepsilon\) [m²/s³]. Dimensional analysis then uniquely determines the energy spectrum:

\[ E(k) = C_K\, \varepsilon^{2/3} k^{-5/3} \]

where \(C_K \approx 1.5\) is the Kolmogorov constant. This \(-5/3\) power law is one of the most celebrated results in theoretical physics: confirmed over more than five decades in wavenumber in atmospheric boundary layer data, oceanic measurements, wind tunnel experiments, and DNS. The Kolmogorov microscale — the scale at which viscous dissipation dominates — is:

\[ \eta = \!\left(\frac{\nu^3}{\varepsilon}\right)^{1/4} \]

The ratio of energy-containing scale \(L\) to dissipation scale \(\eta\) grows as \(L/\eta \sim Re^{3/4}\), so the number of degrees of freedom in 3D turbulence scales as \(Re^{9/4}\). For ocean currents (\(Re \sim 10^{10}\)), resolving all scales requires \(\sim 10^{22.5}\) grid points — far beyond any foreseeable computer.

Turbulence Closure

The closure problem (§7.2 of AMATH 463) is the central unsolved problem of turbulence theory: the RANS equations for the mean flow involve the Reynolds stress tensor \(\overline{u'_i u'_j}\), for which no exact equation in terms of mean quantities exists. Modern approaches:

• \(k\)-\(\varepsilon\) models: solve transport equations for TKE \(k\) and dissipation rate \(\varepsilon\), relating Reynolds stresses to \(k\), \(\varepsilon\), and mean strain via the Boussinesq hypothesis \(\overline{u'_i u'_j} = -2\nu_T S_{ij} + \frac{2}{3}k\delta_{ij}\), where \(\nu_T \propto k^2/\varepsilon\) is the turbulent viscosity
• Large Eddy Simulation (LES): resolve scales larger than a filter width \(\Delta\), parameterise only subgrid scales; feasible for engineering Reynolds numbers
• Direct Numerical Simulation (DNS): resolve all scales from \(L\) to \(\eta\); computationally exact but limited to \(Re \lesssim 10^4\)

Stratified Turbulence and the Miles–Howard Criterion

Stratification fundamentally alters turbulence. The key parameter is the Richardson number:

\[ Ri = \frac{N^2}{(\partial U/\partial z)^2} \]

The proof — due to Miles (1961) with the Howard semicircle theorem providing the geometric picture — follows by multiplying the Taylor–Goldstein equation (the stratified analogue of the Rayleigh equation) by an appropriate function and integrating. When \(Ri < 1/4\), the kinetic energy of the shear can overcome the work done against buoyancy in vertically displacing fluid, allowing instability. The resulting turbulent mixing (Kelvin–Helmholtz billows in stratified flow) is the primary mechanism driving irreversible diapycnal mixing in the interior of the ocean — essential for the global thermohaline circulation.

At \(Ri \gg 1\), stratification strongly suppresses vertical motions. Turbulence becomes quasi-two-dimensional, with energy cascading to larger horizontal scales (an inverse energy cascade, as in 2D turbulence) rather than to smaller scales as in 3D. This two-dimensionalisation by stratification underlies the large horizontal coherence of oceanic mesoscale eddies and atmospheric jets, providing the dynamical bridge between the instabilities studied throughout this chapter and the organised large-scale structures observed in the geophysical world.

Part II — Engineering Applications of Fluid Mechanics

The preceding parts developed fluid mechanics as a branch of mathematical physics: continuum kinematics, the Navier–Stokes equations as conservation laws for mass and momentum, boundary-layer asymptotics, hydrodynamic stability, geophysical dynamics, and the statistical mechanics of turbulence. That viewpoint is essential, but it is not the only one an engineer needs. When the task is to size a pump, choose a pipe diameter, design a weir, or verify that a scaled wind-tunnel model faithfully reproduces full-scale aerodynamics, the relevant questions are rarely about existence of weak solutions or spectral properties of the Orr–Sommerfeld operator — they are about dimensionless groups, empirical correlations, and similitude.

This part surveys the applied toolkit taught in introductory engineering fluid-mechanics courses at the University of Waterloo (ME 351, CIVE 280, CHE 211, SYDE 383, MTE 352, ENVE 280, AE 280, GEOE 280). The material complements — rather than replaces — the rigorous continuum derivations of Parts I–VII: the Hagen–Poiseuille profile reappears here, but as the basis of the Darcy–Weisbach friction factor; Bernoulli’s equation reappears, but as the working equation behind a pitot tube; boundary-layer theory reappears, but in the form of engineering drag-coefficient charts. Throughout, the style remains mathematical — theorems, definitions, derivations — but the emphasis shifts toward the correlations and design equations that engineers actually use.

Chapter 12: Engineering Dimensional Analysis

The Buckingham Pi Theorem

The \(\Pi\) groups are constructed by choosing \(k\) repeating variables whose dimensions span the fundamental set, and forming dimensionless products with each of the remaining variables in turn. The choice is non-unique, but any two complete sets of \(\Pi\)’s are related by invertible functions.

Standard Dimensionless Numbers

Each group expresses the ratio of two physical effects, and the flow regime is governed by which effect dominates:

• Reynolds number \(Re\) — inertia / viscous. Dominates in all viscous flows; controls transition to turbulence (Ch. 2, 5).
• Froude number \(Fr\) — inertia / gravity. Dominates in free-surface flows: ship waves, open channels, hydraulic jumps.
• Mach number \(Ma\) — flow speed / sound speed. Dominates when \(Ma \gtrsim 0.3\); compressibility becomes essential.
• Weber number \(We\) — inertia / surface tension. Relevant for droplets, sprays, small-scale interfacial flows.
• Euler number \(Eu\) — pressure / inertia. The dimensionless pressure drop; closely related to the loss and drag coefficients of later chapters.
• Strouhal number \(St\) — characterises unsteady vortex shedding (e.g. \(St \approx 0.21\) for a cylinder at \(Re \sim 10^3\)).
• Prandtl number \(Pr\) — momentum diffusivity / thermal diffusivity. Couples fluid to heat transfer.

Similitude

Dynamic similarity is the strongest and what Pi-theorem-based model testing actually demands. Wind tunnels match \(Re\) (and \(Ma\) for high-speed aerodynamics); towing tanks match \(Fr\) (surface waves); open-channel hydraulic models match \(Fr\); scaled heat-exchanger rigs match \(Re\) and \(Pr\).

Worked Example: Ship Resistance Scaling

Total ship resistance splits into frictional drag (boundary-layer, scales with \(Re\)) and wave-making drag (free-surface gravity waves, scales with \(Fr\)). A 1:25 scaled model in a towing tank cannot satisfy both simultaneously with water: matching \(Fr\) requires model speed \(U_m = U_p/\sqrt{25} = U_p/5\); matching \(Re\) would require \(U_m = 25\, U_p\). The Froude–Reynolds conflict is resolved by Froude’s hypothesis:

\[ C_{T} = C_{F}(Re) + C_{W}(Fr), \]\[ C_F = \frac{0.075}{(\log_{10} Re - 2)^2}, \]

leaving the wave coefficient \(C_W\), which is scaled directly to full size. The prototype friction is then re-added at the (much higher) full-scale \(Re\). This methodology, due to William Froude in the 1870s, remains the foundation of ship-resistance testing today.

Chapter 13: Fluid Statics — Hydrostatics and Buoyancy

The Hydrostatic Pressure Distribution

In a fluid at rest the momentum equation reduces to \(\nabla p = \rho \mathbf{g}\). For a constant-density fluid in a uniform gravitational field \(\mathbf{g} = -g\hat{\mathbf{z}}\),

\[ p(z) = p_0 + \rho g (z_0 - z), \]

so pressure increases linearly with depth. For a compressible gas under gravity with ideal-gas equation of state and isothermal atmosphere \(T = T_0\), integrating \(dp/dz = -\rho g = -p g/(R T_0)\) gives the barometric formula \(p(z) = p_0 \exp(-z/H)\) with scale height \(H = R T_0/g \approx 8.4~\mathrm{km}\) for Earth’s atmosphere.

Force on a Plane Submerged Surface

A plane surface of area \(A\) submerged in a fluid with its centroid at depth \(\bar{h}\) experiences a resultant hydrostatic force

\[ F = \rho g\, \bar{h}\, A, \]

acting normal to the surface at the centre of pressure, a distance \(I_{\bar{x}\bar{x}}/(\bar{h} A)\) below the centroid (with \(I_{\bar{x}\bar{x}}\) the second moment of area about a horizontal axis through the centroid). For a curved surface one resolves into horizontal and vertical components: the horizontal force equals that on the vertical projection, and the vertical force equals the weight of the fluid column above (or virtually displaced by) the surface.

Archimedes’ Principle and Stability of Floating Bodies

The net pressure force on any closed surface immersed in a fluid at rest equals the weight of the displaced fluid:

\[ \mathbf{F}_{\text{buoy}} = -\int_{\partial V} p\, \hat{\mathbf{n}}\, dA = \rho_f g V\, \hat{\mathbf{z}}. \]

This is Archimedes’ principle. For a floating body with waterplane second moment \(I_{\mathrm{wp}}\) and submerged volume \(V_s\), the metacentric height is

\[ \overline{GM} = \frac{I_{\mathrm{wp}}}{V_s} - \overline{BG}, \]

where \(\overline{BG}\) is the distance from centre of buoyancy to centre of gravity. The body is stable in roll iff \(\overline{GM} > 0\) — the design criterion for ships, barges, and offshore platforms.

Remark. Hydrostatics underlies manometry, dam design, lock gates, and the calibration of depth gauges. The limit \(\mathbf{u} \equiv 0\) makes it a trivial special case of the Navier–Stokes equations, but engineering practice demands that it be stated explicitly, since pressure-on-submerged-surfaces and stability-of-floating-bodies problems are standard content in CIVE 280, ME 351, ENVE 280 and the equivalent first courses.

---

Chapter 14: Internal Flow — Pipe Flow and Friction

Regimes and Hagen–Poiseuille

For flow in a straight circular pipe of diameter \(D\), the transition Reynolds number \(Re_D = \rho U D/\mu\) is conventionally taken as \(Re_D \lesssim 2300\) (laminar), \(Re_D \gtrsim 4000\) (turbulent), with a transitional range between.

\[ \frac{1}{r}\frac{d}{dr}\!\left(r\frac{du_z}{dr}\right) = \frac{1}{\mu}\frac{dp}{dz}, \]

with boundary conditions \(u_z(R) = 0\) and \(du_z/dr|_{r=0} = 0\). Integrating yields the Hagen–Poiseuille profile:

\[ u_z(r) = \frac{R^2}{4\mu}\!\left(-\frac{dp}{dz}\right)\!\left(1 - \frac{r^2}{R^2}\right), \qquad Q = \frac{\pi R^4}{8\mu}\!\left(-\frac{dp}{dz}\right). \]

Darcy–Weisbach Equation

For engineering purposes, one packages the pressure drop over a length \(L\) of pipe as

\[ f = \frac{64}{Re_D}. \]

The Moody Diagram

For turbulent flow, \(f\) depends on both \(Re_D\) and the relative roughness \(\epsilon/D\). The standard correlations are:

• Blasius (smooth, \(4\times 10^3 \le Re_D \le 10^5\)): \(f = 0.316\, Re_D^{-1/4}\).
• Colebrook–White (implicit, all turbulent regimes): \[ \frac{1}{\sqrt{f}} = -2\log_{10}\!\left(\frac{\epsilon/D}{3.7} + \frac{2.51}{Re_D\sqrt{f}}\right). \]
• Swamee–Jain (explicit approximation, within 1% of Colebrook): \[ f = \frac{0.25}{\left[\log_{10}\!\left(\frac{\epsilon/D}{3.7} + \frac{5.74}{Re_D^{0.9}}\right)\right]^2}. \]
• Haaland (explicit, within 2%): \[ \frac{1}{\sqrt{f}} = -1.8\log_{10}\!\left[\left(\frac{\epsilon/D}{3.7}\right)^{1.11} + \frac{6.9}{Re_D}\right]. \]

The graphical representation of these correlations is the Moody diagram — a log-log plot of \(f\) versus \(Re_D\) with \(\epsilon/D\) as parameter. Three regimes are visible: the laminar line \(f=64/Re_D\); a transition region for smooth turbulent flow where \(f\) follows Blasius; and a fully rough regime in which \(f\) is independent of \(Re_D\) and determined by roughness alone, since the viscous sublayer becomes thinner than the roughness elements.

Minor (Local) Losses

\[ h_L = K\,\frac{U^2}{2g}. \]

Representative values: sharp-edged entrance \(K \approx 0.5\); well-rounded entrance \(K \approx 0.04\); sudden expansion from area \(A_1\) to \(A_2\), \(K = (1 - A_1/A_2)^2\) (Borda–Carnot); 90\(^\circ\) elbow \(K \approx 0.3\)–\(0.9\); fully open globe valve \(K \approx 10\); gate valve \(K \approx 0.15\).

Non-Newtonian Pipe Flow — Power-Law and Bingham Fluids

Polymer melts, slurries, blood, ketchup, drilling muds, and many food products fail the Newtonian assumption that shear stress is proportional to strain rate. Engineering practice in CHE 543 (polymer processing), CHE 564 (food processing) and BME 384 (biofluids) relies on two simple generalised-Newtonian models.

\[ \tau = K\, |\dot\gamma|^{n-1}\, \dot\gamma, \]\[ u(r) = \frac{n}{n+1}\left(\frac{\Delta p}{2 K L}\right)^{1/n}\!\!\left(R^{(n+1)/n} - r^{(n+1)/n}\right), \]\[ Q = \frac{n\pi R^3}{3n+1}\left(\frac{R\, \Delta p}{2 K L}\right)^{1/n}, \]

which reduces to Hagen–Poiseuille at \(n = 1\), \(K = \mu\). A generalised Reynolds number \(Re_{\text{gen}} = \rho U^{2-n} D^n / K\,[8(3n+1)/(4n)]^{1-n}\) extends the Moody diagram by matching the laminar \(64/Re\) line exactly (Metzner–Reed correlation).

\[ \tau = \tau_y + \mu_p\, \dot\gamma \quad\text{for } |\tau| > \tau_y,\qquad \dot\gamma = 0\ \text{otherwise}. \]

In pipe flow this creates an unsheared plug of radius \(r_p = 2 L \tau_y/\Delta p\) moving as a rigid body surrounded by a sheared annulus — the Buckingham–Reiner equation gives the flow rate. The dimensionless Bingham number \(Bn = \tau_y D/(\mu_p U)\) measures the relative importance of yield to viscous stresses and controls whether the plug fills the pipe (no flow) or vanishes (Newtonian limit). Herschel–Bulkley (\(\tau = \tau_y + K\dot\gamma^n\)) combines both features and is the standard food-engineering model for tomato paste, chocolate, and many dairy products.

---

Chapter 15: Pipe Network Analysis

Energy Equation with Losses

\[ \frac{p_1}{\rho g} + \frac{U_1^2}{2g} + z_1 + h_{\mathrm{pump}} = \frac{p_2}{\rho g} + \frac{U_2^2}{2g} + z_2 + h_{\mathrm{turbine}} + \sum h_f + \sum h_L, \]

the Bernoulli equation of Chapter 3 augmented by mechanical work and friction.

Series and Parallel Combinations

Hardy Cross Method

For a looped network with \(n_p\) pipes, \(n_j\) junctions, and \(n_\ell\) independent loops, continuity at each junction (\(\sum Q_\text{in} = \sum Q_\text{out}\)) and head-loss closure around each loop (\(\sum_{\text{loop}} h_f = 0\), counted with sign) give \(n_p\) equations for \(n_p\) unknown pipe flows.

The method is Newton’s method applied loop-by-loop to the nonlinear head-loss equations.

Worked Example: Y-Branch

\[ H_A - r_1 Q_1^2 - r_2 Q_2^2 = z_B, \qquad H_A - r_1 Q_1^2 - r_3 Q_3^2 = z_C. \]

Subtracting, \(r_2 Q_2^2 - r_3 Q_3^2 = z_C - z_B\); combined with continuity, a single nonlinear equation in \(Q_2\) remains, solvable by bisection or Newton iteration.

Chapter 16: Pumps and Turbomachinery Basics

Euler Turbomachine Equation

Applying the angular-momentum form of the Reynolds transport theorem (Chapter 2) to a steady control volume bounded by the inlet (1) and outlet (2) of a rotating impeller:

For a centrifugal pump with radial inflow (\(V_{t,1}=0\)), the ideal head is \(gH = U_2 V_{t,2}\).

Characteristic Curves and Operating Point

\[ H_{\mathrm{sys}}(Q) = \Delta z + \sum (f\tfrac{L}{D} + K)\frac{U^2}{2g} = \Delta z + CQ^2, \]

i.e. a static lift plus quadratic friction term. The operating point is the intersection \(H_{\mathrm{pump}}(Q) = H_{\mathrm{sys}}(Q)\).

Specific Speed and Selection

\(N_s\) is a similarity invariant that identifies pump type: low \(N_s\) (\(< 1\)) — radial/centrifugal; mid \(N_s\) (\(1\)–\(4\)) — mixed flow; high \(N_s\) (\(>4\)) — axial.

NPSH and Cavitation

\[ \mathrm{NPSH_a} = \frac{p_{\mathrm{atm}} - p_v}{\rho g} - z_{\mathrm{suction}} - h_{f,\mathrm{suction}} \]\[ \mathrm{NPSH_a} \geq \mathrm{NPSH_r} + \mathrm{margin}. \]

Pumps in Series and Parallel; Affinity Laws

For pumps in series at common \(Q\), heads add: \(H_{\mathrm{total}}(Q) = H_1(Q) + H_2(Q)\). For pumps in parallel at common \(H\), flows add: \(Q_{\mathrm{total}}(H) = Q_1(H) + Q_2(H)\).

These follow from dimensional analysis of the pump-performance dimensionless groups \(C_H = gH/(\omega D)^2\), \(C_Q = Q/(\omega D^3)\), \(C_P = P/(\rho \omega^3 D^5)\).

Chapter 17: Flow Measurement

Pitot–Static Tube

\[ p_0 = p + \tfrac{1}{2}\rho U^2, \]\[ U = \sqrt{\frac{2(p_0 - p)}{\rho}}. \]

For compressible flow at \(Ma > 0.3\), the isentropic-flow correction \(p_0/p = (1 + \tfrac{\gamma-1}{2}Ma^2)^{\gamma/(\gamma-1)}\) must be applied.

Differential-Pressure Flowmeters

Typical values: Venturi \(C_d \approx 0.98\); flow nozzle \(C_d \approx 0.96\); orifice plate \(C_d \approx 0.60\)–\(0.65\) (depending on \(\beta = d/D\) and tap location). Orifice plates are cheap but dissipative; Venturis are expensive but recover most of the pressure.

Weirs (Open-Channel)

\[ Q = C_d\,\tfrac{2}{3}b\sqrt{2g}\,H^{3/2}, \qquad C_d \approx 0.62. \]\[ Q = C_d\,\tfrac{8}{15}\tan(\theta/2)\sqrt{2g}\,H^{5/2}, \qquad C_d \approx 0.58. \]

The \(H^{5/2}\) dependence of a V-notch gives excellent low-flow sensitivity, which is why V-notches are standard in environmental monitoring.

Electromagnetic and Ultrasonic Meters

An electromagnetic flowmeter applies Faraday’s law: a conductive fluid moving at velocity \(U\) through a transverse magnetic field \(B\) induces a voltage \(V = BUD\) across electrodes separated by the pipe diameter \(D\). Ultrasonic transit-time meters measure the difference in travel time of acoustic pulses sent upstream and downstream along a diagonal path; the difference is proportional to the line-averaged axial velocity.

Uncertainty Propagation

\[ \sigma_Q^2 = \sum_{i} \!\left(\frac{\partial Q}{\partial x_i}\right)^{\!2}\sigma_{x_i}^2, \]

assuming uncorrelated errors. For an orifice meter, \(Q \propto C_d\sqrt{\Delta p}\) implies \((\sigma_Q/Q)^2 = (\sigma_{C_d}/C_d)^2 + \tfrac14(\sigma_{\Delta p}/\Delta p)^2\).

Chapter 18: Open-Channel Hydraulics

Chezy and Manning Equations

For uniform flow in a prismatic channel of slope \(S_0\), hydraulic radius \(R_h = A/P\) (area over wetted perimeter), a momentum balance between gravity and bed friction yields

Specific Energy and Critical Flow

\[ y_c = \!\left(\frac{q^2}{g}\right)^{\!1/3}, \qquad E_{\min} = \tfrac{3}{2}y_c. \]

The \(E(y)\) curve for fixed \(q\) has two branches — the upper (subcritical, large \(y\)) and lower (supercritical, small \(y\)) — meeting at \((y_c, E_{\min})\).

Hydraulic Jump

\[ \rho q(U_2 - U_1) = \tfrac{1}{2}\rho g(y_1^2 - y_2^2), \]\[ \frac{y_2}{y_1} = \tfrac{1}{2}\!\left(\sqrt{1 + 8 Fr_1^2} - 1\right). \]\[ \Delta E = \frac{(y_2 - y_1)^3}{4 y_1 y_2}. \]

Jumps are used deliberately downstream of spillways as energy dissipators.

Gradually Varied Flow

\[ \frac{dy}{dx} = \frac{S_0 - S_f}{1 - Fr^2}, \]

where \(S_f\) is the friction slope (from Manning). The sign of the numerator and denominator, relative to the normal depth \(y_n\) (where \(S_0 = S_f\)) and critical depth \(y_c\), classifies the water-surface profiles. For Mild slopes (\(y_n > y_c\)): M1 (\(y > y_n\)), M2 (\(y_c < y < y_n\)), M3 (\(y < y_c\)). Analogous classifications exist for Steep (\(y_n < y_c\)), Critical (\(y_n = y_c\)), Horizontal (\(S_0 = 0\), no \(y_n\)), and Adverse (\(S_0 < 0\)) slopes.

Chapter 19: External Flow — Lift, Drag, and Boundary-Layer Engineering Data

Drag Coefficient for Standard Shapes

The \(C_D(Re)\) curve for a smooth sphere exhibits several canonical regimes:

• Stokes regime, \(Re \ll 1\): \(C_D = 24/Re\) (exact, from the creeping-flow Stokes solution).
• Intermediate (Schiller–Naumann), \(Re \lesssim 10^3\): \(C_D \approx \tfrac{24}{Re}(1 + 0.15 Re^{0.687})\).
• Newton regime, \(10^3 \lesssim Re \lesssim 2\times 10^5\): \(C_D \approx 0.44\), nearly constant.
• Drag crisis, \(Re \sim 3\times 10^5\): boundary layer transitions from laminar to turbulent on the sphere; separation moves downstream; \(C_D\) drops sharply from \(\sim 0.47\) to \(\sim 0.1\). Roughness (dimples on a golf ball) triggers this transition earlier.

For a circular cylinder the picture is qualitatively similar, with an additional feature: in the range \(50 \lesssim Re \lesssim 2\times 10^5\), a Kármán vortex street is shed with Strouhal number \(St \approx 0.21\). Resonance with structural modes drives aeroelastic phenomena (Tacoma Narrows, marine risers, chimney stacks).

For a flat plate at zero incidence, length \(L\), the laminar-Blasius result (Chapter 4) gives \(C_D = 1.328/\sqrt{Re_L}\), and turbulent-boundary-layer correlations give \(C_D \approx 0.074/Re_L^{1/5}\) (Prandtl, \(5\times 10^5 \lesssim Re_L \lesssim 10^7\)).

Lift–Drag Polar and Stall

\[ C_D = C_{D,0} + \frac{C_L^2}{\pi e\, AR}, \]

where \(e \lesssim 1\) is the Oswald efficiency factor. Maximum \(L/D\) occurs where profile and induced drag are equal.

Reynolds Analogy

\[ \frac{C_f}{2} \approx St_H \equiv \frac{h}{\rho c_p U_\infty}, \]

relating the skin-friction coefficient \(C_f\) to the Stanton number \(St_H\) — a remarkable free lunch for engineers estimating convective heat transfer from friction data.

Boundary-Layer Thickness Estimates

\[ \frac{\delta}{x} \approx \frac{5}{\sqrt{Re_x}}, \qquad Re_x = \frac{U_\infty x}{\nu}, \]\[ \frac{\delta^*}{x} \approx \frac{1.721}{\sqrt{Re_x}}, \qquad \frac{\theta}{x} \approx \frac{0.664}{\sqrt{Re_x}}. \]

Transition from laminar to turbulent typically occurs near \(Re_{x,\mathrm{tr}} \sim 5\times 10^5\) on smooth flat plates; freestream turbulence, pressure gradient, surface roughness, or a trip wire can drive transition forward. In turbulent flow, boundary-layer thickness grows approximately as \(\delta/x \approx 0.37\,Re_x^{-1/5}\), faster than laminar.

Chapter 20: Fluidization and Multiphase Flow

Packed-Bed Flow: Ergun Equation

Consider a fluid of density \(\rho\) and viscosity \(\mu\) flowing at superficial velocity \(u_s = Q/A_{\mathrm{bed}}\) through a stationary packed bed of particles of diameter \(d_p\) with bed voidage \(\varepsilon\) (void volume fraction). The Ergun equation gives the pressure drop per unit bed length \(L\):

The first term (viscous, Blake–Kozeny) dominates at low particle Reynolds number \(Re_p = \rho u_s d_p/\mu(1-\varepsilon) \ll 1\); the second (inertial, Burke–Plummer) dominates at \(Re_p \gg 10^3\).

Minimum Fluidization Velocity

\[ \Delta p = (1 - \varepsilon_{mf})(\rho_s - \rho)gL. \]\[ u_{mf} \approx \frac{(\rho_s - \rho)g\,d_p^2}{150\mu}\cdot\frac{\varepsilon_{mf}^3}{1 - \varepsilon_{mf}}. \]

Fluidization Regimes

Beyond \(u_{mf}\), as \(u_s\) increases further, successive regimes appear:

1. Fixed bed (\(u_s < u_{mf}\)) — particles stationary; Ergun applies.
2. Homogeneous (bubbling) fluidization (\(u_{mf} < u_s < u_{mb}\)) — uniform expansion for Geldart A powders; bubbles form beyond \(u_{mb}\).
3. Slugging — in narrow columns with deep beds, bubbles grow to fill the column cross-section.
4. Turbulent fluidization — bubbles break down; violent mixing.
5. Fast fluidization — substantial particle entrainment; requires solids recycling (circulating fluidized bed).
6. Pneumatic transport (\(u_s > u_t\), the terminal velocity) — dilute-phase gas-solids flow.

Terminal Settling Velocity

\[ \tfrac{\pi}{6}d_p^3(\rho_s-\rho)g = C_D\,\tfrac{1}{2}\rho u_t^2\,\tfrac{\pi}{4}d_p^2, \]\[ u_t = \sqrt{\frac{4(\rho_s-\rho)g d_p}{3\rho C_D(Re_p)}}. \]

Two limiting regimes:

• Stokes regime (\(Re_p < 1\), \(C_D = 24/Re_p\)): \(u_t = \dfrac{(\rho_s-\rho)g d_p^2}{18\mu}.\)
• Newton regime (\(10^3 < Re_p < 2\times 10^5\), \(C_D \approx 0.44\)): \(u_t = 1.74\sqrt{\dfrac{(\rho_s-\rho)g d_p}{\rho}}.\)

For intermediate \(Re_p\), \(C_D(Re_p)\) is implicit in \(u_t\); iteration or a direct correlation (e.g. Haider–Levenspiel) is used.

Gas–Liquid Two-Phase Flow

\[ \!\left(\frac{dp}{dz}\right)_{\!TP} = \phi_L^2\!\left(\frac{dp}{dz}\right)_{\!L,\mathrm{only}}, \qquad \phi_L^2 = 1 + \frac{C}{X} + \frac{1}{X^2}, \]

where \(X^2 = (dp/dz)_L/(dp/dz)_G\) is the Martinelli parameter and \(C\) depends on the laminar/turbulent character of each phase. These correlations are the engineering bread-and-butter of oil-and-gas pipeline, boiler, and condenser design.