Numerical Methods for Solving Navier-Stokes Equations

TL;DR

You'll master the main computational approaches that turn the unsolvable Navier-Stokes equations into solvable engineering problems. We'll cover finite difference, finite element, and finite volume methods with their trade-offs. You'll understand when to use each method for real water and civil engineering applications.

1. The Mental Model

The Navier-Stokes equations describe fluid flow perfectly but can't be solved analytically for real problems. Numerical methods chop up space and time into tiny pieces, then approximate the derivatives using algebra instead of calculus. It's like replacing a smooth curve with thousands of tiny straight line segments—close enough for engineering.

2. The Core Material

2.1 Finite Difference Method (FDM)

The finite difference method is your starting point because it's conceptually simplest. You replace derivatives with difference quotients on a rectangular grid.

For the 1D momentum equation ∂u/∂t + u∂u/∂x = -1/ρ ∂p/∂x + ν∂²u/∂x², you'd discretize like this:

• Time derivative: (u^(n+1) - u^n)/Δt
• Spatial derivative: (u_(i+1) - u_(i-1))/(2Δx) (central difference)
• Second derivative: (u_(i+1) - 2u_i + u_(i-1))/Δx²

The beauty is you can set up these equations at every grid point and solve the resulting system of algebraic equations. The downside? FDM struggles with complex geometries—it works best for rectangular channels, pipes, or simple domains.

Stability is crucial. The Courant-Friedrichs-Lewy (CFL) condition demands that Δt ≤ CΔx/|u|, where C ≈ 1. Violate this and your solution explodes into nonsense.

2.2 Finite Element Method (FEM)

FEM shines when you're dealing with irregular boundaries—think river bends, bridge piers, or complex spillway geometries. Instead of a rigid grid, you divide your domain into triangular or tetrahedral elements.

The key insight: you approximate the velocity field u(x,y,t) as a weighted sum of shape functions N_i(x,y):

u(x,y,t) ≈ Σ u_i(t) N_i(x,y)

These shape functions are usually polynomials that equal 1 at one node and 0 at all others. You substitute this into the weak form of Navier-Stokes (derived using integration by parts) and get a system of ordinary differential equations in time.

FEM's strength is geometric flexibility—you can handle any dam, spillway, or hydraulic structure shape. The weakness? It's computationally expensive and requires more sophisticated programming.

2.3 Finite Volume Method (FVM)

FVM is the workhorse of commercial CFD codes like ANSYS Fluent and OpenFOAM. It's built around conservation—mass, momentum, and energy are conserved exactly across control volume boundaries.

You integrate the Navier-Stokes equations over each control volume:

∫∫∫ (∂ρ/∂t + ∇·(ρv)) dV = 0

Using the divergence theorem, volume integrals become surface integrals. This means fluxes entering one cell exactly equal fluxes leaving adjacent cells—perfect conservation.

flowchart TD
    A["Continuous Navier-Stokes PDE"] --> B["Choose Numerical Method"]
    B --> C["Finite Difference Method"]
    B --> D["Finite Element Method"]
    B --> E["Finite Volume Method"]
    C --> F["Regular Grid<br/>Simple Geometry"]
    D --> G["Irregular Mesh<br/>Complex Boundaries"]
    E --> H["Conservation Properties<br/>Industrial CFD"]
    F --> I["Explicit/Implicit<br/>Time Stepping"]
    G --> I
    H --> I
    I --> J["Solve Linear System<br/>Each Time Step"]
    J --> K["Check Convergence<br/>Update Solution"]
    K --> L["Engineering Results"]

Time stepping schemes matter for all methods:
- Explicit: Use known values to calculate new ones. Fast per step but small time steps required.
- Implicit: New values appear on both sides of equations. Larger stable time steps but expensive linear system solves.
- Semi-implicit: Treat diffusion implicitly, convection explicitly. Good compromise for many flows.

2.4 Handling Nonlinearity and Pressure-Velocity Coupling

The Navier-Stokes equations are nonlinear because velocity appears in its own derivative (u∂u/∂x). Most methods use iteration within each time step:

1. Guess velocity field
2. Solve momentum equations with old velocities
3. Check if continuity equation (∇·u = 0) is satisfied
4. Adjust pressure and velocity until convergence

The SIMPLE algorithm (Semi-Implicit Method for Pressure-Linked Equations) is standard:
- Solve momentum with guessed pressure
- Solve pressure correction equation
- Update velocities and pressure
- Repeat until convergence

3. Worked Example

Let's solve flow over a backward-facing step using finite differences—a classic test case for separated flow that occurs downstream of gates, weirs, and sudden expansions.

Problem Setup:
- Channel height H = 1 m, step height h = 0.5 m
- Inlet velocity U = 1 m/s, Reynolds number Re = 100
- Grid: 200×50 points, Δx = 0.1 m, Δy = 0.02 m

import numpy as np
import matplotlib.pyplot as plt

# Grid parameters
nx, ny = 200, 50
dx, dy = 0.1, 0.02
dt = 0.001
Re = 100.0
nu = 1.0/Re

# Initialize arrays
u = np.zeros((ny, nx))  # x-velocity
v = np.zeros((ny, nx))  # y-velocity
p = np.zeros((ny, nx))  # pressure
u_new = np.zeros((ny, nx))
v_new = np.zeros((ny, nx))

# Boundary conditions
u[25:, 0] = 1.0  # Inlet velocity above step
u[:, -1] = u[:, -2]  # Outlet: zero gradient
# Walls: u = v = 0 (already initialized)

# Main time loop (simplified)
for n in range(1000):
    # Solve momentum equations (explicit scheme)
    for i in range(1, ny-1):
        for j in range(1, nx-1):
            # Skip points inside step geometry
            if i < 25 and j < 10:
                continue

            # x-momentum equation
            dudx = (u[i, j+1] - u[i, j-1])/(2*dx)
            dudy = (u[i+1, j] - u[i-1, j])/(2*dy)
            d2udx2 = (u[i, j+1] - 2*u[i, j] + u[i, j-1])/dx**2
            d2udy2 = (u[i+1, j] - 2*u[i, j] + u[i-1, j])/dy**2
            dpdx = (p[i, j+1] - p[i, j-1])/(2*dx)

            u_new[i, j] = u[i, j] + dt*(-u[i, j]*dudx - v[i, j]*dudy 
                                       - dpdx + nu*(d2udx2 + d2udy2))

    # Apply boundary conditions
    u_new[25:, 0] = 1.0
    u_new[:, -1] = u_new[:, -2]

    # Update solution
    u = u_new.copy()

    # Check convergence (simplified)
    if n % 100 == 0:
        print(f"Time step {n}, max velocity: {np.max(u):.4f}")

# The recirculation length (key engineering parameter)
# Find where u = 0 along centerline behind step
centerline = u[15, 10:]  # Approximate centerline
recirculation_length = 0
for j, vel in enumerate(centerline):
    if vel > 0:
        recirculation_length = (j + 10) * dx
        break

print(f"Recirculation length: {recirculation_length:.2f} m")
print(f"Normalized by step height: {recirculation_length/0.5:.1f}")

Engineering Significance: The recirculation length behind the step is typically 6-8 times the step height for Re = 100. This matters for sediment deposition, mixing, and energy losses in channels with sudden expansions.

4. Key Takeaways

4.1 Most Important Concepts

• Grid resolution determines accuracy: You need at least 10-20 grid points across boundary layers and recirculation zones for meaningful results.
• Time step stability: Explicit schemes require Δt < CΔx/|u| while implicit schemes allow larger steps but cost more per step.
• Conservation is non-negotiable: Mass must be conserved exactly—if your continuity residual grows, your solution is wrong.
• Boundary conditions drive the physics: Inlet profiles, wall functions, and outlet conditions determine whether your solution represents reality.
• Iterative convergence: Nonlinear equations require inner iterations within each time step until residuals drop below tolerance.
• Validation against experiments: No numerical solution is trustworthy without comparison to physical data or analytical solutions.
• Computational cost scales rapidly: 3D problems with fine grids can require millions of grid points and hours of compute time.

4.2 Common Misconceptions

• "Finer grids always give better answers" → Grid-independent solutions exist; beyond that point, you're just burning CPU time without gaining accuracy.
• "Implicit schemes are always more stable" → They're more stable for diffusion-dominated flows, but can introduce numerical diffusion that smears sharp gradients.
• "Commercial codes always work correctly" → Default settings often fail for complex flows; you must understand the numerics to choose appropriate schemes and parameters.
• "Steady-state solutions converge quickly" → Separated flows and turbulent flows can take thousands of iterations to reach steady state, if they ever do.

4.3 Compare & Contrast

5. Now Try It

Set up and solve the driven cavity problem: a square box with the top wall moving at constant velocity while other walls are stationary. Use a 20×20 grid with finite differences, Re = 100, top wall velocity = 1 m/s. Implement the momentum equations with explicit time stepping and track how the center vortex develops. Plot velocity vectors every 50 time steps until steady state. Success looks like: a primary recirculation filling most of the cavity with smaller corner vortices, and your centerline velocities matching published benchmark data within 5%.