Numerical Methods for Solving Navier-Stokes Equations

TL;DR

You'll learn why the Navier-Stokes equations resist exact solutions for real engineering flows and how engineers discretize them instead. You'll walk through the finite volume method, the pressure-velocity coupling problem, and the SIMPLE algorithm step by step. By the end you'll understand what's actually happening inside software like HEC-RAS, OpenFOAM, or ANSYS Fluent when it solves a flood or channel flow problem.

1. The Mental Model

The Navier-Stokes equations are a continuous statement of "mass and momentum are conserved everywhere, at every instant." A computer can't handle "everywhere" and "every instant" — it needs a finite list of numbers. So numerical methods chop space into a grid and time into steps, then solve approximate versions of those conservation laws at each point, feeding results forward in time. Numerically solving Navier-Stokes means trading an infinite continuous problem for a very large but finite algebra problem, then solving that algebra problem fast enough to be useful.

2. The Core Material

Why You Can't Just "Solve" Navier-Stokes

Photo by Brett Jordan on Pexels

The incompressible Navier-Stokes equations for a Newtonian fluid are:

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

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

Two things make this brutal to solve analytically for real civil engineering problems (river bends, spillways, culverts, storm surge): the $(\mathbf{u}\cdot\nabla)\mathbf{u}$ term is nonlinear, and pressure $p$ has no equation of its own — it only appears as a gradient, existing purely to enforce $\nabla \cdot \mathbf{u} = 0$. There's no separate "pressure equation" handed to you; you have to construct one. That's the root of most numerical difficulty here — it's not just "solve four equations," it's "solve four coupled equations where one variable is a Lagrange multiplier in disguise."

Analytical solutions exist only for toy geometries (Poiseuille flow, Couette flow). A flood plain, a stormwater network, or turbulent flow past a bridge pier has no closed form. So you discretize.

Discretization: Turning Calculus into Algebra

Photo by Nothing Ahead on Pexels

Three ingredients go into every numerical scheme:

1. Spatial discretization — chop the domain into cells (finite volume, most common in engineering CFD), points (finite difference), or elements (finite element, common in structural-hydraulic coupling). Each cell gets one value of velocity and pressure representing an average over that cell.

2. Temporal discretization — advance in discrete time steps $\Delta t$. Explicit schemes compute the new state directly from the old one (fast per step, but stability-limited). Implicit schemes solve a system involving both old and new states (more expensive per step, but stable for larger $\Delta t$).

3. Linearization — the nonlinear convection term is approximated using values from the previous iteration or time step, turning the nonlinear PDE into a sequence of linear algebra problems you solve repeatedly until they converge.

The finite volume method (FVM) is what almost all civil/water engineering CFD codes use, because it's built directly on conservation. For a control volume $V$ with surface $S$, mass conservation becomes:

$$\oint_S \mathbf{u}\cdot\mathbf{n}\, dA = 0$$

This says: flow in equals flow out, exactly, for every single cell. No fudging. That's why FVM is preferred over finite difference for engineering flows — it guarantees conservation at the discrete level, which matters when you're tracking a flood volume or a pollutant mass balance.

The Pressure-Velocity Coupling Problem

Photo by Едуард Ковтонюк on Pexels

Here's the classic headache. You have a momentum equation that needs a pressure field to compute velocity, but no direct equation for pressure. If you guess a pressure field, solve momentum for velocity, that velocity usually won't satisfy $\nabla \cdot \mathbf{u} = 0$. You need an iterative fix.

The SIMPLE algorithm (Semi-Implicit Method for Pressure-Linked Equations), developed by Patankar and Spalding, is the workhorse:

flowchart TD
    A["Guess pressure field p*"] --> B["Solve momentum equations for u*, v*, w*"]
    B --> C["Check: does velocity satisfy continuity?"]
    C -->|"No, div(u) != 0"| D["Solve pressure-correction equation for p'"]
    D --> E["Correct pressure: p = p* + alpha_p * p'"]
    E --> F["Correct velocity using p' gradient"]
    F --> G["Update boundary conditions, turbulence terms"]
    G --> H{"Converged? Residuals < tolerance"}
    H -->|"No"| A
    H -->|"Yes"| I["Move to next time step"]

The key insight in the pressure-correction step: you derive a Poisson-like equation for the pressure correction $p'$ by substituting the momentum equation into the continuity equation. Solving that equation tells you exactly how much to nudge the pressure guess so that the resulting velocity field is divergence-free (or closer to it). You repeat this loop — guess, solve, correct, check — until the mass imbalance in every cell drops below a tolerance, typically expressed as a residual.

Under-relaxation factors ($\alpha_p$, $\alpha_u$, typically 0.3–0.7) are essential here — apply the full correction each iteration and the solution oscillates or diverges, because the linearization is only an approximation of a nonlinear problem.

Stability: The CFL Condition

Photo by Wilson Ardila on Pexels

For explicit time-stepping schemes, you can't pick $\Delta t$ freely. The Courant-Friedrichs-Lewy (CFL) condition governs it:

$$C = \frac{u \, \Delta t}{\Delta x} \leq 1$$

Physically: information (a fluid parcel, a pressure wave) shouldn't cross more than one grid cell per time step, or the scheme loses track of it and blows up numerically. In a river model with $\Delta x = 1$ m cells and flow velocity $u = 3$ m/s, your time step must satisfy $\Delta t \leq \Delta x / u \approx 0.33$ s. This is why fine grids near structures (weirs, piers, culvert entrances) force small time steps everywhere unless you use implicit schemes or local time-stepping — a huge practical consideration in flood modelling where you're simulating hours of real time.

Turbulence: The Extra Complication

Real civil engineering flows (rivers, sewers, coastal currents) are turbulent at Reynolds numbers civil engineers actually see (often $10^4$–$10^7$). Directly resolving all turbulent eddies (Direct Numerical Simulation, DNS) is computationally impossible at these scales — DNS cost scales roughly as $Re^{9/4}$. So practical codes solve Reynolds-Averaged Navier-Stokes (RANS) equations instead, which average out turbulent fluctuations and model their effect with a turbulence model (k-ε, k-ω) that adds two extra transport equations to solve alongside momentum and continuity. This is why "solving Navier-Stokes" for a river almost always really means "solving RANS plus a turbulence closure," not the raw equations.

3. Worked Example

Let's set up a 1D simplified version of the pressure-correction idea for open-channel flow through a narrowing section — the kind of local acceleration you see at a bridge contraction. This strips away 3D complexity so you can see the coupling logic clearly.

Setup: Steady, incompressible flow through a channel narrowing from width $b_1 = 4$ m to $b_2 = 2$ m over a short reach, depth $h = 1.5$ m (constant, rectangular channel), discharge $Q = 6$ m³/s.

Step 1 — Continuity gives you velocities directly (1D case, no iteration needed yet):

$$u_1 = \frac{Q}{b_1 h} = \frac{6}{4 \times 1.5} = 1.0 \text{ m/s}, \qquad u_2 = \frac{Q}{b_2 h} = \frac{6}{2 \times 1.5} = 2.0 \text{ m/s}$$

Step 2 — Guess a pressure (here, water surface elevation) field. In open channel flow, "pressure" is hydrostatic and tied to depth. Guess a uniform depth $h^* = 1.5$ m everywhere (i.e., assume no drawdown through the contraction).

Step 3 — Solve momentum (Bernoulli with losses) using guessed depth. Energy conservation between section 1 and 2, ignoring losses for a first pass: