Since 1990, the Fluids Engineering Division of ASME has pursued activities concerning the detection, estimation and control of numerical uncertainty and/or error in computational fluid dynamics (CFD) studies. The first quality-control measures in this area were issued in 1986 (1986, “Editorial Policy Statement on Control of Numerical Accuracy,” ASME J. Fluids Eng., 108, p. 2) and revised in 1993 (1993, “Journal of Fluids Engineering Editorial Policy Statement on the Control of Numerical Accuracy,” ASME J. Fluids Eng., 115, pp. 339–340). Given the continued increase in CFD related publications, and the many significant advancements in computational techniques and computer technology, it has become necessary to revisit the issue and formulate a more detailed policy to further improve the quality of publications in this area. This brief note provides specific guidelines for prospective authors for calculation and reporting of discretization error estimates in CFD simulations where experimental data may or may not be available for comparison. The underlying perspective is that CFD-related studies will eventually aim to predict the outcome of a physical event for which experimental data is not available. It should be emphasized that the requirements outlined in this note do not preclude those already published in the previous two policy statements. It is also important to keep in mind that the procedure recommended in this note cannot possibly encompass all possible scenarios or applications.

## Preliminaries

The computer code used for an application must be fully referenced, and previous code verification studies must be briefly described or cited. The word “verification” is used in this note in its broadest sense, meaning that the computer code is capable of solving a system of coupled differential or integral equations with a properly posed set of initial and/or boundary conditions correctly, and reproduces the exact solution to these equations when sufficiently fine grid resolution (both in time and space) is employed. The formal order of accuracy in time and space for each equation solved should be also stated clearly, with proper references where this information is accessible to the readers. Before any discretization error estimation is calculated, it must be shown that iterative convergence is achieved with at least three (preferably four) orders of magnitude decrease in the normalized residuals for each equation solved. (This commonly used criterion does not always ensure adequate convergence; see Appendix) For time-dependent problems, iterative convergence at every time step should be checked, and sample convergence trends should be documented for selected, critically important, variables. A possible method for assessment of iteration errors is outlined in the Appendix.

It should also be recognized that uncertainty in inlet flow boundary conditions could be a significant contributor to the overall uncertainty. Here we recommend that the degree of sensitivity of the presented solution to small perturbations in the inlet conditions be studied and reported.

The recommended method for discretization error estimation is the Richardson extrapolation (RE) method. Since its first elegant application by its originator (1,2), this method has been studied by many authors. Its intricacies, shortcomings, and generalization have been widely investigated. A short list of references given in the bibliography (3,4,5,6,7,8,9,10,11,12,15) is selected for the direct relevance of these references to the subject and for brevity. The limitations of the RE method are well known. The local RE values of the predicted variables may not exhibit a smooth monotonic dependence on grid resolution, and in a time-dependent calculation, this nonsmooth response will also be a function of time and space. Nonetheless, it is currently the most reliable method available for the prediction of numerical uncertainty. Prospective authors can find many examples in the above references. As new and more reliable methods emerge, the present policy statement will be reassessed and modified as needed.

The Grid Convergence Method (GCI) method (and is based on RE) described herein is an acceptable and a recommended method that has been evaluated over several hundred CFD cases (19,14,2,10,7). If authors choose to use it, the method per se will not be challenged in the paper review process. If authors choose to use another method, the adequacy of their method will be judged in the review process. This policy is not meant to discourage further development of new methods. in fact, JFE encourages the development and statistically significant evaluation of alternative methods of estimation of error and uncertainty. Rather, this policy is meant to facilitate CFD publication by providing practitioners with a method that is straightforward to apply, is fairly well justified and accepted, and will avoid possible review bottlenecks, especially when the CFD paper is an application paper rather than one concerned with new CFD methodology.

## Recommended Procedure for Estimation of Discretization Error

Step 1. Define a representative cell, mesh, or grid size $h$. For example, for three-dimensional calculations,
$h=[1N∑i=1N(ΔVi)]1∕3$
1
For two dimensions,
$h=[1N∑i=1N(ΔAi)]1∕2$
2
where $ΔVi$ is the volume, $ΔAi$ is the area of the $ith$ cell, and $N$ is the total number of cells used for the computations. Equations 1,2 are to be used when integral quantities, e.g., drag coefficient, are considered. For field variables, the local cell size can be used. Clearly, if an observed global variable is used, it is then appropriate to use also an average “global” cell size. The area should be interpreted strictly according to the mesh being used, i.e., the mesh is either 2D (consisting of areas) or 3D (consisting of volumes) irrespective of the problem being solved.

Step 2. Select three significantly different sets of grids and run simulations to determine the values of key variables important to the objective of the simulation study, for example, a variable $ϕ$ critical to the conclusions being reported. It is desirable that the grid refinement factor $r=hcoarse∕hfine$ be greater than 1.3. This value of 1.3 is based on experience and not on formal derivation. The grid refinement should, however, be done systematically, that is, the refinement itself should be structured even if the grid is unstructured. Use of geometrically similar cells is preferable.

Step 3. Let $h1 and $r21=h2∕h1$, $r32=h3∕h2$, and calculate the apparent order $p$ of the method using the expression3

$p=1ln(r21)∣ln∣ε32∕ε21∣+q(p)∣$
3a
$q(p)=ln(r21p−sr32p−s)$
3b
$s=1⋅sgn(ε32∕ε21)$
3c
where $ε32=ϕ3−ϕ2$, $ε21=ϕ2−ϕ1$, and $ϕk$ denotes the solution on the $kth$ grid. Note that $q(p)=0$ for $r=const$. Equation 3 can be solved using fixed-point iteration, with the initial guess equal to the first term. The absolute value in Eq. 3a is necessary to ensure extrapolation toward $h=0$ (Celik and Karatekin (4)). Negative values of $ε32∕ε21<0$ are an indication of oscillatory convergence. If possible, the percentage occurrence of oscillatory convergence should also be reported. The agreement of the observed apparent order with the formal order of the scheme used can be taken as a good indication of the grids being in the asymptotic range; the converse should not necessarily be taken as a sign of unsatisfactory calculations. It should be noted that if either $ε32=ϕ3−ϕ2$ or $ε21=ϕ2−ϕ1$ is “very close” to zero, the above procedure does not work. This might be an indication of oscillatory convergence or, in rare situations, it may indicate that the “exact” solution has been attained. In such cases, if possible, calculations with additional grid refinement should be performed; if not, the results may be reported as such.

Step 4. Calculate the extrapolated values from
$ϕext21=(r21pϕ1−ϕ2)∕(r21p−1)$
4
similarly, calculate $ϕext32$.

Step 5. Calculate and report the following error estimates, along with the apparent order $p$.

In approximate relative error,
$ea21=∣ϕ1−ϕ2ϕ1∣$
5
In extrapolated relative error,
$eext21=∣ϕext12−ϕ1ϕext12∣$
6
In the fine-grid convergence index,
$GCIfine21=1.25ea21r21p−1$
7

Table 1 illustrates this calculation procedure for three selected grids. The data used are taken from Ref. 4, where the turbulent two-dimensional flow over a backward facing step was simulated on nonuniform structured grids with total number of cells $N1$, $N2$, and $N3$. Hence, according to Table 1, the numerical uncertainty in the fine-grid solution for the reattachment length should be reported as 2.2%. Note that this does not account for modeling errors.

Table 1

Sample calculations of discretization error

 $ϕ=dimensionless$reattachment length(with monotonicconvergence) $ϕ=axial$ velocity at$x∕H=8$⁠, $y=0.0526$$(p<1)$ $ϕ=axial$ velocity at$x∕H=8$⁠, $y=0.0526$(with oscillatoryconvergence) $N1$⁠, $N2$⁠, $N3$ 18,000, 8000, 4500 18,000, 4500, 980 18,000, 4500, 980 $r21$ 1.5 2.0 2.0 $r32$ 1.333 2.143 2.143 $ϕ1$ 6.063 10.7880 6.0042 $ϕ2$ 5.972 10.7250 5.9624 $ϕ3$ 5.863 10.6050 6.0909 $p$ 1.53 0.75 1.51 $ϕext21$ 6.1685 10.8801 6.0269 $ea21$ 1.5% 0.6% 0.7% $eext21$ 1.7% 0.9% 0.4% $GCIfine21$ 2.2% 1.1% 0.5%
 $ϕ=dimensionless$reattachment length(with monotonicconvergence) $ϕ=axial$ velocity at$x∕H=8$⁠, $y=0.0526$$(p<1)$ $ϕ=axial$ velocity at$x∕H=8$⁠, $y=0.0526$(with oscillatoryconvergence) $N1$⁠, $N2$⁠, $N3$ 18,000, 8000, 4500 18,000, 4500, 980 18,000, 4500, 980 $r21$ 1.5 2.0 2.0 $r32$ 1.333 2.143 2.143 $ϕ1$ 6.063 10.7880 6.0042 $ϕ2$ 5.972 10.7250 5.9624 $ϕ3$ 5.863 10.6050 6.0909 $p$ 1.53 0.75 1.51 $ϕext21$ 6.1685 10.8801 6.0269 $ea21$ 1.5% 0.6% 0.7% $eext21$ 1.7% 0.9% 0.4% $GCIfine21$ 2.2% 1.1% 0.5%

## Discretization Error Bars

When computed profiles of a certain variable are presented, it is recommended that numerical uncertainty be indicated by error bars on the profile, analogous to the experimental uncertainty. It is further recommended that this be done using the GCI in conjunction with an average value of $p=pave$ as a measure of the global order of accuracy. This is illustrated in Figs.1,2.

Figure 1
Figure 1
Close modal
Figure 2
Figure 2
Close modal

Figure 1 (data taken from Ref. 4) presents an axial velocity profile along the $y$-axis at an axial location of $x∕H=8.0$ for a turbulent two-dimensional backward-facing-step flow. The three sets of grids had 980, 4500, and 18000 cells, respectively. The local order of accuracy $p$ calculated from Eq. 3 ranges from 0.012 to 8.47, with a global average $pave$ of 1.49, which is a good indication of the hybrid method applied for that calculation. Oscillatory convergence occurs at 20% of the 22 points. This averaged apparent order of accuracy is used to assess the GCI index values in Eq. 7 for individual grids, which is plotted in the form of error bars, as shown in Fig. 1b. The maximum discretization uncertainty is 10%, which corresponds to $±0.35m∕s$.

Figure 2 (data taken from Ref. 16) presents an axial velocity profile along the $y$-axis at the station $x∕H=8.0$ for a laminar two-dimensional backward-facing-step flow. The Reynolds number based on step height is 230. The sets of grids used were $20×20$, $40×40$, and $80×80$, respectively. The local order of accuracy $p$ ranges from 0.1 to 3.7, with an average value of $pave=1.38$. In this figure, 80% out of 22 points exhibited oscillatory convergence. Discretization error bars are shown in Fig. 2b, along with the fine-grid solution. The maximum percentage discretization error was about 100%. This high value is relative to a velocity near zero and corresponds to a maximum uncertainty in velocity of about $±0.012m∕s$.

In the not unusual cases of noisy grid convergence, the least-squares version of GCI should be considered (13,14,17).

### Appendix: A Possible Method for Estimating Iteration Error

Following Ferziger (18,19), the iteration error can be estimated by
$εiter,in≅(ϕin+1−ϕin)λi−1$
A1
where $n$ is the iteration number and $λ1$ is the principal eigenvalue of the solution matrix of the linear system, which can be approximated from
$λi≅‖ϕin+1−ϕin‖‖ϕin−ϕin−1‖$
A2
The uncertainty $δiter$ in iteration convergence can then be estimated as
$δiter≅‖εiter,in‖λave−1$
A3
where ‖ ‖ is any appropriate norm, e.g., $L∞$ norm. Here, $λave$ is the average value of $λi$ over a reasonable number of iterations. If $λave≈1.0$, the difference between two consecutive iterations would not be a good indicator of iteration error. In order to build conservatism into these estimates, it is recommended that a limiter of $∣λ∣<2$ be applied in calculating $λave$.

It is recommended that the iteration convergence error calculated as suggested above (or in some other rational way) should be at least one order of magnitude smaller than the discretization error estimates for each calculation (for alternative methods see, e.g., Refs. 20,21).

Ismail B. Celik

Urmila Ghia

Patrick J. Roache

Christopher J. Freitas

Hugh Coleman

1.
Richardson
,
L. F.
, 1910, “
The Approximate Arithmetical Solution by Finite Differences of Physical Problems Involving Differential Equations, With an Application to the Stresses in a Masonary Dam
,”
Philos. Trans. R. Soc. London, Ser. A
0962-8428,
210
, pp.
307
357
.
2.
Richardson
,
L. F.
, and
Gaunt
,
J. A.
, 1927, “
The Deferred Approach to the Limit
,”
Philos. Trans. R. Soc. London, Ser. A
0962-8428,
226
, pp.
299
361
.
3.
Celik
,
I.
,
Chen
,
C. J.
,
Roache
,
P. J.
, and
Scheurer
,
G.
, Editors. (1993), “
Quantification of Uncertainty in Computational Fluid Dynamics
,”
ASME Fluids Engineering Division Summer Meeting
,
Washington, DC
, Jun. 20–24, ASME Publ. No. FED-Vol.
158
.
4.
Celik
,
I.
, and
Karatekin
,
O.
, 1997, “
Numerical Experiments on Application of Richardson Extrapolation With Nonuniform Grids
,”
ASME J. Fluids Eng.
0098-2202,
119
, pp.
584
590
.
5.
Eca
,
L.
, and
Hoekstra
,
M.
, 2002, “
An Evaluation of Verification Procedures for CFD Applications
,”
24th Symposium on Naval Hydrodynamics
,
Fukuoka, Japan
, Jul. 8–13.
6.
Freitas
,
C. J.
, 1993, “
Journal of Fluids Engineering Editorial Policy Statement on the Control of Numerical Accuracy
,”
ASME J. Fluids Eng.
0098-2202,
115
, pp.
339
340
.
7.
Roache
,
P. J.
, 1993, “
A Method for Uniform Reporting of Grid Refinement Studies
,”
Proceedings of Quantification of Uncertainty in Computation Fluid Dynamics
, Edited by
Celik
, et al.
,
ASME Fluids Engineering Division Spring Meeting
,
Washington, D.C.
, June 23–24, ASME Publ. No. FED-Vol.
158
.
8.
Roache
,
P. J.
, 1998, “
Verification and Validation in Computational Science and Engineering
,”
Hermosa Publishers
,
Albuquerque, NM
.
9.
Stern
,
F.
,
Wilson
,
R. V.
,
Coleman
,
H. W.
, and
Paterson
,
E. G.
, 2001, “
Comprehensive Approach to Verification and Validation of CFD Simulations—Part 1: Methodology and Procedures
,”
ASME J. Fluids Eng.
0098-2202,
123
, pp.
793
802
.
10.
,
B. L.
,
Rearden
,
B. T.
,
Hopper
,
C. M.
,
Wagschal
,
J. J.
, and
Parks
,
C. V.
, 2004, “
Sensitivity- and Uncertainty-Based Criticality Safety Validation Techniques
,”
Nucl. Sci. Eng.
0029-5639,
146
, pp.
340
366
.
11.
DeVolder
,
B.
,
Glimm
,
J.
,
Grove
,
J. W.
,
Kang
,
Y.
,
Lee
,
Y.
,
Pao
,
K.
,
Sharp
,
D. H.
, and
Ye
,
K.
, 2002, “
Uncertainty Quantification for Multiscale Simulations
,”
ASME J. Fluids Eng.
0098-2202,
124
(
1
), pp.
29
41
.
12.
Oberkampf
,
W. L.
,
Trucano
,
T. G.
, and
Hirsch
,
C.
, 2003, “
Verification, Validation, and Predictive Capability in Computational Engineering and Physics
,” Sandia Report No. SAND2003-3769.
13.
Eça
,
L.
,
Hoekstra
,
M.
, and
Roache
,
P. J.
, 2005, “
Verification of Calculations: An Overview of the Lisbon Workshop
,”
AIAA Computational Fluid Dynamics Conference
,
, Jun., AIAA Paper No. 4728.
14.
Eça
,
L.
,
Hoekstra
,
M.
, and
Roache
,
P. J.
, 2007, “
Verification of Calculations: an Overview of the 2nd Lisbon Workshop
,”
Second Workshop on CFD Uncertainty Analysis, AIAA Computational Fluid Dynamics Conference
,
Miami, FL
, Jun., AIAA Paper No. 2007-4089.
15.
Roache
,
P. J.
,
Ghia
,
K. N.
, and
White
,
F. M.
, 1986, “
Editorial Policy Statement on Control of Numerical Accuracy
,”
ASME J. Fluids Eng.
0098-2202,
108
, p.
2
.
16.
Celik
,
I.
, and
,
A.
, Jr.
, 2003, “
Verification and Validation of DREAM Code
,” Mechanical and Aerospace Engineering Department, Report No. MAE̱IC03/TR103.
17.
Roache
,
P. J.
, 2003, “
Conservatism of the GCI in Finite Volume Computations on Steady State Fluid Flow and Heat Transfer
,”
ASME J. Fluids Eng.
0098-2202,
125
(
4
), pp.
731
732
.
18.
Ferziger
,
J. H.
, 1989, “
Estimation and Reduction of Numerical Error
,”
ASME Winter Annual Meeting
,
San Fransisco, CA
, Dec.
19.
Ferziger
,
J. H.
, and
Peric
,
M.
, 1996, “
Further Discussion of Numerical Errors in CFD
,”
Int. J. Numer. Methods Fluids
0271-2091,
23
, pp.
1263
1274
.
20.
Eça
,
L.
, and
Hoekstra
,
M.
, 2006, “
On the Influence of the Iterative Error in the Numerical Uncertainty of Ship Viscous Flow Calculations
,”
26th Symposium on Naval Hydrodynamics
,
Rome, Italy
, Sept. 17–22.
21.
Stern
,
F.
,
Wilson
,
R.
, and
Shao
,
J.
, 2006, “
Quantitative V&V of CFD Simulations and Certification of CFD Codes
,”
Int. J. Numer. Methods Fluids
0271-2091,
50
, pp.
1335
1355
.