This article discusses numerical errors in unsteady flow simulations, which may include round-off, statistical, iterative, and time and space discretization errors. The estimation of iterative and discretization errors and the influence of the initial condition on unsteady flows that become periodic are discussed. In this latter case, the goal is to determine the simulation time required to reduce the influence of the initial condition to negligible levels. Two one-dimensional, unsteady manufactured solutions are used to illustrate the interference between the different types of numerical errors. One solution is periodic and the other includes a transient region before it reaches a steady-state. The results show that for a selected grid and time-step, statistical convergence of the periodic solution may be achieved at significant lower error levels than those of iterative and discretization errors. However, statistical convergence deteriorates when iterative convergence criteria become less demanding, grids are refined, and Courant number increased.For statistically converged solutions of the periodic flow and for the transient solution, iterative convergence criteria required to obtain a negligible influence of the iterative error when compared to the discretization error are more strict than typical values found in the open literature. More demanding criteria are required when the grid is refined and/or the Courant number is increased. When the numerical error is dominated by the iterative error, it is pointless to refine the grid and/or reduce the time-step. For solutions with a numerical error dominated by the discretization error, three different techniques are applied to illustrate how the discretization uncertainty can be estimated, using grid/time refinement studies: three data points at a fixed Courant number; five data points involving three time steps for the same grid and three grids for the same time-step; five data points including at least two grids and two time steps. The latter two techniques distinguish between space and time convergence, whereas the first one combines the effect of the two discretization errors.

Introduction

The use of computational fluid dynamics (CFD) for the simulation of flows of engineering interest is no longer restricted to (statistically) steady flows. In fact, the development of computer power and numerical techniques has led to several mathematical models that (partially or not) solve turbulent fluctuations in flows around complex geometries at high Reynolds numbers [1]. This enhancement of the quality of the modeling and numerical solution techniques must be accompanied by a careful assessment of the numerical uncertainty, which is more difficult to handle in unsteady flows than in steady flows.

In the last two decades, a large effort has been dedicated to verification (assessment of numerical accuracy) and validation (assessment of modeling accuracy) with several reference works published [24], and a dedicated symposium was organized every year by ASME since 2012. Furthermore, several associations (ASME [5,6], AIAA [7], ITTC [8]) have published guides and standards for assessing the modeling errors and numerical, experimental, and input uncertainties, as for example, the ASME V&V 20 standard dedicated to verification and validation in computational fluid dynamics and heat transfer [6].

In this paper, we focus on solution verification for unsteady flow simulations, i.e., on the estimation of numerical errors for problems that include variations in both space and time. Clearly, the information available in the literature, as for example, Refs. [24] and [7], already covers different contributions to the numerical error of unsteady flow simulations, like round-off, iterative, and discretization errors. However, many of the available examples refer to steady flow situations, where iterative errors are easier to handle and discretization errors contain only space contributions. In implicit time integration schemes (the main focus of this work), a nonlinear problem must be solved at each time-step and so numerical accuracy of the following time steps will be influenced by the selected iterative convergence criterion. Discretization errors include space and time contributions, making the estimation of the discretization error more challenging than in steady flow simulations.

Moreover, unsteady flow simulations require an initial condition that can also affect the accuracy of the simulations. In principle, any inexactness present in the initial condition should be handled as an input uncertainty, which also applies to the boundary conditions, as described for instance in Ref. [6]. However, there are flows tending to periodic behavior, as for example, vortex shedding from a circular cylinder, which are often simulated with an incorrect initial condition. In those cases, it is important to know how many cycles need to be simulated until the solution becomes numerically periodic, i.e., until it statistically converges. If statistical convergence is achieved, it is even more important to know whether or not the final result is affected by the error in the initial condition.

Taking advantage of the method of manufactured solutions [9], it is possible to address iterative, statistical, and discretization errors2 in unsteady flow simulations, using the same framework as that presented in Ref. [10] for steady flow simulations. An unsteady, one-dimensional, convection–diffusion equation is solved numerically with two implicit second-order time integration techniques: the so-called Cranck–Nicholson scheme (CN) that integrates the space terms in time with a trapezoidal rule based on the current and previous time steps; three-point backward finite difference scheme (FD) approximating the time derivative and the space terms determined at the current time-step.

Two different manufactured solutions with Dirichlet boundary conditions not changing in time are selected to mimic two typical cases of unsteady flow: a periodic solution (MSP) and a solution that exhibits a transient behavior before it reaches a steady-state (MSS). Numerical solutions were obtained for seven equally spaced grids using a wide range of time steps. For each space and time discretization, calculations were performed with iterative convergence criteria ranging from machine accuracy (εit=1014 for the maximum residual) to a maximum residual at each time-step of εit=103, which is the default value of several CFD codes. With the knowledge of the exact solution, it is possible to obtain the following results for a quantity of interest ϕ obtained with a given time and space discretization:

  • The numerical error enum from the difference between the solution obtained with εit as iterative convergence criterion ϕit and the exact solution ϕexact:enum(ϕ)=ϕitϕexact.

  • The discretization error ed from the difference between the solution (iteratively) converged to machine accuracy ϕma and the exact solution: ed(ϕ)=ϕmaϕexact.

  • The iterative error eit from the difference between the solutions obtained with εit as iterative convergence criterion and machine accuracy: eit(ϕ)=ϕitϕma.

In the case of the MSP solution, the previous definitions of ed and eit assume that the statistical error, i.e., the difference between solutions of two consecutive periods, has also been reduced to machine accuracy.

In order to address the effect of the initial condition on periodic flow simulations, the MSP is calculated with the initial condition specified from the exact solution and with a wrong initial condition obtained with linear interpolation between the two boundary values. The goal is to check if the four simulations of each test case (CN and FD time integration) reach a numerically periodic solution. Furthermore, we can check whether the solution obtained with the wrong initial condition differs from that obtained with the exact initial condition.

The remainder of this paper is organized in the following way: Sec. 2 presents the techniques used to determine the different contributions to the numerical error in unsteady simulations; the manufactured solutions and the numerical techniques used to determine the solution are described in Sec. 3, and the results obtained in the grid/time refinement studies are presented and discussed in Sec. 4; final remarks in Sec. 5 summarize this work.

Numerical Errors

The numerical error of unsteady flow simulations includes the same contributions as that of steady flow simulations [24]: round-off, iterative and discretization errors. The handling of round-off errors is as for steady flows and will not be addressed in this work; all simulations are performed in double precision (14 digits) to guarantee that its contribution to the numerical error is negligible. However, iterative and discretization errors have some differences for simulations of unsteady compared to steady flow.

Besides, periodic solutions exhibit one more contribution to the numerical error that we will designate by statistical error. It is simple to understand that a simulation that starts from an exact solution will not produce exactly the same result after the simulation of one period due to the existence of numerical errors. This difference between consecutive cycles should tend to zero with the increase of the number of simulated cycles. It is important to emphasize that the statistical convergence is not a consequence of an error in the initial condition. Nonetheless, it is interesting to investigate to what extent an error in the initial condition affects the numerical error of the simulation of a periodic flow.

Statistical Error.

The determination of the statistical error est is straightforward if the period of the solution T is known 
est(ϕ)=ϕ(t+T)ϕ(t)
(1)

where ϕ stands for any functional or local solution quantity. However, for cases where the determination of the period is part of the flow solution, as for example, the flow past a circular cylinder where the period is determined by the vortex shedding in the wake, the assessment of statistical convergence becomes more troublesome and it is outside the scope of this work.

Iterative Error.

Iterative errors in flow solutions are originated by several features of the solution procedure: linearization of the convective terms, segregated solution procedures, deferred corrections in the space discretization, and solution of linear systems of equations with iterative methods. All these features are similar in steady and unsteady flows. However, the iterative errors are more difficult to address in unsteady than in steady simulations because the iterative error at a given time-step propagates to the next time-step.

As for steady flow simulations, the estimation of the iterative error of an unsteady simulation does not require the knowledge of the exact solution [24,10]. For a given grid refinement and time-step, the iterative errors can be estimated by the difference between a solution ϕit calculated with a convergence criterion εit and a solution obtained with the iterative error reduced to machine accuracy ϕma 
eit(ϕit)=ϕitϕma
(2)
Although Eq. (2) is straightforward to apply, it requires a solution converged to machine accuracy, which may not always be affordable. Therefore, a practical procedure is required for the iterative error estimation. As proposed in Ref. [11], the iterative error can be estimated using εit, the convergence criterion used at each time-step, as the independent variable3 in the expression 
eoit(ϕit)=ϕitϕoit=αitexpH(εit)βit
(3)

where ϕoit is the estimate of the solution without iterative error and αit and βit are constants. The function H(εit) must go to minus infinity when εit goes to zero. Two alternatives are tested4: H=ln(εit) and H=1/εitpit, where pit is assumed to be between 0.1 and 3 with Δpit=0.1. If only three levels of εit are used to obtain ϕit it is not possible to choose pit. Therefore, the proposed procedure should be used only with ϕit obtained with at least four levels of εit. In that case, Eq. (3) is solved in the least-squares sense5 to determine ϕoit, αit, and βit and from the two options for H(εit) we choose the one giving the smallest standard deviation in the fit.

Discretization Error.

In unsteady flow simulations, the discretization error has contributions from the space and time discretizations. Therefore, the approximation of the discretization error with a truncated power series expansion leads to 
ed(ϕi)=ϕiϕo=αx(hih1)px+αt(τiτ1)pt
(4)

where ϕo is the estimate of the exact solution, hi is the typical cell size,6τi is the time step6, px is the order of grid convergence, pt is the order of time convergence, and αx and αt are constants. The exact definition of hi is not trivial, and it requires geometrically similar grids [13]. However, in this work, we will not address this problem.

Equation (4) can be used in three different ways to estimate the discretization error.

Simultaneous Grid and Time Refinement.

Assuming that px and pt are the formal orders, known in advance, it is possible to select the time-step so as to simplify Eq. (4) to the standard formulation for steady-state simulations [24] 
ed(ϕi)=ϕiϕo=α(λiλ1)p
(5)
where 
λiλ1=hih1=(τiτ1)px/pt
In this approach, a minimum7 of three data points (ϕ1ϕ3) is required to estimate the discretization error from Eq. (5) leading to a discretization error of ϕ1 given by 
ed(ϕ1)=ϕ1ϕo=ϕ2ϕ1(λ2λ1)p1
(6)

Once ϕo has been determined from Eq. (6), it is possible to estimate the discretization error for any of the three data points.

The observed order of convergence is obtained from the solution of 
|ϕ3ϕ2ϕ2ϕ1|=(λ2λ1)p(λ3λ2)p1(λ2λ1)p1,λ3λ2=λ2λ1p=ln(|ϕ3ϕ2ϕ2ϕ1|)ln(λ2λ1)
(7)

Independent Grid and Time Refinement.

A technique suggested in Ref. [8] is to perform grid refinement for a fixed time-step and combine the estimated error with time refinement for a fixed grid. In this approach, five data points are required derived from computations on three grids for a fixed time-step (ϕ1τm,ϕ2τm, and ϕ3τm) and with three time steps for a fixed grid ϕ1hn,ϕ2hn, and ϕ3hn.

The grid refinement study determines αx, px and 
ϕox=ϕo+αx(hnh1)px
(8)
using equations similar to Eqs. (6) and (7). Analogously, the time refinement exercise determines αt, pt and 
ϕot=ϕo+αt(τmτ1)pt
(9)
Combining the line fits (8) and (9) to form a surface fit, the obvious choice for the estimate of the discretization error of the data point ϕi of grid hn with time-step τm is 
ed(ϕi)=ϕiϕo=(ϕiϕox)+(ϕiϕot)
(10)

Surprisingly, Ref. [8] suggests that ed(ϕi) is to be obtained from the square root of the sum of the squares of ϕox and ϕot, which is not correct; the errors obtained from the two extrapolations must be summed. Furthermore, to use Eq. (10) correctly, the grid refinement study must determine the error of the data point of grid hn and the time refinement study the error obtained with time-step τm. Alternatively, it is possible to determine ϕo and estimate the discretization error of any of the five data points.

Arbitrary Grid and Time Refinement.

The third alternative is to solve directly Eq. (4) using five data points from solutions with arbitrary choices for grid density and time-step, be it that at least two grids and two time steps must be involved. The main drawback of this approach is that a system of five nonlinear equations must be solved to obtain ϕo, αx, αt, px, and pt. However, fixing px and pt and determining a reasonable initial guess for ϕo, αx, αt solving Eq. (4) in the least squares sense leads to a system of linear algebraic equations. Therefore, it is possible to make a robust algorithm to determine the five unknowns, provided that the choices of grid and time refinement do not lead to an ill-conditioned problem.

Model Problem Definition

Manufactured Solutions.

A simple unsteady, one-dimensional, convection–diffusion equation is used to illustrate the determination and estimation of the numerical errors in unsteady flow simulations 
ut+uux1Re2ux2=f(x,t)for0x1andt0
(11)

The reference quantities Uref and L are used to define the dimensionless velocity u, length x, and time t. The constant8 Re = 10 and the function f(x, t) is determined, in accordance with the method of manufactured solutions, so that Eq. (11) holds for a suitably chosen exact solution for u.

The exact solution is chosen to be 
u=cos(π2x)+g(t)sin(πx)
(12)

This definition has intentionally been selected to have boundary conditions independent of time (u(0, t) = 1 and u(1, t) = 0) while the behavior of the solution in time may be varied by the function g(t). The two selected test cases are

  • (1)
    Periodic solution in time (MSP) 
    g1(t)=cos(2πt)
    (13)
  • (2)
    Transient solution that tends to a steady-state (MSS) 
    g2(t)=texp0.2t
    (14)

The function g2(t) exhibits a maximum for t =5 followed by a decay to zero.9

The initial condition is directly available setting t =0 in Eqs. (13) and (14) and the functions f(x, t) for MSP and MSS are given in the Appendix.

Numerical Solution

Space Discretization.

Equation (11) is discretized with a finite volume technique and so it is rewritten for the Nx segments (“volumes”), covering the domain, as 
xixi+1(ut+12uux1Re2ux2f(x,t))dx=0
(15)

where xi and xi+1 are the coordinates of the points (“faces”) at the edges of segment i, fw, and fe. Picard linearization is applied to Eq. (15), which means that uu=u¯u with u¯ taken from the previous iteration. The Nx unknowns uci are stored at the center of each element xci.

Applying the midpoint rule and the divergence theorem, we obtain the following algebraic equation for each segment of the grid 
ucitΔxi+12(u¯feufeu¯fwufw)1Re[(ux)fe(ux)fw]f(xci,t)Δxi=0
(16)

where Δxi is the length of the segment. u¯ is obtained with linear interpolation and central differences are applied to the first-derivatives at the faces of the segment. Fluxes ufe and ufw are obtained with second-order upwind approximation based on uc and (u/x)c with the selected cell center dependent on the sign of u¯fw and u¯fe. Second-order approximations are used to determine first-derivatives at the center of the segments and at the boundaries of the domain, x =0 and x =1. The contributions of (u/x)c,(1/Re)(u/x)fw at x =0 and (1/Re)(u/x)fe at x =1 are dealt with explicitly (deferred corrections), i.e., they are calculated with values from the previous iteration and added to the right-hand side of the system of linear equations. The formal order of grid convergence of the space discretization is px = 2.

Time Integration.

Two alternatives with second-order accuracy (pt = 2) have been tested for the time integration. In the first approach, Eq. (16) is integrated in time as 
tj1tjucitΔxidt=tj1tjS(t)dt
(17)
where tj1 and tj designate the previous and the current time instants and 
S(t)=f(xci,t)Δxi+1Re[(ux)fe(ux)fw]12(u¯feufeu¯fwufw)
Approximating the right-hand side of Eq. (17) with a trapezoidal rule and using Δt to designate the time-step (Δt=tjtj1) leads to 
(uci(tj)uci(tj1))Δxi=12(S(tj)+S(tj1))Δt
(18)
which corresponds to the so-called Crank–Nicolson scheme (CN).
The second alternative (FD) calculates Eq. (16) at time tj and approximates the time derivative with a three-point backward finite difference scheme that for a constant time-step leads to 
(1.5uci(tj)2uci(tj1)+0.5uci(tj2))Δxi=S(tj)Δt
(19)

It should be pointed out that an error in the initial condition will have different consequences for the two time integration techniques. The FD scheme requires only the solution of the previous time steps at the center of the segments, whereas the CN scheme also uses S(tj1) including the space discretization at t =0.

Solution Procedure.

The two alternatives lead to the solution of a nonlinear problem at each time-step with each nonlinear iteration corresponding to a tridiagonal (Aw, Ac, Ae) system of linear algebraic equations that has a straightforward direct solution using the Thomas algorithm 
Awiuci1+Aciuci+Aeiuci+1=RHSiAwiΔuci1+AciΔuci+AeiΔuci+1=Resi
(20)
with the residual Resi obtained from 
Resi=RHSiAwiu¯ci1Aciu¯ciAeiu¯ci+1
(21)

The iterative error is due to three features of the solution procedure: Picard linearization (u¯); second-order corrections of the upwind scheme as deferred corrections; and diffusion fluxes at the boundaries of the domain added to the right-hand side of the system. To enhance the robustness of the solution process, implicit and explicit under-relaxation techniques are applied to the solution of the system (20), i.e., (after determining Resi) Aci is divided by the implicit under-relaxation coefficient γimp, and the solution is updated as uci=u¯ci+γexpΔuci.

Examples of Error Estimation in Unsteady Flow Simulations

The two selected manufactured solutions were numerically solved (formal orders of convergence are px = pt = 2) in sets of seven, geometrically similar, equally spaced grids where the definition of the grid and time-step refinement ratio is straightforward 
hih1=Nx1Nxiandτiτ1=ΔtiΔt1
(22)

Nx1 = 640 and Δt1 = 9.765625 × 10−4, whereas the coarsest grid has Nx1 = 10 and the largest time-step is Δt = 0.5. For each selected grid and time-step, we tested nine different levels of the iterative convergence criterion based on the maximum value of the normalized Resi, (max(Resi/Aci)<εit). Eight of the selected values of εit range from 10−3 to 10−10, whereas the reference solutions to determine the iterative errors were obtained with εit = 10−14. All simulations were performed with γimp = 0.5 and γexp = 0.5 with which iterative convergence at all time steps is guaranteed. For practical reasons, these values were chosen somewhat conservatively; for many of the cases tested a better iterative convergence rate might have been achievable.

Several norms and flow quantities were evaluated to illustrate the behavior of the different kinds of numerical error:

  • L2 and L norms of the numerical error at each time step, L2 (enum) and L (enum);

  • time history of the solutions at x=0.25(u1),x=0.5(u2) and x=0.75(u3) and of the functional quantity Ix=01u(x,t)dx.

All simulations of the MSP case (period T =1) include 1000 cycles, while 250 time units (0<t250) were simulated for the MST case. Grid and time refinement strategies included four cases at constant Courant number (Comax=8.8,Comax=4.4,Comax=2.2,Comax=1.1) and a few extra combinations with time refinement for fixed grids or vice versa. The grid and time-step ranges used in the four studies with simultaneous refinement are presented in Table 1.

All calculations are performed in double precision (14 digits) and one comparison between double and quadruple precision was performed for the MSP using εit=1014 and Comax = 1.1. The maximum differences between the four quantities of interest of the two simulations after 200 periods of simulation are smaller than 2.5 × 10−13 for all values of λi.

Statistical Error.

The test case MSP was deliberately chosen to allow the study of statistical convergence. For each simulation performed, we have determined the number of cycles Ncycles required to reach |ϕ(t)ϕ(t+T)|<1010, where ϕ stands for any of the monitored quantities. The results are presented in Fig. 1 for simultaneous grid and time refinement using CN and FD time integration with the exact and wrong initial conditions. It should be mentioned that no limit was imposed on the number of iterations for solving the nonlinear problem at each time-step, which means that all simulations satisfied the selected εit at all time steps. The results obtained with Comax = 4.4 and Comax = 2.2 are similar to the simulations performed with Comax = 1.1. Therefore, the discussion will be restricted to the plots corresponding to Comax = 1.1 and Comax = 8.8.

The number of cases that do not satisfy the selected statistical convergence criterion is small for all four sets of simulations performed. Simulations not reaching statistical convergence have typically small values of λ (finest grids and smallest time steps) and large values of εit for both the integration techniques (CN and FD). The number of cycles required to satisfy |ϕ(t)ϕ(t+T)|<1010 is similar for the simulations with the smallest Courant numbers tested, leading to 12 cycles for solutions started from the exact initial condition and approximately five cycles more for simulations started with the wrong initial conditions. On the other hand, the simulations performed with Comax = 8.8 required a much larger number of cycles for the CN option when an error is introduced in the initial condition. It must be emphasized, though, that both initial conditions lead to the same periodic solution. Figure 2 presents L[enum(u)(t=999.5)] for the last calculated cycle as a function of λ. The plots contain data obtained with five of the values of εit tested using the CN and FD time integration and the exact and wrong initial conditions.

After statistical convergence, the results obtained with the exact and wrong initial conditions are indistinguishable for both time integration techniques (the open symbols are covered by the corresponding solid ones). The results obtained with ε = 10−14 show the expected second-order convergence of the two time/space integration techniques for the solutions obtained with the two Courant numbers. The results obtained for the largest value of Comax show that the CN option attains the asymptotic regime for smaller values of λi than the FD technique.

Iterative Error.

That the iterative error, governed by εit, can have a significant effect on the numerical error is already illustrated in Fig. 2. For the smaller value of λ, only ε = 10–14 is sufficient to obtain a numerical error free of contamination by the iterative error. Furthermore, εit = 10−3 leads to a considerable increase of the numerical error for all values of λ except the largest λ = 64. The results of Fig. 2 are confirmed by those depicted in Fig. 3, displaying the L norm of the numerical error at t = 5 (enum(u)t=5) where g2(t) attains its maximum value for MST.

There are two main differences compared to the results obtained for the periodic solution MSP: the error constant α for εit=1014 is identical for the CN and FD techniques, whereas different values were obtained for MSP; there is a much smaller influence of the Courant number on the error level than in the MSP case. For both cases the influence of the time discretization is smaller than that of the space discretization.

The present manufactured solutions permit a good estimate of the iterative error eit using the solution obtained with εit=1014 as the reference. Figure 4 presents eit for x=0.5,t=999.5 of the MSP case and x=0.5,t=5 for the MST solution. The results are presented as a function of εit for the seven levels of λ tested. As for the previous case, the plots correspond to the lowest and highest Courant numbers tested.

The results plotted in Fig. 4 exhibit the same trends for the two time integration techniques. There is a small influence of the Courant number for both cases. For a given convergence criterion εit, there is a significant increase of the iterative error with grid/time refinement. With the exception of the two largest values of λ, eit is significantly larger than εit, the difference reaching nearly 4 orders of magnitude for λ = 1. Therefore, estimating the iterative errors in unsteady flow simulations is essential to assess the numerical accuracy of the simulations.

In practice, it is often not possible to determine a reference solution converged to machine accuracy. To illustrate the effort involved in such simulations, Fig. 5 presents the average number of iterations Niter required to satisfy εit=1014 at each time-step of the last simulated cycle of the MSP. For comparison purposes, the values of Niter obtained for εit=106,εit=108, and εit=1010 are also plotted in the same figure. There is a substantial increase of Niter when λ decreases and/or the Courant number increases. For most of the values of λ tested, there is an increase of roughly 1 order of magnitude in Niter when εit changes from 10−6 to 10−14. Although the values of γimp and γexp could be tuned to improve the iterative convergence rate, more than 104 iterations are required to satisfy εit=1014 for λ = 1 and Comax = 8.8 (Niter ≃ 2 × 103 for λ = 1 and Comax = 1.1). These results suggest that the need for an iterative error estimator for practical simulations does not rely on a reference solution converged to machine accuracy.

The error estimates (eoit) based on Eq. (3), which can be justified by the almost linear behavior illustrated in Fig. 6, were tested for u1, u2, u3, and Ix in all the time steps of the last simulated cycle of MSP and for 0<t20 for MST. eoit was determined for λ16 using the data of four levels of εit, which means that eoit was estimated for solutions corresponding to 1010εit106. The assessment of the error estimates performance is based on the ratio Rit=eoit/eit. Ffail stands for the percentage of cases that exhibit Rit < 1, and Fok designates the percentage of cases that satisfy 1Rit5. Bearing in mind that eoit is an error estimate [2,3], we should expect values of Ffail that can reach 50%. However, the results presented in Fig. 6 exhibit much smaller values than 50% for most of the conditions tested. Furthermore, the values of Fok are mostly higher than 50% and Ffail+Fok is close to 100%, which means that eoit is too conservative for a small percentage of cases. Furthermore, over-conservative estimates are mainly obtained for εit=106 which is probably due to the fact that those estimates refer also to the data obtained with εit=103, being the only set of simulations that does not give a good fit to a straight line in Fig. 4. Therefore, Eq. (3) turns out to be a viable alternative to estimate the iterative errors without requiring a solution iteratively converged to machine accuracy.

Discretization Error.

The three techniques described in Sec. 2.3 to estimate the discretization error were applied to obtain the estimate of the exact solution ϕo of u1, u2, u3, and Ix at t =999.5 for MSP and t = 5 for MST. We will designate the simultaneous grid and time refinement by SR, the independent grid and time refinement by IR and the arbitrary grid and time refinement by AR. All the discretization error estimates were made with the simulations iteratively converged to machine accuracy (εit=1014) and for the last simulated cycle for the case of the periodic MSP. These conditions make sure that the numerical errors are dominated by the discretization error.

Four groups of grids (D1–D4) were defined, as detailed in Table 2, to estimate the discretization error of the quantities of interest. Each group has nine data points. Consequently, the number of estimations of ϕo that can be made is: one for SR; 9 for IR and 126 for AR. Unfortunately, not all the 126 combinations of five data points allowed the solution of the nonlinear system of equations to determine ϕo from Eq. (4). Nonetheless, 81 solutions were obtained for D1–D3 and 70 for D4 (coarsest grids and largest time steps). In order to assess the performance of the error estimators, we have determined the ratio between the estimated discretization error and the “exact” discretization error Rd=(ϕiϕo)/(ϕiϕexact) for the smallest value hi/h1 and τi/τ1 of each group. The minimum and maximum values of Rd for the selected quantities of interest are presented in Fig. 7 together with the minimum and maximum values of the observed orders of grid and time convergence p. The results obtained with the two time integration techniques are similar and so only the CN data are presented in Fig. 7.

There is a significant influence of the grid/time refinement level on the accuracy of the estimated discretization errors. For D1 and D2, representing the finest grids and smallest time steps, the results are in excellent agreement with the exact error for the three techniques tested for the MSP with values of p very close to 2. On the other hand, the MST data show a larger spread of estimated values of p that is reflected in the range of Rd obtained. The observed orders of grid and time refinement deviate significantly from 2 for the D3 and D4 groups, which lead to a wider range of values of Rd.

As an illustration of the estimations of the discretization error in unsteady flows, Fig. 8 presents the fits performed with Eqs. (4) (IR,AR) and (5) (SR) for u1 at t =999.5 of MSP and for u2 at t =5 of MST. At this level of grid/time refinement, the different techniques lead to very similar results. However, for the SR technique, Eq. (5), it is not possible to identify the relation between the time and space contributions to the discretization. On the other hand, using Eq. (4) it is possible to show that for the MSP case the time contribution to the discretization error is larger than the space contribution, whereas the opposite trend is obtained for the MST case.

Interaction Between Iterative and Discretization Errors.

The results presented in Figs. 7 and 8 would be significantly different if the discretization error had been estimated with the data of the simulations performed with εit=108, which means that the estimation of the discretization error ed is not independent on the iterative error eit. Therefore, as in Ref. [10], we have used the same data to evaluate the iterative error estimates (presented in Fig. 6) to check if enumeit2+ed2 (RSS). The number of cases that exhibits a numerical error larger than RSS (RSSfail) is depicted in Fig. 9 for the simulations performed with 1010εit106 as a function of the grid/time refinement ratio λ for Comax = 1.1 and Comax = 8.8. The four quantities of interest (u1, u2, u3, and Ix) are checked at all times steps of the last cycle of MSP and between 0 < t <20 for MST.

The results presented in Fig. 9 confirm the conclusion reported in Ref. [10] that discretization and iterative errors are not independent. Hence, RSS is not a conservative estimate of the numerical error. On the other hand, summing the two contributions gives always a proper estimate of the numerical error.

Summary and Conclusions

In this paper, we have dealt with numerical errors in simulations of unsteady flows (as opposed to the simpler case of steady flows). The three important contributions to these numerical errors are statistical, iterative, and discretization errors.

Two manufactured solutions of an unsteady convection–diffusion equation in one space dimension were constructed to illustrate the behavior of the different components of the numerical error. One solution is periodic and the other represents a transient that tends to a steady-state. The obvious advantage of a manufactured solution is that the exact solution is known. Moreover, the chosen flow problems are simple enough to allow a numerical solution converged to machine accuracy. It is then relatively easy to evaluate error components in numerical solutions obtained with less demanding convergence criteria and to monitor the effects of grid and time-step refinement on the error level.

Being aware that for real flow problems in engineering practice a numerical solution converged to machine accuracy is not readily attainable, we also devised a procedure for estimating the iterative error based on a set of numerical solutions with a suitable set of four levels of the iterative convergence criterion. In addition, we have offered three options for estimating the discretization error, making use of a proper set of computations on a sequence of grids and with different time steps.

By analysis of the full set of numerical simulations carried out in this study, we have drawn the following conclusions:

  • The statistical error can be reduced to machine accuracy for grid/time refinement levels and iterative convergence criteria that leads to iterative and discretization errors many orders of magnitude above machine accuracy. However, statistical convergence tends to deteriorate with the increase of the iterative convergence level and grid/time refinement.

  • For an iterative convergence criterion based on the maximum normalized residual at each time step εit, the level required to obtain a negligible contribution of the iterative error to the numerical error is significantly more demanding than what is usually reported in the open literature.

  • A procedure for estimating the iterative error in case a numerical solution converged to machine accuracy is lacking showed promising results for the two selected manufactured solutions.

  • For these well controlled test cases, the three techniques tested for the estimation of discretization errors lead to similar results if the grids and time steps are sufficiently refined, i.e., if the observed orders of grid and time refinement are close to the formal orders of discretization. However, the three methods have different requirements and different information contents:

  • Combining the grid and time refinement leads to the smallest number of data points required (3). The estimation of the exact solution ϕo requires the solution of one nonlinear equation in the general case or the availability of explicit analytical formulas when the refinement ratio is constant. However, there is no information whether it is space or time discretization that gives the largest contribution to the discretization error.

  • Performing grid refinement with constant time step τm next to time refinement for a fixed grid with typical cell size hn provides not only ϕo but also the separate contributions of space and time discretization. If the errors of the two refinement studies are added, the result is equal to the discretization error estimate for grid hn with time step τm.

  • Choosing five data points including at least two grids and two time steps also provides information about the space and time contributions to the discretization error. However, the determination of ϕo requires the solution of a system of five nonlinear equations which may turn out to be unfeasible for some combinations of grids and time steps.

  • Iterative and discretization errors are not independent contributions to the numerical error. Therefore, the square root of the sum of its squares does not provide a conservative estimate of the numerical error; instead, the sum of the two contributions is to be used.

Of course, the quality of many of these error estimates depends on the conditions which may be hard to meet under the circumstances of engineering application of CFD, as for example: the generation of geometrically similar grids; the ability to iteratively converge the solution at each time-step to the desirable tolerance; the ability to reach grid/time refinement levels that lead to orders of grid/time refinement close to the formal order of the schemes; the absence of limiters and/or switches in the discretization scheme. Therefore, numerical uncertainty estimation for numerical flow simulations in industrial practice presents even more challenges than those discussed in this work.

Appendix: Manufactured Solutions

Solutions: 
u(x,t)=cos(π2x)+g(t)sin(πx)
(A1)
 
MSPg(t)=cos(2πt)MSTg(t)=texp0.2t
(A2)
Source terms: 
f(x,t)=ut+uux1Re2ux2
(A3)
where 
MSPut=2πsin(2πt)sin(πx)MSTut=(10.2t)exp0.2tsin(πx)
(A4)
 
ux=π2sin(π2x)+g(t)πcos(πx)
(A5)
and 
2ux2=π24cos(π2x)g(t)π2sin(πx)
(A6)
2

We are assuming that double-precision (14 digits) is sufficient to obtain a negligible contribution of the round-off error to the numerical error.

3

The convergence criterion must be related to residuals or changes between nonlinear iterations. Number of iterations should not be used, because iterative convergence may stagnate.

4

The first option corresponds to a power series expansion with a single term. Only restriction applied is β>0 for all cases.

5

Solution procedure is identical to that described in Ref. [12] for the estimation of discretization errors in steady flows.

6

The subscript 1 is used to identify the finest grid and the smallest time-step employed.

7

If the equation is solved in the least squares sense as proposed in Ref. [13], more data points can be used. However, in this work, we will restrict ourselves to the minimum number of data points required using ϕ1 for the most refined grid/time-step.

8

If Eq. (11) represents the momentum balance for an incompressible fluid, Re is the Reynolds number.

9

For t =200, g2(t) is already smaller than 10−15.

References

References
1.
EU
,
2019
, “
The European Project TILDA
,” accessed June 20, 2019, www.hifiled-symposium.eu
2.
Roache
,
P. J.
,
1998
,
Verification and Validation in Computational Science and Engineering
,
Hermosa Publishers
,
Albuquerque, NM
.
3.
Roache
,
P. J.
,
2009
,
Fundamentals of Verification and Validation
,
Hermosa Publishers
,
Albuquerque, NM
.
4.
Oberkampf
,
W. L.
, and
Roy
,
C. J.
,
2010
,
Verification and Validation in Scientific Computing
,
Cambridge University Press
,
Cambridge, UK
.
5.
ASME
,
2007
, “
Guide for the Verification and Validation in Computational Solid Mechanics
,” American Society of Mechanical Engineers, New York, Standard No. V V10.
6.
ASME
,
2009
, “
Standard for Verification and Validation in Computational Fluid Dynamics and Heat Transfer
,” American Society of Mechanical Engineers, New York, Standard No. V V20.
7.
AIAA
,
1998
, “
Guide for the Verification and Validation of Computational Fluid Dynamics Simulations
,” American Institute of Aeronautics and Astronautics, Reston, VA, Standard No. AIAA-G077-1998.
8.
ITTC
,
2017
, “ITTC Quality Manual, Uncertainty Analysis in CFD, 22nd ITTC [C], 1999,” ITTC.
9.
Roache
,
P. J.
,
2002
, “
Code Verification by the Method of the Manufactured Solutions
,”
ASME J. Fluids Eng.
,
124
(
1
), pp.
4
10
.
10.
Eça
,
L.
, and
Hoekstra
,
M.
,
2009
, “
Evaluation of Numerical Error Estimation Based on Grid Refinement Studies With the Method of the Manufactured Solutions
,”
Comput. Fluids
,
38
(
8
), pp.
1580
1591
.
11.
Eça
,
L.
,
Vaz
,
G.
, and
Hoekstra
,
M.
,
2018
, “
On the Role of Iterative Errors in Unsteady Flow Simulations
,”
Numerical Towing Tank Symposium
(
NuTTS
),
Cortona, Italy
,
Sept. 30–Oct. 2
.https://www.researchgate.net/publication/328064499_On_the_Role_of_Iterative_Errors_in_Unsteady_Flow_Simulations
12.
Eça
,
L.
, and
Hoekstra
,
M.
,
2014
, “
A Procedure for the Estimation of the Numerical Uncertainty of CFD Calculations Based on Grid Refinement Studies
,”
J. Comput. Phys.
,
262
, pp.
104
130
.
13.
Eça
,
L.
,
Klaij
,
C. M.
,
Vaz
,
G.
,
Hoekstra
,
M.
, and
Pereira
,
F. S.
,
2016
, “
On Code Verification of RANS Solvers
,”
J. Comput. Phys.
,
310
(
C
), pp.
418
439
.