simpleFoam Residuals Not Converging: Complete Fix Guide

simpleFoam uses the SIMPLE algorithm (Semi-Implicit Method for Pressure-Linked Equations) to iterate toward a steady-state solution. When residuals stall, oscillate, or diverge, it almost always comes down to the same handful of settings — all of them in fvSolution and fvSchemes.

By CFDpilot · Updated April 2026

Understanding SIMPLE convergence

Before diving into fixes, it is important to understand what convergence means for simpleFoam. The SIMPLE algorithm solves a coupled pressure-velocity system by splitting it into a sequence of pressure and velocity corrections. Each iteration updates U, then p, then turbulence scalars. Convergence is achieved when these corrections become negligibly small — which is tracked by the residual.

The residual reported in the log is the normalised initial residual — the norm of the source term after assembly, normalised by the flux through the domain. A residual of 1e-4 for U means the remaining imbalance in the momentum equation is four orders of magnitude smaller than the total momentum flux. For engineering accuracy, residuals of 1e-4 to 1e-5 are typically sufficient.

Three distinct failure patterns require different fixes:

• Residuals diverge (climb toward 1 and beyond) — numerical instability, fix schemes and relaxation
• Residuals stall (plateau above 1e-4) — mesh quality or physics issue
• Residuals oscillate (periodic variation) — inherently unsteady flow, consider switching to pimpleFoam

1. Relaxation factors: the most common culprit

SIMPLE requires under-relaxation to converge. Without it, each iteration overshoots the solution and the residuals diverge. A missing or incorrect relaxationFactors block is responsible for the majority of simpleFoam failures.

// fvSolution — SIMPLE block
SIMPLE
{
    nNonOrthogonalCorrectors 0;
    consistent               yes;
}

relaxationFactors
{
    fields
    {
        p       0.3;
    }
    equations
    {
        U       0.7;
        k       0.5;
        epsilon 0.5;
        omega   0.5;
    }
}

If residuals oscillate without converging: reduce U to 0.5 and p to 0.2. If they still diverge: check your scheme first (see below).

SIMPLEC: faster convergence with less effort

Activating consistent yes switches the pressure correction algorithm from SIMPLE to SIMPLEC. The key difference is that SIMPLEC uses a more accurate approximation of the off-diagonal velocity coefficients in the pressure equation. This means the pressure and velocity corrections are better matched, and the solver can use higher relaxation factors without losing stability:

SIMPLE
{
    nNonOrthogonalCorrectors 1;
    consistent               yes;
}

relaxationFactors
{
    equations
    {
        U       0.9;
        k       0.7;
        epsilon 0.7;
    }
}

With SIMPLEC you do not need a separate pressure relaxation factor in the fields block — the algorithm handles it internally. Convergence with SIMPLEC is typically 30–50% faster in iteration count compared to SIMPLE with equivalent settings.

2. Unbounded div scheme

div(phi,U) Gauss linear is the most common scheme mistake. Gauss linear is central differencing — second order but unbounded. In combination with SIMPLE, it will diverge on any mesh that isn't perfectly uniform.

Replace with a bounded scheme:

divSchemes
{
    default         none;
    div(phi,U)      Gauss linearUpwind grad(U);
    div(phi,k)      Gauss upwind;
    div(phi,epsilon) Gauss upwind;
    div(phi,omega)  Gauss upwind;
    div((nuEff*dev2(T(grad(U))))) Gauss linear;
}

Use linearUpwind for accuracy, limitedLinear or upwind for stability on coarse or skewed meshes.

Scheme strategy for different mesh types

The right scheme depends heavily on mesh quality and flow regime:

• Structured hexahedral mesh, low non-orthogonality — Gauss linearUpwind grad(U) is optimal
• Unstructured tetrahedral mesh, moderate skewness — Gauss limitedLinearV 1 for U, Gauss upwind for scalars
• High non-orthogonality (>60°) — Gauss upwind for all convective terms until convergence, then switch to higher-order
• High-Re flow with strong gradients — consider Gauss LUST grad(U) for a blended 75% linear / 25% linearUpwind scheme

Never use default Gauss linear as a catch-all for divSchemes in a RANS simulation. Each term has different stability requirements.

3. Turbulence inlet values: k and epsilon

Incorrect inlet values for turbulent kinetic energy k and dissipation epsilon (or omega) are a frequent cause of convergence problems that look like solver instability but are actually a physics issue.

Estimating k and epsilon at the inlet

For an inlet with mean velocity U and turbulence intensity I (typically 1–5%):

// k = 1.5 * (U * I)^2
// Example: U=10 m/s, I=5% → k = 1.5 * (10*0.05)^2 = 0.375 m²/s²

// epsilon = C_mu^0.75 * k^1.5 / L
// L = turbulent length scale ≈ 0.07 * hydraulic_diameter
// Example: D_h=0.1m → L=0.007m, k=0.375 → epsilon = 0.09^0.75 * 0.375^1.5 / 0.007 = 14.2 m²/s³

// omega = epsilon / (C_mu * k)
// Example: omega = 14.2 / (0.09 * 0.375) = 420 s^-1

Never set k = 0 at the inlet. Zero k causes division by zero in most two-equation turbulence models. A very small value like 1e-10 is equally wrong — use physically realistic values.

Inlet boundary condition types for turbulence

For k and epsilon at a velocity inlet, use turbulentIntensityKineticEnergyInlet and turbulentMixingLengthDissipationRateInlet — these compute k and epsilon automatically from intensity and mixing length, avoiding manual errors:

// 0/k
inlet
{
    type            turbulentIntensityKineticEnergyInlet;
    intensity       0.05;        // 5% turbulence intensity
    value           uniform 0.375;
}

// 0/epsilon
inlet
{
    type            turbulentMixingLengthDissipationRateInlet;
    mixingLength    0.007;       // 0.07 * D_h
    value           uniform 14.2;
}

4. Pressure boundary conditions

For internal flows (pipe, channel, duct), a common mistake is setting fixedValue pressure at both inlet and outlet, or leaving the pressure reference undefined. simpleFoam requires exactly one pressure reference point.

• Inlet: zeroGradient for pressure, fixedValue for velocity
• Outlet: fixedValue uniform 0 for pressure, zeroGradient for velocity
• Or: use totalPressure at inlet + fixedValue at outlet for pressure-driven flows

Outlet boundary conditions for recirculating flows

Standard zeroGradient at the outlet assumes all flow exits the domain. When the outlet is too close to a recirculation zone, backflow enters through the outlet with undefined values and destabilises the simulation. Use inletOutlet instead:

// 0/U — outlet with backflow handling
outlet
{
    type            inletOutlet;
    inletValue      uniform (0 0 0);
    value           uniform (0 0 0);
}

// 0/k — outlet with backflow handling
outlet
{
    type            inletOutlet;
    inletValue      uniform 0.375;
    value           uniform 0.375;
}

The inletValue specifies what value to use when fluid enters through the outlet patch (backflow). This prevents backflow-induced divergence without requiring the outlet to be moved further downstream.

5. Convergence criteria too tight

If residuals reach a plateau around 1e-4 to 1e-5 and stop decreasing, the simulation may actually be converged — the remaining residual is due to numerical noise on the mesh, not a physical imbalance. Check your tolerance settings:

solvers
{
    p
    {
        solver          GAMG;
        tolerance       1e-6;
        relTol          0.01;
    }
    U
    {
        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-6;
        relTol          0.1;
    }
}

GAMG settings for pressure solver

GAMG (Generalised Algebraic MultiGrid) is the recommended solver for pressure in simpleFoam. It is significantly faster than PCG/ICCG for the pressure equation on large meshes. Key settings to tune:

solvers
{
    p
    {
        solver          GAMG;
        smoother        GaussSeidel;
        tolerance       1e-6;
        relTol          0.1;
        nCellsInCoarsestLevel 10;
    }
    "(U|k|epsilon|omega)"
    {
        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-6;
        relTol          0.1;
        nSweeps         1;
    }
}

The relTol (relative tolerance) means the inner solver stops when the residual has dropped by a factor of relTol from its initial value at that outer iteration. Setting relTol 0.1 is typically more efficient than relTol 0 because it avoids over-solving the inner equations at each outer iteration.

6. Monitoring convergence correctly

Do not rely solely on residuals to judge convergence. Residuals measure numerical imbalance, not physical accuracy. The correct approach is to monitor integrated quantities that are physically meaningful for your case:

// Add to controlDict — monitor forces and flow rates
functions
{
    forces
    {
        type            forces;
        libs            (forces);
        patches         (wall);
        rho             rhoInf;
        rhoInf          1.225;
        CofR            (0 0 0);
        writeControl    timeStep;
        writeInterval   10;
    }
    inlet_flowrate
    {
        type            surfaceFieldValue;
        libs            (fieldFunctionObjects);
        regionType      patch;
        name            inlet;
        operation       sum;
        fields          (phi);
        writeControl    timeStep;
        writeInterval   10;
    }
}

When forces and flow rates have stabilised to within 0.1% between successive checks, the simulation is converged regardless of what the residual plots show.