When Pressure Rings in a Checkerboard — SIMPLE and Rhie-Chow Coupling

The velocity field looked fine. The pressure field did not. Every cell jumped against its neighbor. A checkerboard. Residuals kept falling. The pressure still rang like a sawtooth. Not a bug. It was baked in the moment pressure and velocity shared a grid point.

Toggle Naive against Rhie-Chow and the seeded sawtooth either sits still or sinks onto a smooth curve.

On Naive the amplitude does not shrink. On Rhie-Chow the same checkerboard settles. Reseed checkerboard plants the sawtooth again.

Why the tooth survives, and how continuity becomes a pressure equation, comes next. The compressible stretch of the same solver is later.

Continuity cannot see the pressure next door#

On a collocated grid (velocity and pressure both stored at cell centers), the semi-discrete momentum equation reads

where $a_P$ is the diagonal coefficient, $\mathbf{H}(\mathbf{u})$ collects the neighbor contributions, $V_P$ is the cell volume, and $(\nabla p)_P$ is the cell-center pressure gradient.

The trouble is that cell-center gradient. Written with central differences, $(\partial p/\partial x)_P = (p_E - p_W)/2\Delta x$. Notice $p_P$ is absent. A cell never sees its own pressure — only the neighbors two cells away. So odd and even cells split apart (odd-even decoupling). The $(+,-,+,-,\dots)$ sawtooth sits in the null space of the discrete equations and produces no residual at all.

Rhie-Chow: build the face velocity from momentum#

The cure is to stop averaging cell velocities to the face. Reapply the momentum equation at the face itself, so the neighboring pressure difference across that face enters directly.

The overbar is linear interpolation. The second term is the pressure gradient recomputed at the face. Now the face velocity is tied to $p_E - p_P$, the difference between adjacent cells. The sawtooth mode has nowhere left to hide.

As long as the face jump $p_E-p_P$ is not zero, the tooth is caught by that term. Another way to remove pressure from the momentum equation is Chorin's projection method.

SIMPLE: turn continuity into a pressure equation#

The root difficulty of incompressible flow is that there is no equation that solves for pressure directly. SIMPLE (Semi-Implicit Method for Pressure Linked Equations) turns continuity $\nabla\cdot\mathbf{u}=0$ into one. Substituting the Rhie-Chow face velocity into continuity yields an elliptic (Poisson) equation for pressure.

The left side is a Laplacian (a diffusion form); the right side is the divergence of the predicted velocity. The procedure is predictor-corrector.

1. Solve momentum with a guessed $p^*$ to predict $\mathbf{u}^*$.
2. Solve the pressure equation above to update $p$.
3. Correct face and cell velocities with the new pressure gradient so continuity holds.
4. Repeat 1-3 until convergence.

PISO adds two or three more corrector passes per step and is used for transient runs. PIMPLE nests an outer SIMPLE loop and an inner PISO loop for large time steps.

Without under-relaxation it blows up#

There is a catch. If you replace the pressure fully every step ($p \leftarrow p^*$), the coupling over-corrects and diverges. SIMPLE only lets a fraction of the update through.

$\alpha_p$ is the pressure under-relaxation factor (0 to 1). Velocity gets its own $\alpha_u$. Too small and convergence crawls; too large and it oscillates and bursts. A common rule of thumb is $\alpha_p \approx 0.3$ and $\alpha_u \approx 0.7$, chosen so the two roughly sum to one.

Change the relaxation factor yourself in the lab below. It solves the pressure-correction equation with a relaxed Gauss-Seidel sweep.

Drop the factor to 0.3 and the residual curve crawls down gently. Near 1.5 it drops fastest. Push it toward 2 and over-correction piles up each sweep, so the residual jumps back up. SIMPLE's $\alpha_p$ lives on exactly this same balance.

One solver for every speed — the compressible reach#

The real appeal of the pressure-based method is that it does not care about the Mach number. In compressible flow continuity becomes a density equation, and an equation of state (EOS) ties density to pressure. Writing $\psi \equiv \partial\rho/\partial p$ (the compressibility, $=1/c^2$) adds a time term to the pressure equation.

This equation carries both convective and diffusive character. As $M\to 0$, $\psi$ grows so the time term does not dominate and pressure is solved elliptically (instantaneous, global coupling). At large $M$ it turns hyperbolic (acoustic waves travel at finite speed). Density-based solvers stiffen at low Mach because the density-pressure coupling weakens, but the pressure-based method keeps that coupling explicit through the EOS. That is how one code spans subsonic to supersonic.

Reviving and erasing the checkerboard in Python#

Both stencils run on a 1-D periodic grid. The question is whether the checkerboard survives.

import numpy as np
 
def poisson_sweep(p, f, naive):
    """One Gauss-Seidel sweep of -p'' = f on a 1D periodic grid."""
    N = len(p)
    for i in range(N):
        if naive:                       # naive linear interp: skip-a-cell, decoupled stencil
            p[i] = 0.5 * (p[(i - 2) % N] + p[(i + 2) % N] + f[i])
        else:                           # Rhie-Chow: compact 3-point stencil couples neighbors
            p[i] = 0.5 * (p[(i - 1) % N] + p[(i + 1) % N] + f[i])
    p -= p.mean()                       # pressure is free up to a constant -> pin the mean
    return p
 
def checkerboard_metric(p):
    """Size of the (+,-,+,-,...) component. Zero means no sawtooth."""
    signs = (-1.0) ** np.arange(len(p))
    return abs(np.dot(p, signs)) / len(p)
 
N = 48
x = 2 * np.pi * np.arange(N) / N
f = np.sin(x) + 0.4 * np.sin(2 * x)
f -= f.mean()
 
for naive in (True, False):
    p = 0.8 * (-1.0) ** np.arange(N)    # seed a checkerboard
    for _ in range(4000):
        p = poisson_sweep(p, f, naive)
    tag = "naive linear" if naive else "Rhie-Chow  "
    print(f"{tag}  checkerboard = {checkerboard_metric(p):.2e}")
 
# naive linear   checkerboard = 8.00e-01   <- the sawtooth stays
# Rhie-Chow      checkerboard = 3.1e-16    <- gone to machine precision

Same source, same initial guess, same iteration count. Only the stencil changed. Naive cannot remove the sawtooth. Rhie-Chow erases it. Same place as the figure above.

The initial pressure is $0.8\times(-1)^i$. The metric is the absolute inner product of $p$ with that sign pattern, divided by the cell count, so it starts at 0.8. After 4000 sweeps on 48 cells, naive is still $8.00\times 10^{-1}$. Not a single digit moved. Rhie-Chow is $3.1\times 10^{-16}$. Machine precision.

The stencil is why naive keeps 0.8. On a sawtooth, index $i$ and $i\pm 2$ share a sign. $p_i=0.5(p_{i-2}+p_{i+2}+f_i)$ does not average that mode away. The cell-center gradient $(p_E-p_W)/2\Delta x$ does the same thing. $E$ and $W$ are both opposite $P$, so $p_E-p_W=0$. A tooth of amplitude 0.8 has zero gradient at the center. It never shows up in the residual. The face jump $p_E-p_P$ flips sign, so its size is 1.6. That is the term Rhie-Chow keeps. The smooth part of $f=\sin x+0.4\sin 2x$ gets solved by both stencils. The only leftover is the tooth.

Set alpha_p to 1 and the corrector blows first#

Replacing pressure in full, $p\leftarrow p^*$, over-corrects the coupling in one step. The update only lets a fraction $\alpha_p$ through. The usual pair $\alpha_p\approx 0.3$, $\alpha_u\approx 0.7$ sums to 1. On the lab curve, 0.3 crawls, near 1.5 is the fastest drop, and near 2 the residual jumps back up every sweep. SIMPLE keeps $\alpha_p$ between 0 and 1, on the slow side of that curve.

Where changing the relaxation changes the converged field is the face velocity that remembers under-relaxation. On a collocated grid, averaging cell velocities onto the face drops the neighbor pressure jump, and the tooth lives outside the residual. Reviving the time term with $\psi=\partial\rho/\partial p$ comes after that. Add compressibility to a pressure equation that still carries the tooth, and only the sound speed gets faster. The checkerboard stays.