3.4 Navier-Stokes equations

Derivation of Navier-Stokes equations

The Navier-Stokes equations are the fundamental governing equations for viscous fluid motion. They combine conservation of mass and momentum with constitutive relations to form a coupled system of partial differential equations relating velocity, pressure, density, and other fluid properties. Nearly every problem in fluid dynamics starts here.

Assumptions and simplifications

Before deriving the equations, several assumptions narrow the scope:

• Continuum assumption: The fluid is treated as a continuous medium, not as discrete molecules. This holds whenever the mean free path of molecules is much smaller than the characteristic length scale of the flow (i.e., the Knudsen number $Kn \ll 1$).
• Newtonian fluid assumption: Shear stress is linearly proportional to strain rate through a constant dynamic viscosity. This covers most common fluids (water, air, oils) but excludes non-Newtonian fluids like blood or polymer solutions.
• Incompressibility assumption: Fluid density is constant throughout the flow. Valid for liquids and for gases at low Mach numbers ($Ma < 0.3$).
• Isothermal flow: Temperature is uniform and constant, so the energy equation decouples from the momentum equations.

Each assumption you relax adds complexity. The full compressible, non-isothermal form requires solving an energy equation alongside mass and momentum conservation.

Conservation of mass

Mass cannot be created or destroyed within a control volume. The general continuity equation is:

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

For incompressible flow ($\rho = \text{const}$), the density drops out and you get the divergence-free condition:

$\nabla \cdot \mathbf{u} = 0$

This constraint says that the net volumetric flux out of any infinitesimal fluid element is zero. It's deceptively simple but plays a critical role: it couples with the momentum equation to determine the pressure field.

Conservation of momentum

Apply Newton's second law to an infinitesimal fluid element. The rate of change of momentum equals the sum of all forces:

$\rho \left( \frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla)\mathbf{u} \right) = -\nabla p + \nabla \cdot \boldsymbol{\tau} + \mathbf{f}$

Breaking down each term:

1. $\rho \frac{\partial \mathbf{u}}{\partial t}$ is the local (unsteady) acceleration.
2. $\rho (\mathbf{u} \cdot \nabla)\mathbf{u}$ is the convective (advective) acceleration, the nonlinear term that makes these equations so difficult to solve.
3. $-\nabla p$ is the pressure gradient force.
4. $\nabla \cdot \boldsymbol{\tau}$ represents viscous forces (diffusion of momentum).
5. $\mathbf{f}$ captures body forces like gravity ($\rho \mathbf{g}$).

Constitutive equations

The momentum equation contains the stress tensor $\boldsymbol{\tau}$, which needs to be related to the velocity field. This is where the constitutive equation closes the system.

For a Newtonian fluid, the relationship is linear:

$\boldsymbol{\tau} = 2\mu \mathbf{D}$

where $\boldsymbol{\tau}$ is the viscous stress tensor, $\mu$ is the dynamic viscosity, and $\mathbf{D} = \frac{1}{2}(\nabla \mathbf{u} + (\nabla \mathbf{u})^T)$ is the symmetric strain rate tensor.

Substituting this into the momentum equation for an incompressible Newtonian fluid gives the familiar form:

$\rho \left( \frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla)\mathbf{u} \right) = -\nabla p + \mu \nabla^2 \mathbf{u} + \mathbf{f}$

The viscous term simplifies to $\mu \nabla^2 \mathbf{u}$ because $\nabla \cdot \mathbf{u} = 0$ for incompressible flow.

For non-Newtonian fluids, more complex models replace the linear relation. Examples include power-law fluids (where viscosity depends on strain rate) and Bingham plastics (which require a yield stress before flowing).

Components of Navier-Stokes equations

Each field variable in the Navier-Stokes equations has a distinct physical role. Understanding what each one represents helps you interpret solutions and set up problems correctly.

Velocity field

The velocity field $\mathbf{u}(\mathbf{x}, t) = (u, v, w)$ in 3D Cartesian coordinates gives the speed and direction of fluid motion at every point in space and time. It's the primary unknown you're solving for. The velocity field determines how fluid elements are advected (carried along) and deformed.

Pressure field

Pressure $p(\mathbf{x}, t)$ is a scalar field representing force per unit area acting normal to any surface within the fluid. The pressure gradient $\nabla p$ drives fluid from high-pressure to low-pressure regions. In incompressible flow, pressure doesn't have its own evolution equation; instead, it acts as a Lagrange multiplier that enforces the divergence-free constraint $\nabla \cdot \mathbf{u} = 0$.

Viscous stress tensor

The viscous stress tensor $\boldsymbol{\tau}$ is a symmetric second-order tensor that captures both shear and normal stresses arising from viscosity. For Newtonian fluids, it's directly proportional to velocity gradients through the constitutive relation. Physically, viscous stresses represent the diffusion of momentum between adjacent fluid layers moving at different speeds.

Body forces

Body forces $\mathbf{f}(\mathbf{x}, t)$ are external forces acting on every fluid element throughout the volume. The most common example is gravity ($\mathbf{f} = \rho \mathbf{g}$), but electromagnetic forces (in magnetohydrodynamics) and Coriolis forces (in rotating reference frames) also fall into this category. They appear as source terms in the momentum equation.

Incompressible vs compressible flow

The distinction between incompressible and compressible flow determines which form of the Navier-Stokes equations you use. The deciding factor is whether density variations matter for the physics you're trying to capture.

Incompressible Navier-Stokes equations

When density is constant, the equations simplify considerably:

• Continuity: $\nabla \cdot \mathbf{u} = 0$
• Momentum: $\rho \left( \frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla)\mathbf{u} \right) = -\nabla p + \mu \nabla^2 \mathbf{u} + \mathbf{f}$

Pressure here isn't a thermodynamic variable; it's determined entirely by the requirement that velocity stays divergence-free. This is valid for most liquid flows and for gas flows at Mach numbers below about 0.3 (where density changes are less than ~5%).

Compressible Navier-Stokes equations

When density varies significantly, you need the full system:

• Continuity: $\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{u}) = 0$
• Momentum: $\frac{\partial (\rho \mathbf{u})}{\partial t} + \nabla \cdot (\rho \mathbf{u} \otimes \mathbf{u}) = -\nabla p + \nabla \cdot \boldsymbol{\tau} + \mathbf{f}$
• Energy: An additional equation governs temperature and internal energy changes.

You also need an equation of state (e.g., the ideal gas law $p = \rho R T$) to close the system, since pressure, density, and temperature are now all coupled.

Mach number considerations

The Mach number $Ma = U/c$ (ratio of flow velocity $U$ to the local speed of sound $c$) is the key parameter:

• $Ma < 0.3$: Incompressible treatment is accurate.
• $0.3 < Ma < 0.8$: Subsonic but compressible; density variations become non-negligible.
• $0.8 < Ma < 1.2$: Transonic regime with mixed subsonic/supersonic regions and possible shocks.
• $Ma > 1$: Supersonic flow with shock waves and expansion fans.
• $Ma > 5$: Hypersonic flow where real-gas effects (dissociation, ionization) may matter.

All transonic, supersonic, and hypersonic flows require the compressible Navier-Stokes equations.

Boundary conditions

Boundary conditions specify how the fluid behaves at the edges of your computational or analytical domain. Without proper boundary conditions, the Navier-Stokes equations don't have a unique solution.

No-slip condition

At a solid wall, the fluid velocity matches the wall velocity:

$\mathbf{u} = \mathbf{u}_\text{wall}$

For a stationary wall, this means $\mathbf{u} = 0$. This condition arises from the fact that viscous fluids "stick" to surfaces at the molecular level. It's the most common wall boundary condition and is responsible for the formation of boundary layers.

Free-slip condition

The normal velocity component matches the boundary's normal velocity, but the tangential component is unconstrained (no shear stress at the wall):

$\mathbf{u} \cdot \mathbf{n} = \mathbf{u}_\text{wall} \cdot \mathbf{n}, \quad \boldsymbol{\tau} \cdot \mathbf{n} \times \mathbf{n} = 0$

This applies to inviscid flow models, symmetry planes, and free surfaces where tangential stress is negligible.

Inflow and outflow conditions

• Inflow: You typically prescribe velocity (Dirichlet condition) or both velocity and pressure at the inlet. The specific choice depends on whether the flow is subsonic or supersonic.
• Outflow: Conditions here are trickier because you don't want artificial reflections contaminating the interior solution. Common approaches include zero-gradient (Neumann) conditions, advective outflow conditions, or characteristic-based boundary conditions that let waves pass through cleanly.

Symmetry and periodicity

• Symmetry conditions let you solve only half (or a quarter, etc.) of a geometrically symmetric problem, cutting computational cost. The normal velocity and normal gradients of tangential quantities are set to zero on the symmetry plane.
• Periodic conditions assume the flow repeats identically in one or more directions. The solution at one periodic boundary is mapped to the opposite boundary. These are widely used in channel flows, turbomachinery blade passages, and homogeneous turbulence simulations.

Dimensionless form

Non-dimensionalizing the Navier-Stokes equations reduces the number of independent parameters and reveals which physical effects dominate. You scale all variables by characteristic quantities: a reference length $L$, velocity $U$, time $L/U$, and pressure $\rho U^2$.

Reynolds number

$Re = \frac{\rho U L}{\mu}$

The Reynolds number is the ratio of inertial forces to viscous forces and is the single most important dimensionless number in fluid dynamics. It determines the flow regime:

• Low $Re$ (say, $Re < 1$): Viscous forces dominate. Flow is laminar and often reversible (Stokes flow).
• Moderate $Re$: Laminar flow with inertial effects.
• High $Re$ (above a critical threshold that depends on geometry): Flow transitions to turbulence. For pipe flow, the critical $Re \approx 2300$.

Strouhal number

$St = \frac{f L}{U}$

The Strouhal number relates the frequency $f$ of unsteady phenomena (like vortex shedding) to the convective time scale $L/U$. For a circular cylinder, $St \approx 0.2$ over a wide range of Reynolds numbers, meaning the shedding frequency is predictable from the flow speed and cylinder diameter.

Froude number

$Fr = \frac{U}{\sqrt{g L}}$

The Froude number compares inertial forces to gravitational forces. It governs free-surface flows (rivers, ship wakes, dam breaks) and stratified flows. When $Fr < 1$ (subcritical flow), gravity waves can propagate upstream; when $Fr > 1$ (supercritical), they cannot.

Dimensionless Navier-Stokes equations

After non-dimensionalization, the incompressible Navier-Stokes equations become:

$St \frac{\partial \mathbf{u}^*}{\partial t^*} + (\mathbf{u}^* \cdot \nabla^*)\mathbf{u}^* = -\nabla^* p^* + \frac{1}{Re} \nabla^{*2} \mathbf{u}^* + \frac{1}{Fr^2} \hat{\mathbf{g}}$

where starred quantities are dimensionless. The power of this form is that two flows with the same $Re$, $St$, and $Fr$ are dynamically similar regardless of their physical size, speed, or fluid. This is the theoretical basis for wind tunnel testing and scale-model experiments.

Numerical methods for solving

The nonlinearity of the convective term and the pressure-velocity coupling make analytical solutions rare. Most practical problems require numerical methods that discretize the equations on a computational grid.

Finite difference methods

Derivatives are approximated by finite differences on a structured grid (e.g., central differences, upwind schemes). These methods are straightforward to implement and computationally efficient for simple, regular geometries. The main drawbacks are numerical diffusion (which smears sharp gradients) and difficulty handling complex geometries without coordinate transformations.

Finite volume methods

The equations are integrated over discrete control volumes, and fluxes across cell faces are computed to ensure conservation of mass, momentum, and energy at the discrete level. Finite volume methods handle unstructured grids and complex geometries well, making them the backbone of most commercial CFD codes (ANSYS Fluent, OpenFOAM).

Finite element methods

The solution is approximated using basis functions within each element, and a weighted residual (weak) formulation is solved. Finite element methods offer high-order accuracy and geometric flexibility. However, they require careful stabilization for convection-dominated flows (e.g., SUPG or GLS stabilization) and tend to be more computationally expensive per degree of freedom than finite volume methods.

Spectral methods

The solution is expanded as a sum of global basis functions (Fourier modes for periodic domains, Chebyshev or Legendre polynomials otherwise). For smooth solutions, spectral methods converge exponentially fast, giving very high accuracy with relatively few degrees of freedom. The trade-off is that they're largely restricted to simple geometries and struggle with discontinuities (Gibbs phenomenon).

Turbulence modeling

Turbulent flows contain eddies spanning a huge range of spatial and temporal scales. Resolving all of them directly is often impractical, so various modeling strategies trade off between accuracy and computational cost.

Direct numerical simulation (DNS)

DNS solves the full Navier-Stokes equations on a grid fine enough to capture every scale of turbulence, from the largest energy-containing eddies down to the smallest dissipative (Kolmogorov) scales. No modeling is involved, so DNS provides the most accurate results. The cost scales roughly as $Re^3$ in 3D, which limits DNS to moderate Reynolds numbers (currently up to about $Re \sim 10^4$ for simple geometries on large supercomputers).

Large eddy simulation (LES)

LES explicitly resolves the large, energy-carrying eddies and models only the small-scale (subgrid-scale) motions using a subgrid-scale model (e.g., Smagorinsky model, dynamic model). A spatial filter is applied to the Navier-Stokes equations to separate resolved from unresolved scales. LES is much cheaper than DNS but still demands significant computational resources, especially near walls where turbulent structures become very small.

Reynolds-averaged Navier-Stokes (RANS) models

RANS decomposes every flow variable into a time-averaged mean and a fluctuation (Reynolds decomposition: $\mathbf{u} = \bar{\mathbf{u}} + \mathbf{u}'$). Substituting into the Navier-Stokes equations and averaging produces equations for the mean flow, but with extra unknown terms: the Reynolds stresses $\overline{u_i' u_j'}$. RANS is by far the cheapest approach and is the workhorse of industrial CFD, but its accuracy depends heavily on the turbulence model chosen.

Closure problem and turbulence models

The Reynolds stresses that appear in the RANS equations are unknowns, creating more unknowns than equations. This is the closure problem. Turbulence models provide the missing relationships:

• Algebraic models (e.g., mixing length model): Simple, zero-equation models that specify an eddy viscosity based on local flow properties. Limited to simple shear flows.
• One-equation models (e.g., Spalart-Allmaras): Solve one transport equation for a modified eddy viscosity. Popular in aerospace applications.
• Two-equation models (e.g., $k$-$\epsilon$, $k$-$\omega$, SST): Solve transport equations for turbulent kinetic energy $k$ and a second variable (dissipation rate $\epsilon$ or specific dissipation rate $\omega$). The SST model blends $k$-$\omega$ near walls with $k$-$\epsilon$ in the freestream and is widely used.
• Reynolds stress models (RSM): Solve transport equations for each component of the Reynolds stress tensor. More accurate for flows with strong anisotropy (swirling flows, separation) but more expensive and harder to converge.

Applications and examples

The Navier-Stokes equations underpin analysis across a wide range of engineering and scientific problems.

Pipe flow and pressure drop

Viscous flow through pipes is one of the few cases with an exact analytical solution (Hagen-Poiseuille flow for laminar conditions). The velocity profile is parabolic, and the pressure drop is proportional to the flow rate and inversely proportional to $r^4$ (where $r$ is the pipe radius). For turbulent pipe flow, empirical correlations like the Moody chart relate the friction factor to $Re$ and surface roughness. These results are essential for designing oil and gas pipelines, water distribution networks, and HVAC systems.

Flow over airfoils and lift generation

When fluid flows over an airfoil, the curved upper surface accelerates the flow, lowering pressure (by Bernoulli's principle in the inviscid region), while the lower surface sees higher pressure. This pressure difference produces lift. The Navier-Stokes equations capture the full picture, including viscous effects in the boundary layer, the Kutta condition at the trailing edge, and stall behavior at high angles of attack. Aircraft wing design, wind turbine blades, and propeller optimization all rely on solving these equations.

Boundary layer theory and separation

Near a solid surface, viscous effects are confined to a thin boundary layer where the velocity transitions from zero (at the wall, due to no-slip) to the freestream value. The boundary layer can be laminar or turbulent, and its thickness grows along the surface. When the pressure gradient becomes sufficiently adverse (pressure increasing in the flow direction), the boundary layer can separate from the surface, creating a recirculation zone that dramatically increases drag and can cause loss of lift on airfoils.

Vortex shedding and wake dynamics

When flow passes a bluff body (like a cylinder), the boundary layers on either side separate and roll up into alternating vortices that shed periodically into the wake. This is the von Kármán vortex street. The shedding frequency is characterized by the Strouhal number ($St \approx 0.2$ for a cylinder). The oscillating forces from vortex shedding can cause structural vibrations, which is a critical concern in the design of bridges, chimneys, offshore platforms, and heat exchanger tubes.