Numerical Methods for Solving Navier-Stokes Equations

TL;DR

You'll learn why the Navier-Stokes equations almost never have exact solutions for real civil engineering flows, and how engineers discretize space and time to solve them numerically instead. You'll walk through the finite volume method, the pressure-velocity coupling problem, and how a solver actually marches a flood or river simulation forward in time. You'll finish knowing how to read a CFD setup and judge whether its numerics are trustworthy.

1. The Mental Model

You can't solve Navier-Stokes with algebra for a river bending around a bridge pier — there's no formula for that. So instead you chop the river into thousands of tiny boxes, write approximate versions of the equations for each box, and solve all those approximate equations together, over and over, as time ticks forward. Every "CFD simulation" you've heard of is just this process done fast on a computer. Numerical methods turn one impossible continuous problem into millions of easy discrete ones.

2. The Core Material

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

The incompressible Navier-Stokes equations are:

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

The troublemaker is $(\mathbf{u}\cdot\nabla)\mathbf{u}$ — the convective term. It's nonlinear: velocity multiplies its own derivative. Nonlinear PDEs generally don't have closed-form solutions except in toy geometries (pipe flow, flow between plates). A river with a bridge pier, an irregular channel bed, or turbulent flow past a spillway gate has none of that symmetry. So you approximate.

There's a second headache specific to incompressible flow: pressure has no equation of its own. It's not a variable you evolve in time — it exists purely to enforce $\nabla \cdot \mathbf{u} = 0$ (mass conservation, no compression). This creates the "pressure-velocity coupling" problem discussed below.

Discretization: Turning Space into a Grid

The first job is to convert continuous space into a mesh of discrete points or cells. Three main approaches exist:

Finite Difference Method (FDM): Approximate derivatives directly using nearby grid point values. For a 1D velocity field:

$$\frac{\partial u}{\partial x} \approx \frac{u_{i+1} - u_{i-1}}{2\Delta x}$$

Simple, but struggles with irregular boundaries — bad for a natural river channel with a wandering bank.

Finite Volume Method (FVM): This is the workhorse in civil/water engineering (used by OpenFOAM, ANSYS Fluent, and most river/coastal solvers). Instead of approximating derivatives at points, you divide the domain into small control volumes and enforce conservation of mass and momentum across each volume's faces. This directly respects physics — what flows in must flow out or accumulate — which makes it forgiving of messy, irregular geometry like a floodplain.

Finite Element Method (FEM): Represents the solution as a sum of shape functions over elements. Common in structural-fluid coupling problems, less common for pure CFD in civil hydraulics, but you'll see it in groundwater flow and seepage problems.

For open-channel and river modeling, FVM dominates because mass conservation at each cell is exactly what you need to track water depth correctly.

Pressure-Velocity Coupling: SIMPLE and PISO

Here's the actual numerical trick that makes incompressible flow solvable. Since pressure doesn't have its own transport equation, you can't just march it forward like velocity. Instead, solvers use an iterative correction procedure.

The SIMPLE algorithm (Semi-Implicit Method for Pressure-Linked Equations) works like this:

1. Guess a pressure field $p^*$.
2. Solve the momentum equation using $p^*$ to get a velocity field $u^*$ — but this $u^*$ generally does NOT satisfy $\nabla \cdot \mathbf{u} = 0$.
3. Derive a pressure correction equation from the continuity error: essentially, wherever mass isn't conserved, that's telling you the pressure guess was wrong.
4. Correct pressure and velocity: $p = p^* + p'$, $u = u^* + u'$.
5. Repeat until the continuity residual (mass imbalance) drops below a tolerance, say $10^{-5}$.

PISO (Pressure-Implicit with Splitting of Operators) does the same idea but with extra correction loops per time step, making it more stable for transient problems — which is why it's preferred for unsteady flows like a dam-break wave or storm surge, where SIMPLE alone would need very small time steps to stay accurate.

flowchart TD
    A["Guess pressure field p*"] --> B["Solve momentum equation for u*"]
    B --> C["Check continuity: is div(u*) = 0?"]
    C -->|No, mass imbalance found| D["Solve pressure correction equation p'"]
    D --> E["Update p = p* + p', u = u* + u'"]
    E --> C
    C -->|Yes, residual below tolerance| F["Advance to next time step"]

Time Stepping and Stability: The CFL Condition

Once space is discretized, you still need to march forward in time. The simplest approach, explicit Euler, computes tomorrow's velocity directly from today's:

$$u^{n+1} = u^n + \Delta t \cdot (\text{RHS terms at } t^n)$$

This is cheap per step but has a strict stability limit called the Courant-Friedrichs-Lewy (CFL) condition:

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

where $C_{max}$ is typically around 1 for explicit schemes. Physically: information (a fast-moving parcel of water) shouldn't cross more than one grid cell in a single time step, or your numerical scheme loses track of it and the solution blows up. This is the single most common cause of a diverging CFD simulation — someone picked $\Delta t$ too large relative to the local flow speed and mesh size.

Implicit schemes (like backward Euler or Crank-Nicolson) solve for $u^{n+1}$ using equations that reference $u^{n+1}$ itself, requiring a matrix solve at each step. They're more expensive per step but unconditionally stable (or nearly so), so you can take much larger time steps — valuable for slow-evolving processes like groundwater seepage or reservoir sedimentation over months.

Turbulence: Why You Rarely Solve "Real" Navier-Stokes

For turbulent flows — nearly every civil engineering flow of interest, from river rapids to stormwater pipes — fully resolving every eddy (Direct Numerical Simulation, DNS) is computationally impossible at real scale. The number of grid points needed scales as $Re^{9/4}$, where $Re$ is the Reynolds number. A river might have $Re$ in the millions; DNS at that scale needs more grid points than exist in any current supercomputer.

So engineers use turbulence models that solve a modified, time-averaged version of Navier-Stokes:

• RANS (Reynolds-Averaged Navier-Stokes): Averages out turbulent fluctuations, adds a turbulence model (k-ε, k-ω) to approximate their effect. Cheap, standard for most practical river and pipe flow design.
• LES (Large Eddy Simulation): Resolves large eddies directly, models only the small ones. More accurate, much more expensive — used for research-grade sediment transport or scour studies.

For a bridge scour study or spillway design, RANS with a k-ε model is the typical industry default because it balances accuracy against the computing budget of a real engineering project.

3. Worked Example

Let's set up a simplified 1D scenario: flow in a straight, uniform channel segment with a small pressure-driven acceleration, ignoring viscosity for a moment to focus purely on the discretization and stability logic (this is the kind of sanity check you'd do before trusting a full 2D/3D CFD run).

Setup: Channel velocity $u = 2 \text{ m/s}$, cell size $\Delta x = 0.5\text{ m}$. We want to know the maximum time step allowed by CFL with $C_{max} = 1$.

$$\Delta t \leq \frac{C_{max}\, \Delta x}{u} = \frac{1 \times 0.5}{2} = 0.25\text{ s}$$

So if your solver output shows $\Delta t = 0.4\text{ s}$ for this mesh and flow speed, you should immediately expect instability — the Courant number is $C = 2 \times 0.4/0.5 = 1.6$, above 1, and the simulation is likely to diverge or produce spurious oscillations.

Now let's code the actual explicit finite-difference update for 1D advection of velocity (a simplified stand-in for the convective term in Navier-Stokes) to see the CFL limit bite in practice:

import numpy as np
import matplotlib.pyplot as plt

# Domain setup
L = 10.0          # channel length (m)
nx = 40           # number of grid points
dx = L / (nx - 1)
u_speed = 2.0      # advection speed (m/s), like channel velocity

# Two time steps to compare: one CFL-safe, one CFL-violating
def run_advection(dt, nt):
    x = np.linspace(0, L, nx)
    # initial condition: a velocity "bump" (like a surge entering the channel)
    u = np.exp(-((x - 2.0)**2) / 0.5)
    C = u_speed * dt / dx
    print(f"dt={dt}, Courant number C={C:.2f}")

    for n in range(nt):
        u_new = u.copy()
        for i in range(1, nx - 1):
            # upwind explicit scheme for du/dt + u_speed*du/dx = 0
            u_new[i] = u[i] - C * (u[i] - u[i - 1])
        u = u_new
    return x, u

x1, u_stable = run_advection(dt=0.1, nt=20)   # C = 0.8, stable
x2, u_unstable = run_advection(dt=0.3, nt=20) # C = 2.4, unstable

plt.plot(x1, u_stable, label="dt=0.1 (C=0.8, stable)")
plt.plot(x2, u_unstable, label="dt=0.3 (C=2.4, unstable)")
plt.legend(); plt.xlabel("x (m)"); plt.ylabel("u")
plt.show()

Running this, dt=0.1 gives a smooth bump that shifts downstream cleanly — physically sensible. dt=0.3 prints C=2.4 and the resulting u_unstable array will show wild oscillating values, some far outside the initial [0,1] range — the numerical signature of CFL violation. This tiny example is exactly the same failure mode that crashes full 3D river CFD models when someone sets too coarse a mesh with too large a time step.

4. Key Takeaways

4.1 Most Important Concepts

• Nonlinearity forces numerics: the convective term $(\mathbf{u}\cdot\nabla)\mathbf{u}$ has no general