[CFD Q&A] When simpleFoam Pressure Residual Fails to Converge

simpleFoam ate half a day last week. The pressure residual would not go under 1e-3. It climbed. The mesh was the first suspect. The knob that mattered was the relaxation factor.

Push the αp slider past 0.6 and the pressure-residual curve bends upward again near iteration 30.

Park it near 0.3 and the same curve keeps falling out to iteration 200. Time=80 in the log below, initial residual 0.891, sits in that same bend.

Problem Situation#

OpenFOAM v2312. A 2D rectangular channel, inlet to outlet. Turbulence model, standard k-ε. Solver, steady simpleFoam.

The residual fell for the first 20 to 30 iterations. Then only the pressure jumped.

Time = 50
 
smoothSolver:  Solving for Ux, Initial residual = 0.0023, Final residual = 1.2e-5, No Iterations 3
smoothSolver:  Solving for Uy, Initial residual = 0.0019, Final residual = 9.8e-6, No Iterations 3
GAMG:  Solving for p, Initial residual = 0.412, Final residual = 0.038, No Iterations 8
time step continuity errors : sum local = 4.2e-5, global = 1.1e-5, cumulative = 3.4e-4
...
Time = 80
 
GAMG:  Solving for p, Initial residual = 0.891, Final residual = 0.21, No Iterations 20

Velocity settled. Pressure did not. That is the face of a pressure-velocity coupling that is rocking.

Root Cause Analysis Process#

Step 1: Check Mesh Quality#

checkMesh came first.

checkMesh -case .

Max non-orthogonality 32 degrees, max skewness 0.8. Too quiet a pair of numbers to blame the mesh.

Step 2: Review Boundary Conditions#

Inlet pressure was zeroGradient, outlet fixedValue 0. A common pair for a channel. Nothing wrong there.

The initial p field still held a large value near the inlet. Leftovers from a funkySetFields trial that was supposed to be gone.

Step 3: Check fvSolution#

// fvSolution (original settings)
relaxationFactors
{
    fields
    {
        p               0.7;
    }
    equations
    {
        U               0.9;
        k               0.7;
        epsilon         0.7;
    }
}

U=0.9 is a fast setting. Raise the velocity relaxation and the residual looks quicker. If the pressure update cannot keep up, the loop rocks.

Two causes stacked. The relaxation was high. Garbage remained in the initial p. How the pressure correction stays tied to the velocity field is written out in SIMPLE when pressure starts to checkerboard.

Solution#

Lowering Relaxation Factors#

// fvSolution (after modification)
relaxationFactors
{
    fields
    {
        p               0.3;
    }
    equations
    {
        U               0.7;
        k               0.5;
        epsilon         0.5;
    }
}

The first run starts low. The factors go up only after the residual has actually sat down.

Initializing the p Field#

# Cleanly initialize the p file in the 0 directory
cp 0.orig/p 0/p

Opening 0/p and setting internalField uniform 0; does the same job.

// 0/p
dimensions      [0 2 -2 0 0 0 0];
internalField   uniform 0;
 
boundaryField
{
    inlet
    {
        type            zeroGradient;
    }
    outlet
    {
        type            fixedValue;
        value           uniform 0;
    }
    walls
    {
        type            zeroGradient;
    }
}

Verification#

The same case, run again, changes the log to this.

Time = 50
 
smoothSolver:  Solving for Ux, Initial residual = 0.0018, Final residual = 8.3e-6, No Iterations 3
GAMG:  Solving for p, Initial residual = 0.089, Final residual = 0.0041, No Iterations 5
time step continuity errors : sum local = 2.1e-6, global = 4.3e-7, cumulative = 8.9e-6
 
Time = 200
 
GAMG:  Solving for p, Initial residual = 3.2e-4, Final residual = 1.1e-5, No Iterations 3

Near iteration 200 the pressure residual had come down to about 1e-4. Inner iterations were back in single digits.

The Ratios Sitting in the Log#

Before the fix, Time=50 dropped Ux from 0.0023 to 1.2e-5 in three inner iterations, about 192 times. Uy went from 0.0019 to 9.8e-6, about 194 times. Pressure only went from 0.412 to 0.038, about 11 times in eight iterations. Velocity fell nearly 200-fold. Pressure managed about 11.

Time=80 is thirty steps on. The initial pressure residual rose from 0.412 to 0.891, about 2.2 times. Inner iterations grew from 8 to 20, a factor of 2.5. The cut, 0.891 divided by 0.21, shrank to about 4.2. More iterations, a smaller drop.

After the fix, Time=50 pressure falls from 0.089 to 0.0041, about 22 times in five iterations. Time=200 goes from 3.2e-4 to 1.1e-5, about 29 times in three. Cumulative continuity shrank from 3.4e-4 to 8.9e-6, about 38 times. p at 0.3 is about 43% of the old 0.7. U moved from 0.9 to 0.7, and k and epsilon from 0.7 to 0.5. Cutting the pressure correction below half brought the inner count down from 20 into single digits.

Set the final residual at Time=200, 1.1e-5, next to the old Time=80 final of 0.21. Divide 0.21 by 1.1e-5 and the ratio is about 19100. Same channel, same mesh. Only the pressure relaxation and the initial field changed. That ratio is the number to read before regenerating a mesh.

The Knob to Turn Before Opening the Mesh#

Non-orthogonality of 32 degrees and skewness of 0.8 are not a reason to rebuild the mesh first. The value to lower is the p relaxation. At 0.7 the pressure correction could not keep up with velocity at 0.9. The first run uses p 0.3 and U 0.7. The factors go up only after the Time=200 residual has sat near 1e-4.

The funkySetFields leftover in the initial p wrecked the same run. Leave it, and a lower relaxation still carries a junk gradient in from the inlet. A lower relaxation is not automatically the same physics. The Rhie–Chow case where one relaxation change moved the converged solution is where to check whether 0.7 down to 0.3 stays in the residual curve or leaks into the final field.