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 a handful of systematic root causes — convection schemes, relaxation, Courant number, mesh, boundary conditions, turbulence, the linear solver and thermophysical bounds — and each has a clear fix. This guide walks through all of them, plus how to read the divergence pattern in the log to find which one is yours.

By CFDpilot · Updated June 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.

div scheme comparison: stability vs accuracy

When a case is blowing up, move down this table toward "more bounded"; once it runs cleanly, move back up for accuracy:

Scheme	Order	Bounded?	Use it when
Gauss linear	2nd	No	Only on very regular meshes / DNS-LES. Diverges on most RANS.
Gauss linearUpwind grad(U)	2nd	Mostly	Default for div(phi,U) in steady RANS. Best accuracy/stability balance.
Gauss limitedLinear 1	1st–2nd	Yes	Scalars (k, ω, ε, T) and tough cases. Blends central/upwind.
Gauss upwind	1st	Yes	Rescue scheme. Fully stable but diffusive — use to get a first solution, then refine.

The classic recovery move: switch everything to upwind, get the case running and roughly converged, then step div(phi,U) back up to linearUpwind for the final accurate solution. See the fvSchemes guide for a full recommended set per solver type.

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.

6. Turbulence model instability

When momentum and pressure look stable but the case still blows up, the turbulence fields are often the cause. If k, omega or epsilon go negative, OpenFOAM clips them and prints warnings in the log:

bounding k, min: -3.4e+02 max: 9.1e+02 average: 2.7e+01
bounding omega, min: -1.2e+04 max: ...

A few bounding lines early on are normal. Persistent or growing bounding means the turbulence field is unstable and nut (= k/ω for k-ω models) will eventually explode and take the whole solution with it. The usual causes and fixes:

• Unbounded scalar scheme — div(phi,k) on a central scheme produces the negative values. Use Gauss upwind or Gauss limitedLinear 1 for all turbulence scalars.
• Unrealistic inlet values — k=0 or a tiny omega gives a near-zero or huge nut. Set k = 1.5·(U·I)² and a consistent omega = k^0.5/(Cμ^0.25·L).
• Wall-function / y+ mismatch — high-Re wall functions on a y+≈1 mesh (or vice-versa) corrupt the near-wall turbulence. Match them in the wall functions guide.
• Relaxation too high — drop k and omega relaxation to 0.3–0.5.

See turbulence model selection if you are unsure whether the model itself fits your flow regime.

7. Linear (pressure) solver settings

Sometimes the divergence is inside the linear solve itself: the pressure residual climbs instead of dropping within a single outer iteration. Two settings matter most.

The right solver for the pressure matrix. For standard incompressible and most compressible runs the pressure matrix is symmetric, so GAMG (or PCG) is correct. But the moment you set transonic yes in a density-based steady run (rhoSimpleFoam), the pressure equation gains a convective term and becomes asymmetric — a symmetric solver or preconditioner (GAMG, DIC, FDIC) is then invalid and will diverge. Switch to an asymmetric pair:

p
{
    solver          PBiCGStab;
    preconditioner  DILU;
    tolerance       1e-08;
    relTol          0.01;
}

Tolerances too loose. A relTol near 1 means the pressure is barely solved each outer step, so error accumulates and the coupling drifts apart. Keep relTol around 0.01–0.1 for p. More on solver choice in the GAMG and linear solvers guide.

8. Compressible and thermophysical bounds

Compressible and buoyant solvers (rhoPimpleFoam, rhoSimpleFoam, buoyantSimpleFoam) carry an energy equation and a thermophysical model with a valid temperature range. If T or p drifts outside that range — usually from an aggressive source term, a bad initial field, or a shock — the thermo lookup returns nonsense and the run crashes with a floating point exception.

Bound the fields explicitly while you stabilise the case. In fvSolution's SIMPLE/PIMPLE block:

SIMPLE
{
    rhoMin          0.1;
    rhoMax          10;
    pMinFactor      0.5;
    pMaxFactor      2;
}

And clamp temperature directly with an fvOption:

limitT
{
    type       limitTemperature;
    min        250;
    max        2500;
    selectionMode all;
}

Also relax the enthalpy/energy equation hard (h or e ≈ 0.1) at the start. See the rhoPimpleFoam setup guide for the full compressible configuration.

9. Solver-specific gotchas

• interFoam / VOF — if alpha.water goes outside [0,1], tighten MULES (nAlphaSubCycles, cAlpha) and check the interFoam VOF setup; spurious velocities at the interface often come from surface-tension/Courant coupling.
• rhoSimpleFoam (transonic) — outlet pressure waves reflecting upstream is a classic crash; add a stretched outlet extension and use a less reflective back-pressure, on top of the asymmetric solver above.
• buoyant solvers — divergence usually traces to p_rgh being over-constrained on open boundaries; use fixedFluxPressure / prghPressure rather than a plain fixedValue.

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/omega: are turbulence fields non-zero and realistic, and is the log bounding them?
7. Compressible? Confirm T stays in the thermo range and the pressure solver matches the matrix (asymmetric if transonic yes).

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;
}