The most relevant errors from a practical standpoint are solution errors^[2]. They are the difference between the exact solution of the model equations and the numerical solution. The relative solution error can be formally defined as:

(6–1)

Equation 6–1 is valid for every grid point for which the numerical solution exists. A global number can be defined by applying suitable norms, as:

(6–2)

The goal of a numerical simulation is to reduce this error below an acceptable limit.

Obviously, this is not a straightforward task, as the exact solution is not known and the error can therefore not be computed. Exceptions are simple test cases for code verification where an analytical solution is available.

Given a grid spacing , and the truncation error order of a consistent discretization scheme, , a Taylor series can be written to express the exact solution as:

(6–3)

In other words, the numerical solution converges towards the exact solution with the power of the grid spacing. Analogous definitions are available for time discretization errors.

Spatial discretization errors are the result of replacing the analytical derivatives or integrals in the exact equations by numerical approximations that have a certain truncation error. The truncation error can be obtained by inserting a Taylor series expansion of the numerical solution into the different terms of the discretized equations:

(6–4)

where is the derivative of the exact solution at a given location. An example is a central difference for a spatial derivative:

(6–5)

This formulation has a truncation error of order 2 and is therefore second-order accurate. The overall truncation error order of the spatial discretization scheme is determined by the lowest order truncation error after all terms have been discretized.

In the term of Equation 6–5, the leading term is proportional to . First-order upwind differencing of the convective terms yields truncation errors with leading term proportional to . This term then contributes to the diffusion term (numerical/false diffusion), which is most dangerous in 3D problems with grid lines not aligned to the flow direction. These schemes enhance the dissipation property of the numerical algorithm (see for example, Ferziger and Peric [141]) and are not desirable in high-quality CFD simulations.

From a practical standpoint, it is important to understand that for a first-order method, the error is reduced to 50% by a doubling of the grid resolution in each spatial direction. For a second-order method, it is reduced to 25% for the same grid refinement.

Time adds another dimension to a CFD simulation. The definition of time discretization errors is therefore similar to the definition of the spatial discretization errors. The spatial discretization usually results in a system of nonlinear algebraic equations of the form:

(6–6)

The error in the time discretization can again be obtained by a Taylor series expansion of the numerical formulation of this equation. With the example of a backward Euler integration:

(6–7)

the discretization error is:

(6–8)

The error is therefore first-order for the time derivative.

An additional complication for implicit methods comes from the inclusion of the unknown in the right hand side of Equation 6–7. In order to benefit from an implicit method, a linearization of has to be included:

(6–9)

The resulting discretized equation is therefore:

(6–10)

This constitutes an implicit formulation with first-order accuracy. A second-order time differencing is not compatible with this linearization of the right hand side, as the linearization introduced a first-order error in . In order to be able to satisfy the implicit dependency of the right hand side on the time level more closely, inner iterations (or coefficient loops) are frequently introduced:

(6–11)

where an additional iteration over the index is carried out. This equation can be reformulated as:

(6–12)

This equation can be converged completely (left hand side goes to zero) in in order to solve the original exact implicit formulation given by Equation 6–7. It is obvious that it is not necessary to converge the coefficient loop to zero, while the right hand side has a finite (first-order) error in . It can be shown that for a first-order time integration, one coefficient loop is consistent with the accuracy of the method. In a case where a second-order accurate scheme is used in the time derivative, two coefficient loops will ensure overall second-order accuracy of the method. Note, however, that this is correct only if the coefficient loops are not under-relaxed in any way.

For explicit methods, no coefficient loops are required and the time discretization error is defined solely from a Taylor series expansion.

The iteration error is similar to the coefficient loop error described above. It occurs in a case where a steady-state solution is sought from an iterative method. In most CFD codes, the iteration is carried out via a (pseudo-) timestepping scheme as given in this example, which also appears above:

(6–13)

Zero iteration error would mean that the left hand side is converged to zero, leading to the converged solution . However, in practical situations, the iterative process is stopped at a certain level, in order to reduce the numerical effort. The difference between this solution and the fully converged solution defines the iteration error.

The iteration error is usually quantified in terms of a residual or a residual norm. This can be the maximum absolute value of the right hand side, , for all grid points, or a root mean square of this quantity. In most CFD methods, the residual is nondimensionalized to enable a comparison between different applications with different scaling. However, the non-dimensionalization is different for different CFD codes, making general statements as to the required absolute level of residuals impractical. Typically, the quality of a solution is measured by the overall reduction in the residual, compared to the level at the start of the simulation.

The iteration error should be controlled with the use of the target variables. The value of the target variable can be plotted as a function of the convergence level. In case of iterative convergence, the target variable should remain constant with the convergence level. It is desirable to display the target variable in the solver monitor during the simulation.

Another numerical error is the round-off error. It results from the fact that a computer only solves the equations with a finite number of digits (around 8 for single-precision and around 16 for double-precision). Due to the limited number of digits, the computer cannot differentiate between numbers that are different by an amount below the available accuracy. For flow simulations with large-scale differences (for instance, extent of the domain vs. cell size), this can be a problem for single-precision simulations. Round-off errors are often characterized by a random behavior of the numerical solution.

The most practical method to obtain estimates for the solution error is systematic grid refinement or timestep reduction. In the following, the equations for error estimation are given for grid refinement. The same process can be used for timestep refinement.

If the asymptotic range of the convergence properties of the numerical method is reached, the difference between solutions on successively refined grids can be used as an error estimator. This allows the application of Richardson extrapolation to the solutions on the different grids (Roache [139]). In the asymptotic limit, the solution can be written as follows:

(6–14)

In this formulation, is the grid spacing (or a linear measure of it) and the are functions independent of the grid spacing. The subscript, , refers to the current level of grid resolution. Solutions on different grids are represented by different subscripts.

The assumption for the derivation of an error estimate is that the order of the numerical discretization is known. This is usually the case. Assuming a second-order accurate method, the above expansion can be written for two different grids:

(6–15)

Neglecting higher-order terms, the unknown function can be eliminated from this equation. An estimate for the exact solution is therefore:

(6–16)

The difference between the fine grid solution and the exact solution (defining the error) is therefore:

(6–17)

For an arbitrary order of accuracy, , of the underlying numerical scheme, the error is given by:

(6–18)

In order to build the difference between the solutions and , it is required that the coarse and the fine grid solution is available at the same location. In the case of a doubling of the grid density without a movement of the coarse grid nodes, all information is available on the coarse grid nodes. The application of the correction to the fine-grid solution requires an interpolation of the correction to the fine grid nodes (Roache [139]). In the case of a general grid refinement, the solutions are not available on the same physical locations. An interpolation of the solution between the different grids is then required for a direct error estimate. It has to be ensured that the interpolation error is lower than the solution error in order to avoid a contamination of the estimate.

Richardson interpolation can also be applied to integral quantities (target variables), such as lift or drag coefficients. In this case, no interpolation of the solution between grids is required.

Note that the above derivation is valid only if the underlying method has the same order of accuracy everywhere in the domain and if the coarse grid is already in the asymptotic range (the error decreases with the order of the numerical method). In addition, the method magnifies round-off and iteration errors.

The intention of the Richardson interpolation was originally to improve the solution on the fine grid. This requires an interpolation of the correction to the fine grid and introduces additional inaccuracies into the extrapolated solution, such as errors in the conservation properties of the solution. A more practical use of the Richardson extrapolation is the determination of the relative solution error, :

(6–19)

An estimate, , of this quantity can be derived from Equation 6–16:

(6–20)

It can be shown (Roache [139]) that the exact relative error and the approximation are related by:

(6–21)

Equation 6–20 can also be divided by the range of or another suitable quantity in order to prevent the error to become infinite as goes to zero.

In order to arrive at a practical error estimator, the following definitions are proposed:

Field error:

(6–22)

Maximum error:

(6–23)

RMS error:

(6–24)

Target variable error:

(6–25)

where is the defined target variable (list, drag, heat transfer coefficient, maximum temperature, mass flow, and so on).

Similar error measures can be defined for derived variables, which can be specified for each test case. Typical examples would be the total mass flow, the pressure drop, or the overall heat transfer. This will be the recommended strategy, as it avoids the interpolation of solutions between the coarse and the fine grid.

For unstructured meshes, the above considerations are valid only in cases of a global refinement of the mesh. Otherwise, the solution error will not be reduced continuously across the domain. For unstructured refinement the refinement level, , can be defined as follows:

(6–26)

where is the number of grid points and is the dimension of the problem.

It must be emphasized that these definitions do not impose an upper limit on the real error, but are estimates for the evaluation of the quality of the numerical results. Limitations of the above error estimates are:

For three-dimensional simulations, the demand that the coarse grid solution be in the asymptotic range is often hard to ensure. It is therefore required to compute the error for three different grid levels, to avoid fortuitous results. If the solution is in the asymptotic range, the following indicator should be close to constant:

(6–27)