0
Announcements

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

J. Fluids Eng 130(7), 078001 (Jul 22, 2008) (4 pages) doi:10.1115/1.2960953 History: Published July 22, 2008

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.

FIGURES IN THIS ARTICLE
<>

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-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 1. Define a representative cell, mesh, or grid size h. For example, for three-dimensional calculations,Display Formula

h=[1Ni=1N(ΔVi)]13
(1)
For two dimensions,Display Formula
h=[1Ni=1N(ΔAi)]12
(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=hcoarsehfine 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=h2h1, r32=h3h2, and calculate the apparent order p of the method using the expressionDisplay Formula

p=1ln(r21)lnε32ε21+q(p)
(3a)
Display Formula
q(p)=ln(r21psr32ps)
(3b)
Display Formula
s=1sgn(ε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. 3 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 fromDisplay Formula

ϕext21=(r21pϕ1ϕ2)(r21p1)
(4)
similarly, calculate ϕext32.

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

In approximate relative error,Display Formula

ea21=ϕ1ϕ2ϕ1
(5)

In extrapolated relative error,Display Formula

eext21=ϕext12ϕ1ϕext12
(6)

In the fine-grid convergence index,Display Formula

GCIfine21=1.25ea21r21p1
(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.

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.  12.

Figure 1 (data taken from Ref. 4) presents an axial velocity profile along the y-axis at an axial location of xH=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. 1. The maximum discretization uncertainty is 10%, which corresponds to ±0.35ms.

Figure 2 (data taken from Ref. 16) presents an axial velocity profile along the y-axis at the station xH=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. 2, 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.012ms.

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

Following Ferziger (18-19), the iteration error can be estimated byDisplay Formula

εiter,in(ϕin+1ϕin)λi1
(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 fromDisplay Formula
λiϕin+1ϕinϕinϕin1
(A2)
The uncertainty δiter in iteration convergence can then be estimated asDisplay Formula
δiterεiter,inλave1
(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 λave1.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

Peter E. Raad

Copyright © 2008 by American Society of Mechanical Engineers
View article in PDF format.

References

Figures

Grahic Jump Location
Figure 1

(a) Axial velocity profiles for a two-dimensional turbulent backward-facing-step flow calculation (16); (b) Fine-grid solution, with discretization error bars computed using Eq. 7.

Grahic Jump Location
Figure 2

(a) Axial velocity profiles for a two-dimensional laminar backward-facing-step flow calculation (16); (b) Fine-grid solution, with discretization error bars computed using Eq. 7

Tables

Table Grahic Jump Location
Table 1
Sample calculations of discretization error

Errata

Discussions

Some tools below are only available to our subscribers or users with an online account.

Related Content

Customize your page view by dragging and repositioning the boxes below.

Related Journal Articles
Related eBook Content
Topic Collections

Sorry! You do not have access to this content. For assistance or to subscribe, please contact us:

  • TELEPHONE: 1-800-843-2763 (Toll-free in the USA)
  • EMAIL: asmedigitalcollection@asme.org
Sign In