Learning #26 – Fluid Dynamics & Turbulence: Navier-Stokes, Reynolds Decomposition, Boundary Layers and the Kolmogorov Cascade

The Navier-Stokes equations were written down in the 1820s–1840s and describe all of classical fluid mechanics. A million-dollar Clay Mathematics Institute Millennium Prize awaits whoever proves (or disproves) that smooth solutions always exist for three-dimensional incompressible flow. Yet turbulence — which requires such solutions — has already been understood well enough to design aircraft, submarines and gas turbines. This productive tension between mathematical incompleteness and engineering success is what makes fluid dynamics so fascinating.

1. The Navier-Stokes Equations

The Navier-Stokes equations express Newton’s second law for a fluid parcel: the rate of change of momentum equals the sum of pressure gradients, viscous stresses, and body forces. For an incompressible Newtonian fluid (constant density ρ, dynamic viscosity μ), they take the compact vector form:

Continuity (incompressibility):
  ∇ · u = 0    (ρ = const)

Momentum:
  ρ(∂u/∂t + (u·∇)u) = −∇p + μ∇²u + ρg
  ┌┐               ┌┐        ┌┐      ┌┐         ┌┐
  inertia       advection   pressure viscous   body force

Dimensionless form (scale by U, L):
  ∂u*/∂t* + (u*·∇*)u* = −∇*p* + (1/Re)∇*²u*

Reynolds number Re = ρUL/μ = UL/ν    (ν = μ/ρ = kinematic viscosity)
  Re << 1: viscous forces dominate → laminar, creeping flow (Stokes flow)
  Re ~ 100: laminar with steady recirculation (cylinder wake)
  Re ~ 10^3: periodic vortex shedding (Kármán street, St = fD/U ≈ 0.2)
  Re > 10^4: transitional;
  Re > 10^5: fully turbulent

Bernoulli equation (inviscid, along streamline):
  p + ½ρu² + ρgz = const   (energy per volume)
  Origin: ∇(½u²) = (u·∇)u + u×(∇×u)

The non-linear advection term (u·∇)u is the source of almost all complexity in fluid mechanics. It transfers kinetic energy between scales in a way that the linear viscous term cannot undo, giving rise to eddies, wakes and ultimately turbulence. At low Reynolds numbers this term is negligible compared to viscosity and flows are smooth and predictable. At high Reynolds numbers viscosity can dissipate only the finest scales; the large-scale motions are essentially inviscid and their energy must pass down through progressively smaller eddies before it can be converted to heat.

2. Boundary Layer Theory

When a fluid flows over a solid surface, the no-slip condition requires that the fluid velocity matches the wall velocity at the surface. The rapid transition from zero (at the wall) to the free-stream velocity U occurs across a thin boundary layer, first analysed by Ludwig Prandtl in 1904. This insight unlocked practical aerodynamics: outside the thin boundary layer, flow is approximately inviscid; inside, viscous effects dominate.

Blasius laminar boundary layer (flat plate, zero pressure gradient):
  δ(x) / x = 5.0 / √Re_x       (99% thickness)
  δ*(x) / x = 1.72 / √Re_x      (displacement thickness: missing mass flux)
  C_f(x) = 0.664 / √Re_x            (local skin friction coefficient)

Transition to turbulence (flat plate):
  Re_x,crit ≈ 5×10^5 (freestream turbulence level ~0.1%)
  Higher freestream turbulence or adverse pressure gradient → earlier transition

Turbulent boundary layer (power-law profile):
  u/U ≈ (y/δ)^{1/7}    (1/7 power law, Re > 10^5)
  C_f = 0.0592 Re_x^{-1/5}

Law of the wall (inner layer):
  u^+ = y^+                    (viscous sublayer, y^+ < 5)
  u^+ = (1/κ) ln(y^+) + B   (log-law region, 30 < y^+ < 200)
  y^+ = y u_τ / ν,   u_τ = √(τ_w/ρ) (friction velocity)
  κ ≈ 0.41 (von Kármán constant),  B ≈ 5.0

Boundary layer separation:
  Occurs when pressure gradient dP/dx > 0 (adverse pressure gradient)
  Near-wall fluid decelerates to zero; reversal → separation bubble → wake
  Separation control: suction, vortex generators, turbulators (golf ball dimples)

3. Hydrodynamic Instabilities

A laminar flow is unstable if small perturbations grow over time rather than decay. Linear stability analysis determines whether a given base flow is stable by inserting a small perturbation and asking whether the linearised equations admit growing solutions. Two instabilities are especially important for turbulence transition and mixing layer dynamics: Kelvin-Helmholtz and Rayleigh-Taylor.

Kelvin-Helmholtz instability (shear across interface):
  Two inviscid fluids with velocities U_1, U_2 and densities ρ_1, ρ_2 meeting at y=0
  Growth rate for wavenumber k:
    σ^2 = −gk(ρ_1−ρ_2)/(ρ_1+ρ_2) + k²ρ_1ρ_2(U_1−U_2)²/(ρ_1+ρ_2)²
  Pure shear (same density): always unstable for any k > 0
  Stabilised by gravity & density stratification (Richardson number Ri = N²/S²)
  Ri < 1/4 (Miles-Howard criterion): shear overcomes stratification → KH rolls

Rayleigh-Taylor instability (heavy fluid over light):
  Heavy fluid (ρ_2) resting on lighter fluid (ρ_1) in gravitational field g
  Growth rate: σ = √[Atkg]    (A = Atwood number = (ρ_2−ρ_1)/(ρ_2+ρ_1))
  All wavelengths unstable, but short waves grow faster; surface tension stabilises k > k_c
  Nonlinear: mushroom-shaped plumes of heavy fluid penetrating downward (spikes) and light fluid rising (bubbles)
  Applications: supernova ejecta, inertial confinement fusion capsule implosion, ocean mixing

Rayleigh-Bénard convection (buoyancy-driven):
  Fluid layer heated from below, cooled above
  Rayleigh number: Ra = gβΔT L³ / (να)  (β = thermal expansion, α = diffusivity)
  Onset (conduction → convection rolls): Ra_c = 1708  (Chandrasekhar 1961)
  Ra > 10^4: oscillatory; Ra > 10^6: turbulent plumes
  Nusselt number Nu = h L / λ ~ Ra^{1/3} (turbulent regime)

4. Reynolds Decomposition and RANS

Fully turbulent flows contain eddies across an enormous range of scales. Directly simulating every length and timescale (Direct Numerical Simulation, DNS) requires grid points proportional to Re^{9/4} — for a jet engine turbine at Re ~ 10^7, that is ~10^{15 } points, far beyond present computing power. Engineering models instead decompose the velocity into a mean and fluctuating part: the Reynolds decomposition.

Reynolds decomposition:
  u_i(x,t) = U_i(x) + u_i'(x,t)
  <u_i> = U_i  (time average),  <u_i'> = 0

Reynolds-averaged momentum equation:
  ρ U_j ∂U_i/∂x_j = −∂P/∂x_i + ∂/∂x_j [μ∂U_i/∂x_j − ρ<u_i'u_j'>]
                                                └──────────────────────└
                                                Reynolds stress tensor τ^R_ij

Closure problem:
  τ^R_ij has 6 independent unknowns per point → more unknowns than equations
  Boussinesq hypothesis: τ^R_ij = −2ν_t S_ij  (ν_t = turbulent eddy viscosity)
  Analogy with molecular viscosity; fails in strongly anisotropic flows

k-ε model (Launder & Spalding, 1974 — most widely used RANS model):
  Transport equation for turbulent kinetic energy: k = ½<u_i'u_i'>
  Transport equation for dissipation rate: ε = ν<∂u_i'/∂x_j ∂u_i'/∂x_j>
  Eddy viscosity: ν_t = C_μ k² / ε  (C_μ = 0.09)
  Free-stream turbulence length scale: L_t = k^{3/2}/ε

k-ω SST model (Menter 1994 — standard for aerodynamics):
  Blends k-ω near wall (accurate in sublayer) with k-ε in freestream
  Better prediction of adverse-pressure-gradient separation, airfoil stall
  Used in ANSYS Fluent, OpenFOAM, SU2

5. Kolmogorov’s Theory of the Energy Cascade

In 1941, Andrei Kolmogorov formulated a statistical theory of fully developed turbulence that remains the foundation of the field despite its simplicity. The central idea is the Richardson energy cascade: energy is injected at large scales (integral scale L), transferred to progressively smaller eddies without significant dissipation, and finally dissipated as heat at the viscous Kolmogorov scale η. In the intermediate inertial subrange, the statistics of velocity fluctuations are universal.

Kolmogorov microscales (only depend on ν and ε):
  Length: η = (ν³/ε)^{1/4}
  Time:   τ_η = (ν/ε)^{1/2}
  Velocity: u_η = (νε)^{1/4}

Scale separation:
  L / η ~ Re^{3/4}       (integral scale to Kolmogorov scale)
  For atmospheric boundary layer (Re ~ 10^8): L/η ~ 10^6

First similarity hypothesis (local isotropy for small scales):
  At scales r << L, statistics of velocity increments are universal
  (determined only by ν and ε)

Second similarity hypothesis (inertial subrange, η << r << L):
  Statistics depend only on ε (viscosity irrelevant)
  Structure function: S_2(r) = <|u(x+r) − u(x)|²> = C_2 (εr)^{2/3}
  C_2 ≈ 2.0 (Kolmogorov constant)

Energy spectrum (Kolmogorov −5/3 law):
  E(k) = C_K ε^{2/3} k^{−5/3}   (inertial subrange)
  C_K ≈ 1.5 (Kolmogorov spectral constant)
  First verified experimentally by Grant et al. (1962) in tidal channel

Intermittency corrections (K62, Obukhov 1962):
  Real turbulence is spatially intermittent; dissipation is fractal
  Scaling exponent: ζ_p = p/3 − μp(p−3)/18   (μ ≈ 0.25)
  Multifractal formalism (Parisi & Frisch 1985)

DNS capabilities (2024):
  Maximum DNS Re_λ ~ 2000 (Taylor-scale Reynolds number)
  Grid: 12 288³ ≈ 2×10^12 points  (runs on top-10 supercomputers)
  Goal: Re_λ ≈ 10^4 (comparable to atmospheric surface layer) within ~15 yr

6. Large-Eddy Simulation and Applications

Between the brute-force accuracy of DNS and the speed of RANS sits Large-Eddy Simulation (LES): directly compute eddies larger than the grid, and model only the subgrid-scale effects. LES captures unsteady mixing, coherent structures and aeroacoustics that RANS misses entirely, at a fraction of the cost of DNS. Modern industrial CFD increasingly uses hybrid RANS-LES methods such as Detached Eddy Simulation (DES) and Scale-Adaptive Simulation (SAS).

Grid filter (bar denotes filtered quantity):
  ū(x) = ∫ G(x − x') u(x') dx'    (G = box or Gaussian filter, width Δ)
  ū contains scales > Δ;  subgrid stress τ_ij^sgs = ​u_i u_j – ū_i ū_j

Smagorinsky model (1963):
  τ_ij^sgs = −2(C_s Δ)² |S̄| S̄_ij
  C_s ≈ 0.1–0.18 (universal for isotropic turbulence; must be reduced near walls)

Dynamic Smagorinsky (Germano 1991):
  C_s computed locally from ratio of test-filter to grid-filter stresses
  Removes need for manual tuning; handles transition automatically

Wall-modelled LES (WMLES):
  Model near-wall region with RANS or algebraic law-of-wall → coarser grid near surface
  N_grid ~ Re^{13/7} (WMLES) vs Re^{37/14} (wall-resolved LES) vs Re^{9/4} (DNS)
  WMLES of full aircraft at Re ~ 10^8 feasible with ~10^9 points

Applications:
  Combustion:     LES resolves turbulent flame brushing, premixed swirl burners
  Wind energy:    LES of atmospheric boundary layer + turbine wakes (SOWFA, AMR-Wind)
  Aeroacoustics:  pressure fluctuations from LES + FW-H equation → far-field noise
  Climate:        ocean mesoscale eddies (~30 km) now resolved in CMIP7 hi-res models

Related Simulations