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 [2–4], 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. [2–4] 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 errors^{2} 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 ($\epsilon it=10\u221214$ for the maximum residual) to a maximum residual at each time-step of $\epsilon it=10\u22123$, 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

*e*_{num}from the difference between the solution obtained with*ε*_{it}as iterative convergence criterion*ϕ*_{it}and the exact solution*ϕ*_{exact:}$enum(\varphi )=\varphi it\u2212\varphi exact$.The discretization error

*e*from the difference between the solution (iteratively) converged to machine accuracy_{d}*ϕ*_{ma}and the exact solution: $ed(\varphi )=\varphi ma\u2212\varphi exact$.The iterative error

*e*_{it}from the difference between the solutions obtained with*ε*_{it}as iterative convergence criterion and machine accuracy: $eit(\varphi )=\varphi it\u2212\varphi ma$.

In the case of the MSP solution, the previous definitions of *e _{d}* and

*e*

_{it}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 [2–4]: 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.

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.

*ϕ*

_{it}calculated with a convergence criterion

*ε*

_{it}and a solution obtained with the iterative error reduced to machine accuracy

*ϕ*

_{ma}

*ε*

_{it}, the convergence criterion used at each time-step, as the independent variable

^{3}in the expression

where $\varphi oit$ is the estimate of the solution without iterative error and *α*_{it} and *β*_{it} are constants. The function $H(\epsilon it)$ must go to minus infinity when *ε*_{it} goes to zero. Two alternatives are tested^{4}: $H=ln(\epsilon it)$ and $H=\u22121/\epsilon itpit$, where *p*_{it} is assumed to be between 0.1 and 3 with $\Delta pit=0.1$. If only three levels of *ε*_{it} are used to obtain *ϕ*_{it} it is not possible to choose *p*_{it}. 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 sense^{5} to determine $\varphi oit$, *α*_{it}, and *β*_{it} and from the two options for $H(\epsilon it)$ we choose the one giving the smallest standard deviation in the fit.

### Discretization Error.

where *ϕ _{o}* is the estimate of the exact solution,

*h*is the typical cell size,

_{i}^{6}

*τ*is the time step

_{i}^{6},

*p*is the order of grid convergence,

_{x}*p*is the order of time convergence, and

_{t}*α*and

_{x}*α*are constants. The exact definition of

_{t}*h*is not trivial, and it requires geometrically similar grids [13]. However, in this work, we will not address this problem.

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

#### Simultaneous Grid and Time Refinement.

^{7}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

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

#### 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 ($\varphi 1\tau m,\u2009\varphi 2\tau m$, and $\varphi 3\tau m$) and with three time steps for a fixed grid $\varphi 1hn,\u2009\varphi 2hn$, and $\varphi 3hn$.

Surprisingly, Ref. [8] suggests that $ed(\varphi i)$ is to be obtained from the square root of the sum of the squares of $\varphi ox$ and $\varphi 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 *h _{n}* and the time refinement study the error obtained with time-step

*τ*. Alternatively, it is possible to determine

_{m}*ϕ*and estimate the discretization error of any of the five data points.

_{o}#### 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}*p*, and

_{x}*p*. However, fixing

_{t}*p*and

_{x}*p*and determining a reasonable initial guess for

_{t}*ϕ*,

_{o}*α*,

_{x}*α*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.

_{t}## Model Problem Definition

### Manufactured Solutions.

The reference quantities *U*_{ref} and *L* are used to define the dimensionless velocity *u*, length *x*, and time *t*. The constant^{8} 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*.

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

The function *g*_{2}(*t*) exhibits a maximum for *t *=* *5 followed by a decay to zero.^{9}

### Numerical Solution

#### Space Discretization.

*Nx*segments (“volumes”), covering the domain, as

where *x _{i}* and $xi+1$ are the coordinates of the points (“faces”) at the edges of segment

*i*,

*f*, and

_{w}*f*. Picard linearization is applied to Eq. (15), which means that $uu=u\xafu$ with $u\xaf$ taken from the previous iteration. The

_{e}*Nx*unknowns $uci$ are stored at the center of each element $xci$.

where $\Delta xi$ is the length of the segment. $u\xaf$ is obtained with linear interpolation and central differences are applied to the first-derivatives at the faces of the segment. Fluxes *u _{fe}* and

*u*are obtained with second-order upwind approximation based on

_{fw}*u*and $(\u2202u/\u2202x)c$ with the selected cell center dependent on the sign of $u\xaffw$ and $u\xaffe$. Second-order approximations are used to determine first-derivatives at the center of the segments and at the boundaries of the domain,

_{c}*x*=

*0 and*

*x*=

*1. The contributions of $(\u2202u/\u2202x)c,\u2009(1/Re)(\u2202u/\u2202x)fw$ at*

*x*=

*0 and $(1/Re)(\u2202u/\u2202x)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*

*p*= 2.

_{x}#### Time Integration.

*p*= 2) have been tested for the time integration. In the first approach, Eq. (16) is integrated in time as

_{t}*t*to designate the time-step $(\Delta t=tj\u2212tj\u22121)$ leads to

*t*and approximates the time derivative with a three-point backward finite difference scheme that for a constant time-step leads to

_{j}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(tj\u22121)$ including the space discretization at *t *=* *0.

#### Solution Procedure.

*Aw*,

*Ac*,

*Ae*) system of linear algebraic equations that has a straightforward direct solution using the Thomas algorithm

*obtained from*

_{i}The iterative error is due to three features of the solution procedure: Picard linearization ($u\xaf$); 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 Res* _{i}*)

*Ac*is divided by the implicit under-relaxation coefficient

_{i}*γ*

_{imp}, and the solution is updated as $uci=u\xafci+\gamma \u2009exp\u2009\Delta uci$.

## Examples of Error Estimation in Unsteady Flow Simulations

*p*=

_{x}*p*= 2) in sets of seven, geometrically similar, equally spaced grids where the definition of the grid and time-step refinement ratio is straightforward

_{t}*Nx*_{1} = 640 and Δ*t*_{1} = 9.765625 × 10^{−4}, whereas the coarsest grid has *Nx*_{1} = 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 Res* _{i}*, ($max(Resi/Aci)<\epsilon it$). Eight of the selected values of

*ε*range from 10

_{it}^{−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:

*L*_{2}and*L*norms of the numerical error at each time step,_{∞}*L*_{2}(*e*_{num}) and*L*(_{∞}*e*_{num});time history of the solutions at $x=0.25(u1),\u2009x=0.5(u2)$ and $x=0.75(u3)$ and of the functional quantity $Ix=\u222b01u(x,t)dx$.

All simulations of the MSP case (period *T *=* *1) include 1000 cycles, while 250 time units ($0<t\u2264250$) were simulated for the MST case. Grid and time refinement strategies included four cases at constant Courant number ($Comax=8.8,\u2009Comax=4.4,\u2009Comax=2.2,\u2009Comax=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 $\epsilon it=10\u221214$ and Co_{max} = 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 *N*_{cycles} required to reach $|\varphi (t)\u2212\varphi (t+T)|<10\u221210$, 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 Co_{max} = 4.4 and Co_{max} = 2.2 are similar to the simulations performed with Co_{max} = 1.1. Therefore, the discussion will be restricted to the plots corresponding to Co_{max} = 1.1 and Co_{max} = 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 $|\varphi (t)\u2212\varphi (t+T)|<10\u221210$ 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 Co_{max} = 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\u221e[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 Co_{max} 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

*g*

_{2}(

*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 $\epsilon it=10\u221214$ 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 *e*_{it} using the solution obtained with $\epsilon it=10\u221214$ as the reference. Figure 4 presents *e*_{it} 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 *λ*, *e*_{it} 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 *N*_{iter} required to satisfy $\epsilon it=10\u221214$ at each time-step of the last simulated cycle of the MSP. For comparison purposes, the values of *N*_{iter} obtained for $\epsilon it=10\u22126,\u2009\epsilon it=10\u22128$, and $\epsilon it=10\u221210$ are also plotted in the same figure. There is a substantial increase of *N*_{iter} 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 *N*_{iter} 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 10^{4} iterations are required to satisfy $\epsilon it=10\u221214$ for *λ* = 1 and Co_{max} = 8.8 (*N*_{iter} ≃ 2 × 10^{3} for *λ* = 1 and Co_{max} = 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 *u*_{1}, *u*_{2}, *u*_{3}, and *I _{x}* in all the time steps of the last simulated cycle of MSP and for $0<t\u226420$ for MST. $eoit$ was determined for $\lambda \u226416$ using the data of four levels of

*ε*

_{it}, which means that $eoit$ was estimated for solutions corresponding to $10\u221210\u2264\epsilon it\u226410\u22126$. The assessment of the error estimates performance is based on the ratio $Rit=eoit/eit$.

*F*

_{fail}stands for the percentage of cases that exhibit

*R*

_{it}< 1, and

*F*

_{ok}designates the percentage of cases that satisfy $1\u2264Rit\u22645$. Bearing in mind that $eoit$ is an error estimate [2,3], we should expect values of

*F*

_{fail}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

*F*

_{ok}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 $\epsilon it=10\u22126$ which is probably due to the fact that those estimates refer also to the data obtained with $\epsilon it=10\u22123$, 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

*u*

_{1},

*u*

_{2},

*u*

_{3}, and

*I*at

_{x}*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 ($\epsilon it=10\u221214$) 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

*ϕ*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=(\varphi i\u2212\varphi o)/(\varphi i\u2212\varphi exact)$ for the smallest value $hi/h1$ and $\tau i/\tau 1$ of each group. The minimum and maximum values of

_{o}*R*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

_{d}*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 *R _{d}* 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

*R*.

_{d}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 *u*_{1} at *t *=* *999.5 of MSP and for *u*_{2} 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 $\epsilon it=10\u22128$, which means that the estimation of the discretization error *e _{d}* is not independent on the iterative error

*e*

_{it}. Therefore, as in Ref. [10], we have used the same data to evaluate the iterative error estimates (presented in Fig. 6) to check if $enum\u2264eit2+ed2$ (RSS). The number of cases that exhibits a numerical error larger than RSS (RSS

_{fail}) is depicted in Fig. 9 for the simulations performed with $10\u221210\u2264\epsilon it\u226410\u22126$ as a function of the grid/time refinement ratio

*λ*for Co

_{max}= 1.1 and Co

_{max}= 8.8. The four quantities of interest (

*u*

_{1},

*u*

_{2},

*u*

_{3,}and

*I*) are checked at all times steps of the last cycle of MSP and between 0 <

_{x}*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

*ϕ*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._{o}Performing grid refinement with constant time step

*τ*next to time refinement for a fixed grid with typical cell size_{m}*h*provides not only_{n}*ϕ*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_{o}*h*with time step_{n}*τ*._{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

*ϕ*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._{o}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

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

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.

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

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

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

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.

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

For *t *=* *200, *g*_{2}(*t*) is already smaller than 10^{−15}.