When a good mesh gives the wrong answer: OpenFOAM vs code_saturne on the backward-facing step
How much can a mesh go wrong before a RANS solver stops converging? And when it does, does a gradient limiter really fix it, or does it just hide the problem? We ran about 95 simulations of the same 2D test case in OpenFOAM and code_saturne on cloudHPC to find out, with four levels of numerical stabilisation. The whole study cost less than 20 โฌ.
The test case
We used the backward-facing step of Driver & Seegmiller, as set up in the NASA Turbulence Modeling Resource: step height H, Re_H = 36,000, SpalartโAllmaras model, steady solution. The quantity we track is the reattachment length x_r/H. The experiment gives 6.26 ยฑ 0.10. CFL3D, solving the same SA model on a very fine grid, gives 6.06.
Both codes ran the same meshes, converted from the same gmsh file, with the same boundary conditions and relaxation. We tested 18 meshes in four families:
| Setting | OpenFOAM | code_saturne |
|---|---|---|
| Algorithm | steady SIMPLE (incompressibleFluid) | steady algorithm |
| Turbulence | SpalartโAllmaras, y+ โ 1 | SpalartโAllmaras, y+ โ 1 |
| Inlet | U = 1, ฮฝฬ = 3ฮฝ | U = 1, ฮฝฬ = 3ฮฝ |
| Relaxation | U, ฮฝฬ 0.7 ยท p 0.3 | U, ฮฝฬ 0.7 ยท p 0.3 |
| Pressure solver | GAMG | algebraic multigrid |
| Max iterations | 20,000 | 5,000โ20,000 |
- Refinement: the NASA structured grid and two refinements (20k to 320k cells).
- Distortion: the same grid with interior nodes randomly moved by up to 45ยฐ.
- Hybrid: structured wall layers plus an unstructured quad or triangle core.
- Hybrid variants: different size jumps between the last wall layer and the core, first-cell heights (y+ โ 1, 5, 30) and growth ratios (1.12, 1.2, 1.4).
Result 1: on good meshes the two codes agree
On the structured and distorted meshes both codes converged with plain second-order numerics and no limiters. Distortion up to 45ยฐ changed almost nothing.
| Mesh | Cells | x_r/H OpenFOAM | x_r/H code_saturne |
|---|---|---|---|
| NASA grid (T3) | 20k | 6.04 | 6.01 |
| T3, refined 2ร | 80k | 6.03 | 5.99 |
| T3, refined 4ร | 319k | 6.05 | 5.98 |
| Distorted 15ยฐ / 30ยฐ | 80k | 6.03 | 5.99 |
| Distorted 45ยฐ | 80k | 5.97 | 5.99 |
Both codes sit within 2.5% of CFL3D (6.06) and 3โ6% below the experiment. That gap is the SA model, not the numerics: every code solving SA on a good mesh lands there.
The velocity profiles of the two codes are indistinguishable and follow CFL3D through the shear layer. Near the wall at x/H = 1 and 4 all three codes predict stronger reverse flow than measured: a known SA behaviour.
Result 2: hybrid meshes are where things break
The critical spot is the transition from the last wall layer to the unstructured core. Without limiters:
- OpenFOAM diverged on 6 hybrid meshes. With a 40-layer boundary layer it diverged even when the size jump to the core was only about 3.7.
- code_saturne diverged on the two triangle-core meshes and converged on 4 of the 6 hybrid variants.
- Growth ratio 1.2 was the only hybrid variant on which both codes converged.
| Hybrid mesh | What changes | OpenFOAM | code_saturne |
|---|---|---|---|
| G1 quad / tri | core up to 4โ5 H far from the step | diverges | no convergence / diverges |
| G1b quad | core capped at 0.3 H (jump โ 10) | diverges | converges (5.97) |
| G1b tri | same, triangle core | diverges | diverges |
| Jump โ 3.7 | core 0.1 H, 297k cells | diverges | converges (5.97) |
| Jump โ 22 | core 0.6 H | diverges | converges (5.91) |
| y+ โ 5 | first cell 3.5e-3 | stalls | stalls |
| y+ โ 30 | wall functions | stalls (5.78) | stalls (5.70) |
| Growth 1.2 | 28 layers | converges (5.97) | converges (5.96) |
| Growth 1.4 | 17 layers | stalls | converges (5.91) |
The checkMesh aspect ratio does not explain the failures. The NASA grid reaches about 15,700 and converges; the hybrids stay around 1,500 and do not. Wall functions at y+ โ 30 shorten the bubble by 4โ5% on top of the convergence problem.
Result 3: limiters stop the divergence, not the problem
We reran the five failed cases (four in OpenFOAM, one in code_saturne) with four stabilisation levels. Each level changes one thing on top of L1, so the effect of each knob can be read on its own.
| Level | OpenFOAM | code_saturne | What it tests |
|---|---|---|---|
| L1 | cellLimited gradients, limited corrected 0.5, 1 non-orthogonal corrector | slope test + cell gradient limiter | gradient limiting alone |
| L2 | L1 + bounded Gauss upwind on ฮฝฬ, limited 0.33, 2 correctors | L1 + upwind, no reconstruction on ฮฝฬ | first-order turbulence transport |
| L3 | L1 + limitedLinear 1 on ฮฝฬ | L1 + 50% SOLU / 50% upwind blend on ฮฝฬ | a bounded second-order alternative to L2 |
| L4 | L1 + PCG/DIC pressure solver instead of GAMG | L1 + Jacobi-preconditioned CG on pressure | whether the pressure solver drives the failure |
| Setting | Cases converged | x_r/H range |
|---|---|---|
| L1 | 0 of 5 (all stall, bounded) | 6.00โ6.05 |
| L2 | 1 of 5 | 6.04โ6.48 |
| Failed case | L1 | L2 | L3 | L4 |
|---|---|---|---|---|
| OpenFOAM, hybrid quad | stalls, 6.05 | stalls, 6.11 | stalls, 6.05 | stalls, 6.05 |
| OpenFOAM, hybrid tri | stalls, 6.00 | stalls, 6.48 | stalls, 6.08 | stalls, 6.00 |
| OpenFOAM, jump โ 3.7 | stalls, 6.00 | stalls, 6.04 | stalls, 6.00 | stalls, 6.00 |
| OpenFOAM, jump โ 22 | stalls, 6.04 | converges, 6.28 | stalls, 6.09 | stalls, 6.04 |
| code_saturne, hybrid tri | stalls, 6.04 | stalls, 6.37 | stalls, 6.20 | stalls like L1 (1,900 it, 6ร slower) |
Values are x_r/H. โStallsโ means bounded but not converged after 20,000 iterations (OpenFOAM) or 5,000 (code_saturne). Good meshes give 5.91โ6.05.
What the two new levels tell us:
- L4: the pressure solver is not the culprit. Swapping GAMG for a Krylov solver changed nothing: same stall, same x_r as L1 to the second decimal. We had suspected multigrid struggling with the high-aspect-ratio wall layers. That hypothesis is now ruled out, and the mesh transition is the remaining suspect.
- L3: a bounded second-order scheme is no free lunch either. limitedLinear on ฮฝฬ did not converge anywhere, and its ฮฝฬ residual stalled about four orders of magnitude higher than with L1 (2e-4 to 9e-4 instead of 1e-8 to 1e-6). x_r drifted up by up to 1.5% in OpenFOAM and 3% in code_saturne: less than L2, in the same direction.
- Rule of thumb: every step away from pure second-order transport of ฮฝฬ (L3, then L2) lengthens the bubble. On these meshes the limiter level is a numerical parameter you must verify, just like the grid.
L1 is the useful one. It removed every divergence, and although the residuals stalled, x_r stayed in the same band as the good meshes.
L2 is the trap. The upwind scheme on the turbulence variable lengthened the bubble by up to 8%. In the one case where it converged, it gave x_r/H = 6.28: almost exactly the experiment. That is not accuracy. It is a numerical error that happens to compensate the known shortfall of the SA model.
Result 4: a perfect mesh that gives the wrong answer
Our first mesh family was structured and perfectly orthogonal, and passed checkMesh with zero non-orthogonality. The coarsest level gave x_r/H = 6.06: a great result. Refining it made things worse at every level, down to 2.10 in OpenFOAM and 4.88 in code_saturne.
| Level | Cells | x_r/H OpenFOAM | x_r/H code_saturne |
|---|---|---|---|
| Coarse | 18k | 6.06 | 5.94 |
| Medium | 74k | 4.36 | 5.41 |
| Fine | 296k | 2.10 | 4.88 |
The cause: the generator kept wall-type clustering along the step-lip height all the way downstream, where there is no wall but a free shear layer. Every refinement added more thin cells there and distorted the eddy viscosity.
Three lessons:
- Mesh-quality checks look at cell shapes, not at whether the cells suit the flow.
- One mesh level is never enough. Only the refinement study exposed the problem.
- Both codes were wrong in the same direction. Comparing solvers does not replace mesh verification.
Takeaways
- Keep the transition from the wall layers to the core gradual, with a growth ratio of 1.2 or less.
- Use gradient limiting (L1) as a safety net, and check that x_r does not move when you switch it on. If L1 stalls, fix the mesh: in our tests neither a different convection scheme (L2, L3) nor a different pressure solver (L4) did.
- Avoid first-order upwind on the turbulence variable as a default. It converges more easily and moves the answer.
- Separate verification from validation. Matching the experiment better with a cheaper scheme is not a better simulation.
All runs used hypercore machines on cloudHPC with 2 to 8 vCPU. A typical run cost between 0.04 and 0.30 โฌ.
Run it yourself on cloudHPC
Each case is a small 2D model, so it fits on the smallest machines:
- Pack the OpenFOAM case folder (or the code_saturne study folder) as a tar.gz and upload it to your cloudHPC storage.
- Launch it on a hypercore machine: 2 vCPU are enough up to about 100k cells, 8 vCPU for about 300k.
- Download the results archive and compute x_r from the zero crossing of the wall shear stress on the lower wall.
A full mesh study like this one (refinement, distortion and limiter checks) costs a few euros and runs in parallel in an afternoon.
References
CloudHPC is a HPC provider to run engineering simulations on the cloud. CloudHPC provides from 1 to 224 vCPUs for each process in several configuration of HPC infrastructure - both multi-thread and multi-core. Current software ranges includes several CAE, CFD, FEA, FEM software among which OpenFOAM, FDS, Blender and several others.
New users benefit of a FREE trial of 300 vCPU/Hours to be used on the platform in order to test the platform, all each features and verify if it is suitable for their needs
