CFD Mesh Independence Study (Grid Independence Study)

The principle is the same. Both are tasks to confirm the state where "the solution hardly changes even if the mesh is made finer." Estimating the solution at mesh size $h \to 0$ using Richardson extrapolation and quantifying the grid error with GCI (Grid Convergence Index) are also common.

There are a few different points to note. In CFD, there are three additional conditions. First, whether y+ is consistent with the turbulence model. If this is not correct, the solution will not converge to the correct one no matter how fine the mesh is. Second, because CFD has strong nonlinearity in the Navier-Stokes equations, the observed convergence order $p$ is often lower than the theoretical value ($p=2$ for a second-order accurate scheme). Third, if the residual for iterative convergence is not sufficiently reduced, discretization error and iterative error mix, making GCI meaningless.

Yes, there is. The paper by Celik et al. (2008) "Procedure for Estimation and Reporting of Uncertainty due to Discretization in CFD Applications" (Journal of Fluids Engineering) is the de facto standard procedure for applying GCI in CFD. It has also become a review criterion for ASME Journal submissions. If you master this 5-step procedure, you can use it in both practical work and papers.

If we Taylor expand the discrete solution $f(h)$, its relationship with the exact solution $f_{\text{exact}}$ is as follows:

Here, $h$ is the representative cell size, $p$ is the accuracy order of the discretization scheme, and $C$ is an unknown constant. For example, with a theoretically second-order accurate central difference, $p=2$, so halving the mesh reduces the error to about $1/4$. The basic idea of Richardson extrapolation is to estimate $f_{\text{exact}}$ using two or more solutions with different mesh sizes based on this relationship.

Exactly. With three mesh levels (fine $h_1$, medium $h_2$, coarse $h_3$, where $h_1 < h_2 < h_3$) and solutions $f_1, f_2, f_3$, we can estimate the observed convergence order $p$ by solving simultaneous equations. With only two levels, we have to assume $p$, which leads to a subjective judgment like "the change was small when the mesh was refined, so it's OK."

The formula itself is the same. The GCI between the fine and medium meshes is:

The meaning of each variable is as follows:

• $e_a^{21} = (f_1 - f_2)/f_1$: Relative difference between the fine mesh solution $f_1$ and the medium mesh solution $f_2$
• $r_{21} = h_2 / h_1$: Grid refinement ratio ($r > 1$)
• $p$: Observed convergence order (calculated using the formula mentioned later)
• $F_s$: Safety factor. $F_s = 1.25$ when using three-level extrapolation; $F_s = 3.0$ when assuming $p$ with two levels

Yes. $F_s = 3$ is a penalty for the "uncertainty of assuming $p$." Just doing one extra level reduces the safety factor to $1/2.4$, so there's no reason not to do three levels. For example, if the GCI for the drag coefficient $C_d$ in external automotive aerodynamics is below 1.5%, it can be considered accurate enough for comparison with wind tunnel tests.

Let me explain Celik et al.'s (2008, JFE 130(7)) 5-step procedure.

Step 1. Definition of representative cell size $h$

For 3D:

$N$ is the total number of cells, and $\Delta V_i$ is the volume of each cell. For 2D, use the square root of the area.

Step 2. Obtain solutions with three mesh levels

The ideal refinement ratio $r = h_{\text{coarse}}/h_{\text{fine}}$ is $r \geq 1.3$. $r = 2$ (8 times the number of cells) is best, but if computational cost is tight, $r = 1.3$ to $1.5$ (2-3 times the cells) is OK. Obtain three solutions $f_1, f_2, f_3$.

Step 3. Determine the observed convergence order $p$

Here, $\varepsilon_{21} = f_2 - f_1$, $\varepsilon_{32} = f_3 - f_2$. If $r_{21}$ and $r_{32}$ are equal (geometric refinement), then $q=0$, and the simplified formula can be used:

Step 4. Obtain the Richardson extrapolated value $f_{\text{ext}}^{21}$

Step 5. Report the error estimates

• Relative approximate error: $e_a^{21} = |(f_1 - f_2)/f_1|$
• Relative extrapolated error: $e_{\text{ext}}^{21} = |(f_{\text{ext}}^{21} - f_1)/f_{\text{ext}}^{21}|$
• $\text{GCI}_{\text{fine}}^{21} = \frac{1.25 \, e_a^{21}}{r_{21}^{p} - 1}$

Solve it iteratively. Start with an initial value for $p$ (e.g., from the simplified formula), update $q(p)$, recalculate $p$, and repeat a few times until convergence. In Python, you can use scipy.optimize.fsolve to solve it in one go. For geometric refinement ($r_{21} = r_{32}$), $q=0$, so no iteration is needed; you can use the simplified version directly.

Good question. y+ is the dimensionless distance of the first cell layer from the wall, defined as:

$y$ is the distance from the wall to the center of the first cell layer, $u_\tau$ is the friction velocity, and $\nu$ is the kinematic viscosity. The key point is that each turbulence model has a specific required range for y+:

• Using wall functions (standard $k$-$\varepsilon$, etc.): $y^+ = 30\text{ to }300$
• Wall-resolved ($k$-$\omega$ SST, SA, LES): $y^+ < 1$ (ideal is $y^+ \approx 1$)
• Blended wall treatment (SST's automatic wall treatment): $y^+ < 5$ is desirable

Wall functions assume the logarithmic law (region where $y^+ > 30$), so $y^+ = 5$ is right in the middle of the buffer layer, a "physically inconsistent" zone where neither the logarithmic law nor the linear law of the viscous sublayer holds. There, the wall shear stress will be significantly wrong. As the mesh is refined, y+ decreases and moves further outside the applicable range of wall functions, causing the solution to show divergent behavior or non-monotonic convergence. As a result, the premise of GCI (monotonic convergence) breaks down.

Exactly. A common pattern in practice is "calculating with wall functions but y+ is around 3 on some walls, not noticing it, doing a mesh study, and struggling with 'it won't converge!'" It's like saying "that's strange" while measuring without properly inserting a thermometer.

Celik's procedure uses the global average. For 3D, it's the cube root of the average volume of all cells:

Practically, "if you know the total number of cells $N$, you can approximately calculate it from the total volume of the computational domain $V_{\text{total}}$ as $h = (V_{\text{total}}/N)^{1/3}$," so you can calculate it quickly if the solver output includes the cell count.

Sharp observation. Strictly speaking, the global $h$ is appropriate only when the mesh in the region on which the evaluation quantity depends is uniformly refined. For example, if only the near-wall region is refined while the far field remains coarse, the effective refinement ratio $r$ changes. The best practice is to refine the entire domain similarly (uniformly increasing the total cell count while maintaining the geometry's proportions). In Fluent, you can scale all parameters in Mesh > Controls > Sizing at once.

Ideally yes, but in practice, $p$ often varies widely, around $p = 0.5$ to $3.0$. For nonlinear problems, $p < 2$ is common. As a guideline:

• $p \approx 1.5$ to $2.5$: Normal. Can be judged to be in the asymptotic range for second-order accuracy.
• $p < 0.5$: Mesh is still too coarse to be in the asymptotic range, or there is a bug.
• $p > 3.0$: Possibly due to accidental cancellation. Should add another level to confirm.
• $p < 0$ ($\varepsilon_{32}/\varepsilon_{21} < 0$): Non-monotonic convergence. The standard procedure cannot be applied.

Absolutely no problem. The advantage of Celik's procedure is that by using the observed value for $p$ instead of fixing it to a theoretical value, GCI becomes a reliable error estimate even for nonlinear problems. It gives a much more honest estimate than calculating by assuming $p=2$.