OpenFOAM Simulation Diverging: Causes and Fixes

Divergence is the most common failure mode in OpenFOAM. Your residuals spike, the solver crashes with a floating point exception, and the log shows values going to infinity. There are four systematic root causes — and each has a clear fix.

By CFDpilot · Updated April 2026

1. Unbounded convection scheme in fvSchemes

The single most frequent cause of divergence in steady-state simulations is using Gauss linear for the convective term div(phi,U). Gauss linear is unbounded — it can produce values outside the physical range of the solution, which compounds with each iteration until the field blows up.

Fix: Use a bounded scheme. For most incompressible RANS cases:

divSchemes
{
    div(phi,U)      Gauss linearUpwind grad(U);
    div(phi,k)      Gauss upwind;
    div(phi,epsilon) Gauss upwind;
}

linearUpwind is second-order accurate and bounded. For high-Re flows with coarse meshes, pure upwind is more stable (first-order but guaranteed bounded).

Understanding scheme boundedness

A numerical scheme is bounded when the computed cell value is guaranteed to fall within the range of its neighbouring values. Bounded schemes prevent the local extrema that trigger divergence. The trade-off is that unbounded schemes can be more accurate on very regular meshes, but this accuracy is irrelevant if the simulation diverges.

The limitedLinear scheme offers a middle ground: it blends between central differencing (accurate, unbounded) and upwind (diffusive, bounded) based on a limiter function. It is a good choice when you want more accuracy than pure upwind without the instability of central differencing:

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

The 1 argument to limitedLinearV controls the limiting aggressiveness: 1 is maximum limiting (closer to upwind), 0 is minimum limiting (closer to central). Start with 1 for stability and reduce only after the simulation runs cleanly.

2. Missing or wrong relaxation factors in fvSolution

In steady-state solvers (simpleFoam, incompressibleFluid), under-relaxation is the primary stability mechanism. If relaxationFactors are missing or set to 1.0, the SIMPLE algorithm has no damping and diverges immediately.

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

Standard starting values: U = 0.7, p = 0.3. If still diverging, reduce both to 0.5 and 0.2 respectively.

How under-relaxation works in SIMPLE

The SIMPLE algorithm solves for a velocity and pressure correction at each iteration. Without damping, the correction applied at each step can overshoot the solution and push the field further from equilibrium. Under-relaxation scales the correction by a factor between 0 and 1: a factor of 0.7 for U means only 70% of the computed correction is applied.

The price of aggressive under-relaxation (low factors) is slower convergence — more iterations are needed. The benefit is stability. For a new case, always start with the conservative values above and increase them only once the simulation is running stably.

Note that the consistent keyword in the SIMPLE block activates the SIMPLEC algorithm variant, which allows higher relaxation factors for pressure. With consistent yes, you can often use p = 0.9 instead of 0.3:

SIMPLE
{
    nNonOrthogonalCorrectors 1;
    consistent               yes;
}

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

3. Courant number too high (transient simulations)

In transient solvers (pimpleFoam, icoFoam, interFoam), the Courant number Co = U·Δt/Δx must stay below 1 for numerical stability. A Co of 5 or 10 means the solver is trying to advance information faster than the mesh can resolve it, leading to divergence.

Fix: Use adjustable time stepping in controlDict:

adjustTimeStep  yes;
maxCo           0.8;
maxDeltaT       0.01;

Or manually reduce deltaT until the Courant max stays below 1. Check the log for lines like Courant Number mean: X max: Y.

Estimating a safe deltaT before running

You can estimate a safe initial time step before the simulation starts. From your mesh, identify the minimum cell size (available in checkMesh output under minimum face area, or estimate from your refinement level). From your boundary conditions, identify the maximum expected velocity. Then:

// Safe deltaT estimate
// deltaT = Co_target * dx_min / U_max
// Example: U_max=10 m/s, dx_min=0.5mm=0.0005m, Co_target=0.8
// deltaT = 0.8 * 0.0005 / 10 = 4e-5 s

Use this as your starting deltaT with adjustable time stepping enabled. The solver will then refine it automatically.

4. Mesh quality issues

High non-orthogonality, severe skewness, or negative cell volumes can cause divergence regardless of how well the numerics are configured. Run checkMesh and look for:

• Max non-orthogonality > 70° — add nNonOrthogonalCorrectors 2 or 3 in fvSolution
• Max skewness > 4 — the mesh needs to be rebuilt
• Negative cell volumes — fatal, the mesh is invalid

For non-orthogonal meshes, also use corrected or limited Laplacian schemes instead of uncorrected.

Interpreting the checkMesh report

The key section to examine in the checkMesh output is the mesh quality metrics block. A clean mesh produces this output:

Mesh OK.

Mesh non-orthogonality Max: 38.5  average: 12.1
*   Max skewness = 1.2 — OK.
    Max aspect ratio = 8.3 — OK.

A problematic mesh looks like this:

***High aspect ratio cells found, Max aspect ratio: 2340
***Number of severely non-orthogonal (> 70 degrees) faces: 127.
***Max skewness = 5.8 — Too many highly skewed faces.

Every line starting with *** is a warning. Even a single cell with a severe defect can cause divergence because that cell will blow up first and propagate NaN values through the field. Fix these warnings before adjusting any numerical settings.

Non-orthogonality correctors

The nNonOrthogonalCorrectors setting tells the pressure solver to perform additional correction sweeps to account for the non-orthogonal component of the Laplacian operator. Each corrector adds computational cost but improves the accuracy of the pressure equation on non-orthogonal meshes:

SIMPLE
{
    nNonOrthogonalCorrectors 2;
}

laplacianSchemes
{
    default  Gauss linear corrected;
}

Use 0 correctors for max non-orthogonality below 40°, 1–2 for 40–70°, and 2–3 for 70–85°. Above 85° the mesh must be repaired — correctors alone are not enough.

5. Boundary condition errors

Incorrectly specified boundary conditions can drive divergence by injecting unphysical values into the domain. Common mistakes include:

• Setting fixedValue pressure at both inlet and outlet simultaneously (creates an over-constrained pressure problem)
• Using pressure in Pascal (Pa) when the incompressible solver expects kinematic pressure in m²/s²
• Setting k=0 or epsilon=0 at turbulent inlets (causes immediate division by zero in the turbulence model)
• Listing a patch in 0/U that does not exist in the constant/polyMesh/boundary file

Always verify that every patch in your boundary files corresponds to an actual mesh boundary. Mismatched patch names cause OpenFOAM to silently use a default type, which is often wrong for the physics.

Systematic diagnosis checklist

1. Check fvSchemes: is div(phi,U) bounded?
2. Check fvSolution: are relaxationFactors present and reasonable?
3. Check the log: what is the max Courant number?
4. Run checkMesh: any warnings on non-orthogonality or skewness?
5. Check 0/: are boundary conditions physically consistent?
6. Check 0/k and 0/epsilon: are turbulence fields non-zero and physically realistic?

Reading the divergence pattern in the log

The pattern of how residuals behave before the crash gives strong diagnostic clues:

• Diverges immediately at iteration 1 — mesh issue (negative volumes) or zero turbulence initial conditions
• Diverges after 5–20 iterations — relaxation factors too high or unbounded div scheme
• Diverges after 100+ iterations — subtle boundary condition issue or turbulence model instability
• Residuals oscillate but never diverge — the flow is not steady-state; switch to transient solver
• Only pressure residual diverges — pressure boundary condition problem, missing reference cell

When the pressure residual diverges while velocity stays bounded, it almost always means there is no unique pressure solution — either both inlet and outlet have fixed pressure, or the domain is fully enclosed with no pressure reference. Add pRefCell 0; pRefValue 0; to the SIMPLE block to anchor the pressure:

SIMPLE
{
    nNonOrthogonalCorrectors 1;
    pRefCell                 0;
    pRefValue                0;
}