Verification of Numerical Solutions of the Advection-Diffusion and Burgers Equations

Where E is caused by four sources of error [1]: truncation, iteration, round-off andprogramming. When the other sources do not exist or are very small compared to thetruncation error, E can also be called by discretization error. In practical situations, a numerical solution is obtained because the analytical solution isunknown. Hence, the true value of the numerical error is also unknown and must therefore beestimated. There are two estimators for the discretization error that are widely used with themethods of finite difference and finite volume, both of them based on Richardsonextrapolation. One of them is Richardson estimator, which is given by [2]


Introduction
Accurate and reliable numerical solutions depend on the estimation of the numerical error (E), which can be defined as the difference between the exact analytical solution (Φ) of a variable of interest and its numerical solution (φ), i.e., ( ) -ϕ ϕ = Φ E (1) Where E is caused by four sources of error [1]: truncation, iteration, round-off andprogramming. When the other sources do not exist or are very small compared to thetruncation error, E can also be called by discretization error. In practical situations, a numerical solution is obtained because the analytical solution isunknown. Hence, the true value of the numerical error is also unknown and must therefore beestimated. There are two estimators for the discretization error that are widely used with themethods of finite difference and finite volume, both of them based on Richardsonextrapolation. One of them is Richardson estimator, which is given by [2] Ri r (2) where 1 2 φ φ and =numerical solutions obtained, respectively, with fine (h 1 ) and coarse (h 2 ) grids, h=grid spacing or distance between every two consecutive nodes of the grid (in this work, h=1/N), N+1=number of nodes in the grid, r=h 2 /h 1 =grid refinement ratio, p L =asymptotic order [3] of the error predicted for each variable of interest based on a priori analyses using a Taylor series [4].
The accurate estimation (Ri/E→1) of E by means of Eq. (2) is only possible if the correct value of p L is used, and h is sufficiently small (h → 0). One way of checking in practice, i.e., through numerical experiments, if the valuededuced a priori for p L is correct is to use the concept of effective order (p E ) of the true error, defined by [5] as: As can be seen in Eq. (3), the effective order (p E ) is a function of the true error of the variable of interest. Thus, for problems whose analytical solution is known, it can be used to verify a posteriori if, as h → 0, one obtains the asymptotic order (p L ) of the discretization error, whichis a theoretical result obtained a priori. In problems whose analytical solution is unknown, the concept of apparent order [1] can be used to obtain a posteriori the p L of each φ.
In this work, two one-dimensional (1D) problems are solved by the finite difference method with uniform grids. First and second-order accurate numerical approximations in space, and a mixture of the two (hybrid), are used. The objectives of this work are: (1) based on numerical experiments, verify [3] the true value of the discretization error (E) as a function of the size (h) of the grid for four variables of interest (φ) in each problem; (2) deduce the value of the asymptotic order (p L ) of the discretization error of each φ, be it primary or secondary; (3) also based onnumerical experiments, verify if the value of the effective order (p E ) of the discretization error tends toward p L when h → 0; and (4) over the discretization error and its order, show the effect caused by the Peclet number, by the numerical approximation used, by the mix factor of hybrid schemes, and by the source term. Another objective is to clarify inconsistent statements found in the literature. For example, Celik and Zhang [6] state that the asymptotic order of a hybrid scheme is variable. On the other hand, Roache [2] suggests using the lower order between the two pure schemes. Even in onedimensional problem, Roy [7] considers that the error is reduced non-monotonically when atleast two schemes with different asymptotic orders are used.
The importance of this work lies in confirming the correct value of the asymptotic order (p L ) of the discretization error for some numerical approximations that are very common in the method of finite difference. This will allow the Richardson estimator and its variants to be used correctly, since they depend directly on the value of p L , as can be seen in Eq. (2). The correct value of the asymptotic order is also important to use the repeated Richardson extrapolation efficiently [8]. Another contribution is to show the behavior of the error with the deferred correction scheme [9][10][11]. In addition, we intend to clarify the aforementioned open questions in the literature. One-dimensional problems are used here for the following reasons: (1) the possibility of using highly refined grids in one direction up to an order of millions of nodes, which enables one to verify asymptotic behaviors and thus unequivocally prove the value of practical p L ; (2) due to the speed in obtaining each solution, a large number of tests can be performed, facilitating systematic studies of several parameters; and (3) it is presumed that the onedimensional results are applicable to two and three dimensions. The work is divided as follows: sections II and III address two one-dimensional problems in a steady state, the advectiondiffusion and Burgers equations, and section IV reports the conclusion of this work.

Mathematical model
Starting from the thermal energy conservation equation, considering a steady state 1D incompressible fluid flow without heat generation and viscous dissipation, and consideringthat the properties and velocities are constant in a continuous medium, one obtains theadvection-diffusion equation: where P e =the Peclet number, x=coordinate direction, and u=temperature. The boundary conditions are of the Dirichlet type: The variables of interest, i.e., the variables for which the numerical solution is obtainedand its discretization error and effective order are verified, are: (a) Temperature (u) at x=½: the principal variable of the problem, which is obtained fromthe solution of Eq. (4). Its analytical solution is (b) Mean temperature (U): a global variable obtained from the following definition Its analytical solution is Pe e (8) (c) Inclination (I): a local variable obtained from the following definition Its analytical solution is (d) The mean of the l 1 norm of the discretization error (L), defined mathematically by where i represents each of the nodes of the grid. Its analytical solution has a null value.

Numerical model
The numerical model is characterized by the use of the finite difference method [4], uniform grids and the TDMA (Tri Diagonal Matrix Algorithm) method to solve the system ofequations. The programming language used is FORTRAN 95 with quadruple precision. The diffusive term (second-order derivative) of Eq. (4) was approximated with the centraldifferencing scheme (CDS). In the case of the advective term (first-order derivative), three approximations were used: (1) CDS, (2) first-order upwind differencing scheme (UDS), and(3) β scheme, which is a hybrid UDS and CDS scheme by means of the deferred correction of [10,11], i.e., Where Φ=exact value of the advective term, φ=numerical approximation in the currentiteration, φ*=numerical approximation in the previous iteration, and β=mix factor whosevalue varies from 0 (UDS) to 1 (CDS). In the deferred correction scheme, all the terms that involve β are considered source terms, remaining in the independent term of the system of equations. The coefficients of the system of equations are the same as those of the pure UDS scheme. The matrix of coefficients is of the tridiagonal type. For the UDS and CDS schemes, the TDMA solves the system directly. In the case of the β scheme, due to the second term of Eq. (12), the solution is iterative; in this case, the initial estimation is equal to the analyticalsolution. Details about the above described numerical model and those of the next sections aregiven by [4,11]. The variable of mean temperature (U) was obtained through the trapezoidal rule [12,13], andthe variable of inclination (I) was obtained with the UDS-2 scheme, i.e., the second-order UDS [4].

Estimation of the asymptotic order (p L )
Based on the Taylor series and following the procedure of [4], the truncation error (ε) of the discretized differential equation (EDD) at each i node of the grid is By definition, the asymptotic order (p L ) is the lowest exponent of the error, whose termprevails when h → 0. For uniform grids, it is known [3] that the p L of the discretization errorof the unknown of the differential equation, E(u), is equal to the p L of ε (EDD). Therefore, for u, one has For U, obtained by means of the trapezoidal rule, one has [11] p L =2. This is the sameresult as for I, obtained with the UDS-2 scheme [4]. But these values are valid for thesituation in which u i has no error, i.e., using the analytical solution in each node. In the case ofsecondary variables (φ), i.e., variables that depend on the primary variable (u), in other words, the unknown in the differential equation, due to error propagation, the p L of the discretizationerror (E) is given by where Min=minimum value between the two arguments, and ε(φ)=truncation error of a numerical approximation φ. Hence, the result of Eq. (14) is also valid for U, I and L.

Numerical results
The numerical solution of the four variables of interest was obtained with grids of 3, 5, …up to 67,108,865 nodes, which correspond to h=½, ¼, ... down to ≈ 1.49×10 -8 m. In thesolution of Eq. (4), the schemes UDS, CDS and β=0.999 were used for Pe=1 and 10. Acomputer program was implemented in FORTRAN 95 language, version 9.1 of Intel, usingquadruple precision (Real*16). The simulations were performed on a microcomputerequipped with an Intel processor (Xeon Quad Core X5355 2.66 GHz), 16 GB RAM andWindows XP 64 bits operational system. The maximum CPU time was 9 minutes for the β scheme with eight iterations, a sufficient number to reach the machine round-off error. Figure 1 shows the modulus of the discretization error (E) of the four variables as a function of the h grid used. Considering the modulus of E, one can see that for a relativelylarge h, in general |E(β)| ≈ |E(CDS)| < |E(UDS)|. And for h → 0, |E(CDS)| < |E(β)| < |E(UDS)|. Moreover, for the same scheme, in general |E(Pe=1)| < |E(Pe=10)|. Figure 2 shows the effetive order (p E ) of the error of the four variables as a function ofthe h grid used. All the results of pE presented in this work were obtained with Eq. (3) and a refinement ratio (r) equal to 2. Note, in Figure 2, that: 1) In the coarser grids, as expected [5], the values of the effective order (pE) can besignificantly different from those of the asymptotic order (p L ), presenting negative or evenundefined values.
2) For h → 0, p E → p L as predicted by Eq. (14) for the three schemes (UDS, CDS and β) and the four variables of interest, even for β scheme with its value equal to 0.999.   significantly smallerthan that of the UDS, depending on the value of β (Figure 1).
The literature contains the following inconsistent statements concerning the numerical solution of the advection-diffusion equation: (i) Celik and Zhang [6] state that the asymptotic order of a hybrid scheme is variable; (ii) Roache [2] suggests using the lower order between the two pure schemes; and (iii) even in one-dimensional problems, Roy [7] considered that the error is reduced non-monotonically when at least two schemes with different asymptotic orders are used. Therefore, based on the theoretical and numerical results presented in this section, the statement of Roache [2] is correct, that of Celik et al. [6] is incorrect, and that of Roy [7] is correct depending on the variable of interest.

Mathematical model
Starting from the equation of conservation of the momentum, considering a continuous medium of incompressible fluid with constant properties and steady state 1D laminar flow, one obtains the Burgers equation with a source term given by where Re=Reynolds number, x=coordinate direction, and u=velocity. The boundaryconditions are of the Dirichlet type, given by Eq. (5).
The variables of interest, i.e., the variables for which the numerical solution is obtained and its discretization error and effective order are verified, are: 4) For the mean velocity and the inclination, in the case of β scheme, p E ≈ p L (CDS) in thecoarsest grids. Upon reducing h, there is an interval in which p E varies monotonicallydown to h → 0, p E → p L (UDS). 5) Although for β scheme, p E → p L (UDS) for h → 0, its error may be significantly smallerthan that of the UDS, depending on the value of β (Figure 3).
The same comments made in the last paragraph of the section about the advectiondiffusionequation also apply to the results of the problem in this section with respect toBurgers equation.

Conclusion
This work found principally that: 1) For h → 0, in all the cases, p E → p L , as predicted by Eq. (14), and monotonically for asufficiently small h. 2) The p L of a hybrid scheme is equal to the p L of the pure scheme of the lowest order.
3) For hybrid schemes, the value of the error modulus falls within those of the pure schemes,except in very coarse grids. The proximity of the value of the error modulus between thehybrid scheme and the scheme of the highest order depends on the value employed for the β factor.
4) The value of the parameters Pe and β can affect the value of the error significantly.