Procedure for Estimation and Reporting of Uncertainty Due to Discretization in CFD Applications

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.

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<h2<h3$ and $r21=h2∕h1$⁠, $r32=h3∕h2$⁠, and calculate the apparent order $p$ of the method using the expression

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 5. Calculate and report the following error estimates, along with the apparent order $p$⁠.

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.

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 (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).

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

Peter E. Raad