The SIMPLE Algorithm on a Staggered Grid

Most steady-state incompressible flow solvers — including OpenFOAM's simpleFoam — use the SIMPLE algorithm proposed by Patankar and Spalding (1972). Understanding SIMPLE is essential for both writing your own solvers and interpreting convergence behaviour in commercial codes.

The Core Problem

For steady incompressible flow, we need to solve:

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

The difficulty is that $p$ and $\mathbf{u}$ are coupled: the momentum equation needs $p$ to solve for $\mathbf{u}$, but there is no explicit equation for $p$. SIMPLE resolves this through an iterative predictor–corrector loop.

The MAC Staggered Grid

Before describing SIMPLE, the grid arrangement is crucial. On a uniform 2-D grid of cells with width $\Delta x$ and height $\Delta y$:

• $u_{i+1/2,j}$ is stored at the east face of cell $(i,j)$
• $v_{i,j+1/2}$ is stored at the north face of cell $(i,j)$
• $p_{i,j}$ is stored at the cell centre

This eliminates pressure–velocity decoupling. The discrete continuity equation naturally involves face velocities:

$$\frac{u_{i+1/2,j} - u_{i-1/2,j}}{\Delta x} + \frac{v_{i,j+1/2} - v_{i,j-1/2}}{\Delta y} = 0$$

The SIMPLE Loop

Step 1: Solve Momentum with Guessed Pressure

Given an initial (or previous-iteration) pressure field $p^*$, solve the momentum equations for a preliminary velocity field $\mathbf{u}^*$:

$$A_u\,u^*_{i+1/2,j} = \sum_{\text{nb}} a_{\text{nb}}\,u^*_{\text{nb}} - \frac{p^*_{i+1,j} - p^*_{i,j}}{\Delta x}$$

where $A_u$ is the central coefficient from the finite-volume discretisation of convection and diffusion, and $a_{\text{nb}}$ are the neighbour coefficients.

The velocity field $\mathbf{u}^*$ will generally not satisfy continuity.

Step 2: Derive the Pressure Correction Equation

Define corrections: $p = p^* + p'$ and $\mathbf{u} = \mathbf{u}^* + \mathbf{u}'$.

The SIMPLE approximation drops the neighbour velocity correction terms (this is the "semi-implicit" approximation — the source of the algorithm's name and its main approximation):

$$u'_{i+1/2,j} \approx -\frac{\Delta y}{A_u}\,(p'_{i+1,j} - p'_{i,j})$$

Substituting into the discrete continuity equation gives the pressure correction equation:

$$a_P\,p'_{i,j} = a_E\,p'_{i+1,j} + a_W\,p'_{i-1,j} + a_N\,p'_{i,j+1} + a_S\,p'_{i,j-1} + b_{i,j}$$

where $b_{i,j}$ is the mass imbalance (residual of continuity) at cell $(i,j)$.

Step 3: Correct Pressure and Velocity

$$p^{n+1} = p^* + \alpha_p\,p'$$ $$u^{n+1}_{i+1/2,j} = u^*_{i+1/2,j} - \frac{\Delta y}{A_u}\,(p'_{i+1,j} - p'_{i,j})$$

The relaxation factor $\alpha_p < 1$ (typically 0.3) stabilises the iteration. Under-relaxation is essential — without it, SIMPLE diverges.

Step 4: Check Convergence and Repeat

Compute the residuals (sum of absolute mass and momentum imbalances). If they are below a tolerance, stop. Otherwise, set $p^* = p^{n+1}$, $\mathbf{u}^* = \mathbf{u}^{n+1}$, and repeat from Step 1.

Connection to OpenFOAM

In OpenFOAM's simpleFoam, the SIMPLE loop appears in UEqn.H and pEqn.H. The pressure correction equation is solved with GAMG (algebraic multigrid), and velocity is solved with PBiCGStab. The number of SIMPLE iterations is controlled by nNonOrthogonalCorrectors in fvSolution.

// pEqn.H (simplified)
volScalarField rAU(1.0/UEqn.A());
surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));

U = rAU * UEqn.H();
phi = fvc::flux(U) + rAUf * fvc::ddtCorr(U, phi);

fvScalarMatrix pEqn(fvm::laplacian(rAUf, p) == fvc::div(phi));
pEqn.solve();
phi -= pEqn.flux();
U -= rAU * fvc::grad(p);

Why SIMPLE Converges (and When It Does Not)

SIMPLE converges when the pressure correction equation is diagonally dominant — which is guaranteed by the staggered grid arrangement and the dropping of neighbour corrections in the velocity correction step. However:

• High Reynolds number flows require smaller relaxation factors, slowing convergence.
• Highly skewed meshes break the assumption that face normals align with cell-centre lines, introducing non-orthogonal correction terms.
• Pressure–velocity coupling is not unique — the reference pressure must be fixed at one cell (e.g., the outlet) to make the system non-singular.

What SIMPLEC and PISO Add

SIMPLEC (Consistent) retains the neighbour correction terms in a simplified form, allowing higher relaxation factors and faster convergence.

PISO (Pressure Implicit with Splitting of Operators) performs multiple pressure corrections per time step, making it suitable for transient flows. OpenFOAM's pisoFoam uses this.

Understanding the differences matters when setting up OpenFOAM cases: for steady problems, simpleFoam; for transient, pisoFoam or pimpleFoam (PISO + SIMPLE hybrid).