Mesh independence and the Grid Convergence Index: how far to trust a CFD number

"The result is mesh independent" appears in most CFD reports, often supported by two runs that happen to agree. That is a weak basis for a design decision. A structured refinement study gives something better: an estimate of the discretisation error and an uncertainty band you can quote beside the result. This article walks through the standard procedure and the ways it goes wrong.

What is being estimated

A CFD solver replaces differential equations with algebraic ones on a finite mesh. The difference between the exact solution of the model equations and the solution on your mesh is the discretisation error. It shrinks as the mesh is refined, at a rate set by the order of accuracy of the numerical scheme.

It is one of several errors. Iterative error comes from stopping the solver early. Modelling error comes from the turbulence model, the boundary conditions and the simplified geometry. A mesh study addresses only the first of these. It is verification, which asks whether the equations are solved correctly, and not validation, which asks whether they are the right equations.

The three-mesh study

Build three meshes, fine (1), medium (2) and coarse (3), by refining systematically: the same topology and the same relative distribution of cells, changing only the density. For each, define a representative cell size h. For an unstructured three-dimensional mesh with N cells in a volume V, h = (V / N)^1/3.

The refinement ratio is r = h_coarse / h_fine. The widely used procedure of Celik and co-workers recommends r greater than 1.3, so that the differences between solutions stand clear of the noise. In three dimensions r = 1.3 already means about 2.2 times as many cells, and r = 2 means eight times.

Take a quantity of interest f, with values f₁, f₂ and f₃ on the three meshes. For a constant refinement ratio the observed order of accuracy is

p = ln[(f₃ − f₂) / (f₂ − f₁)] / ln r

If r differs between the two pairs, p has to be found iteratively. With p known, Richardson extrapolation estimates the value at zero cell size:

f_ext ≈ f₁ + (f₁ − f₂) / (r^p − 1)

The Grid Convergence Index

Roache's Grid Convergence Index (GCI) converts the same information into an uncertainty band on the fine-mesh result:

GCI_fine = F_s |(f₂ − f₁) / f₁| / (r^p − 1)

F_s is a safety factor. The usual values are 1.25 when three meshes are used and p is calculated from them, and 3.0 when only two meshes are available and the formal order of the scheme has to be assumed. The larger factor reflects how little two meshes tell you.

The formulae assume that the solutions lie in the asymptotic range, where error falls steadily at the rate h^p. Two checks are in common use. The observed p should be reasonably close to the formal order of the scheme, which is 2 for most industrial finite-volume settings. And the indices for the two mesh pairs should satisfy GCI₃₂ / (r^p GCI₂₁) ≈ 1, where GCI₃₂ is the same expression evaluated for the medium and coarse meshes.

A worked example

Suppose three meshes with r = 2 give a reattachment length behind a step of 5.80, 6.10 and 6.20 step heights, from coarse to fine. These are illustrative numbers, not results from a project. Then (f₃ − f₂) / (f₂ − f₁) = 3, so p = ln 3 / ln 2 ≈ 1.58. The extrapolated value is 6.20 + 0.10 / 2 = 6.25, and GCI_fine = 1.25 × (0.10 / 6.20) / 2 ≈ 1.0 %. The asymptotic check gives 1.02. The result can be reported as 6.20 step heights with a numerical uncertainty of about 1 %.

Practical pitfalls

• Iterative convergence comes first. If the solver has not converged on each mesh, the differences between meshes include iterative error and the analysis is meaningless. Residuals should have fallen by several orders of magnitude and the monitored quantity should be flat. The iterative error needs to be well below the mesh-to-mesh differences.
• Monitor the quantity you care about. Integrated quantities such as drag or pressure drop usually converge sooner than local ones such as peak wall temperature, wall shear or a reattachment point. Compute a GCI for each quantity you intend to report.
• Oscillatory convergence. If (f₃ − f₂) and (f₂ − f₁) have opposite signs, the solution is oscillating with refinement and the formula for p has no real solution. If the differences grow with refinement, the solution is diverging. In either case the meshes are not in the asymptotic range. Report the spread, refine further, and do not quote a GCI as though it were reliable.
• An implausible order. An observed p well above the formal order is a warning, not good news, because it makes the GCI unrealistically small. A common conservative practice is to cap p at the formal order.
• Non-systematic refinement. Changing only a global size while leaving inflation layers or local refinements untouched does not produce a consistent family of meshes. With wall functions there is a further conflict: the first cell must stay in the log layer, so it cannot be refined freely. See y⁺ and near-wall meshing.
• Transient runs. The time step is a discretisation parameter too, and needs refining along with the mesh.

Model choice is a separate question, covered in RANS, LES or DES? Note that in an implicitly filtered LES the mesh also sets the filter width, so the procedure above does not transfer directly.

How CFD Pro can help

CFD Pro supports clients from problem definition through to verified results and reporting, and mesh-convergence studies of this kind are part of that verification. If you need a simulation result you can defend in a design review, send us a project brief.