Harish Ragul — CFD & Computational Engineering

The incompressible Navier–Stokes equations present a fundamental challenge: the pressure appears in the momentum equation but has no independent evolution equation of its own. Chorin's projection method (1968) resolves this by splitting the time advance into two fractional steps.

The Equations

Find the velocity field $\mathbf{u}(\mathbf{x}, t)$ and pressure $p(\mathbf{x}, t)$ satisfying:

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

\nabla \cdot \mathbf{u} = 0

The divergence-free constraint (continuity) is the difficulty: $p$ acts as a Lagrange multiplier to enforce it, but solving for $p$ and $\mathbf{u}$ simultaneously requires a coupled system.

Step 1: Intermediate Velocity (Ignore Pressure)

Advance the velocity field from time $n$ to an intermediate state $\mathbf{u}^*$ using only advection and diffusion — ignoring the pressure gradient entirely:

\frac{\mathbf{u}^* - \mathbf{u}^n}{\Delta t} = -(\mathbf{u}^n \cdot \nabla)\mathbf{u}^n + \nu\,\nabla^2\mathbf{u}^n

This is a simple explicit update. The result $\mathbf{u}^*$ will generally not satisfy $\nabla \cdot \mathbf{u}^* = 0$ .

Step 2: Pressure Poisson Equation

We want to find a pressure correction $\phi$ such that the corrected velocity:

\mathbf{u}^{n+1} = \mathbf{u}^* - \Delta t\,\nabla \phi

satisfies $\nabla \cdot \mathbf{u}^{n+1} = 0$ . Taking the divergence:

\nabla \cdot \mathbf{u}^{n+1} = \nabla \cdot \mathbf{u}^* - \Delta t\,\nabla^2 \phi = 0

This gives the pressure Poisson equation:

\nabla^2 \phi = \frac{1}{\Delta t}\,\nabla \cdot \mathbf{u}^*

The right-hand side is computable from $\mathbf{u}^*$ . Solving this elliptic equation yields $\phi$ , which acts as (or is proportional to) the pressure.

Step 3: Project onto Divergence-Free Space

Correct the velocity:

\mathbf{u}^{n+1} = \mathbf{u}^* - \Delta t\,\nabla\phi

Update the pressure: $p^{n+1} = \phi$ (in the simplest formulation).

This projection step is why the method is called the projection method — it projects $\mathbf{u}^*$ onto the space of divergence-free vector fields.

Python Implementation (Lid-Driven Cavity)

import numpy as np

def build_rhs(u, v, dx, dy, dt, rho):
    """RHS of pressure Poisson equation."""
    b = np.zeros_like(u)
    b[1:-1, 1:-1] = (rho / dt) * (
        (u[1:-1, 2:] - u[1:-1, :-2]) / (2 * dx) +
        (v[2:, 1:-1] - v[:-2, 1:-1]) / (2 * dy)
    )
    return b

def pressure_poisson(p, b, dx, dy, nit=50):
    """Solve pressure Poisson with Gauss-Seidel."""
    for _ in range(nit):
        pn = p.copy()
        p[1:-1, 1:-1] = (
            (pn[1:-1, 2:] + pn[1:-1, :-2]) * dy**2 +
            (pn[2:, 1:-1] + pn[:-2, 1:-1]) * dx**2 -
            b[1:-1, 1:-1] * dx**2 * dy**2
        ) / (2 * (dx**2 + dy**2))
        # Neumann BC: dp/dn = 0 on walls
        p[:, -1] = p[:, -2]
        p[0, :] = p[1, :]
        p[:, 0] = p[:, 1]
        p[-1, :] = 0   # fixed pressure reference
    return p

Why the Staggered Grid Matters

In the implementation above, storing $u$ , $v$ , and $p$ all at the same grid nodes leads to the checkerboard pressure instability — a spurious oscillation where pressure alternates sign on neighbouring cells without affecting the velocity.

The fix is a MAC (Marker-and-Cell) staggered grid:

$u$ stored at left/right cell faces (staggered in $x$ )
$v$ stored at top/bottom cell faces (staggered in $y$ )
$p$ stored at cell centres

This arrangement means the discrete divergence of velocity and the discrete gradient of pressure are naturally paired — eliminating the instability. The SIMPLE algorithm, which I cover in a companion note, uses this grid by default.

Convergence and Stability

The explicit scheme above is conditionally stable. The CFL condition requires:

\Delta t \leq \frac{1}{2}\min\!\left(\frac{\Delta x^2}{\nu},\,\frac{\Delta y^2}{\nu}\right)

For flows at higher Reynolds number (smaller $\nu$ ), this forces very small time steps. The SIMPLE algorithm avoids this by treating the pressure implicitly.

What I Learned

Building this solver from scratch — rather than using a black-box tool — exposed something important: the pressure Poisson solve is the computational bottleneck. The elliptic solver (Gauss–Seidel above) converges slowly; in production codes, multigrid methods or conjugate gradient with preconditioning are used instead. This is one of the reasons OpenFOAM uses algebraic multigrid (GAMG) for pressure.

Projection Method for Incompressible Navier–Stokes: Step by Step