Abstract

We perform a spectral analysis of the dispersive and dissipative properties of two time-splitting procedures, namely, locally one-dimensional (LOD) Lax-Wendroff and LOD (1, 5) [9] for the numerical solution of the 2D advection-diffusion equation. We solve a 2D numerical experiment described by an advection-diffusion partial differential equation with specified initial and boundary conditions for which the exact solution is known. Some errors are computed, namely, the error rate with respect to the norm, dispersion and dissipation errors. Lastly, an optimization technique is implemented to find the optimal value of temporal step size that minimizes the dispersion error for both schemes when the spatial step is chosen as 0.025, and this is validated by numerical experiments.

1. Introduction

The advection-diffusion equation is a parabolic partial differential equation combining the diffusion and advection (convection) equations, which describes physical phenomena where particles, energy, or other physical quantities are transferred inside a physical system due to two processes: diffusion and advection [1]. The numerical solution of advection-diffusion equation plays an important role in many fields of science and engineering. These include the transport of air and groundwater pollutants, oil reservoir flow [2], heat transfer in draining film, flow through porous media, the dispersion of pollutants in rivers and streams, long range transport of pollutants in the atmosphere, thermal pollution in river systems, and dispersion of dissolved salts in groundwater [3].

The 3D advection-diffusion equation is given by where , , and are the velocity components of advection in the directions of , , and , respectively, and , , and are the coefficients of diffusivity in the -, -, and -directions, respectively.

This study deals with the 2D advection-diffusion equation, where , in the domain , , with specified initial and boundary conditions.

Dehghan [3] proposed two time-splitting procedures for the solution of the two-dimensional transport equation. The time-splitting procedure used is the locally-one dimensional (LOD) in proceeding from one time step to the next step. LOD replaces the complicated multidimensional partial differential equations by a sequence of solutions of simpler one-dimensional partial differential equations. The originality in this work is that we perform a spectral analysis of the dispersion and dissipation properties of the two schemes at some values of the temporal and spatial step sizes. 3D plots of the relative phase error per unit time step (RPE) and the modulus of the amplification factor (AFM) versus phase angles in - and -directions are obtained. We then use optimization strategies to compute the optimal values of the temporal step size for the two schemes when the spatial step size is 0.025. We then validate this result by performing a 2D numerical experiment with specified initial and boundary conditions. Lastly, we have the conclusion and references.

2. Numerical Dispersion and Dissipation

Dissipation reduces the amplitude of sinusoids in a Fourier series. This is caused by the presence of derivatives like and in the modified equation [4]. On the other hand, the amplitude of sinusoids in a Fourier series is increased by antidissipation. Derivatives like and in the modified equation are generally antidissipative. Dispersion affects the speed of sinusoids in a Fourier series causing phase lag or phase lead and is caused due to the presence of odd-order derivatives in the modified equation [4].

The modulus of the amplification factor (AFM) is a measure of the stability of a scheme and it is also used to measure the dissipative characteristics of the scheme. If the modulus of the amplification factor is equal to one, a disturbance neither grows nor damps [5]. If the modulus of the amplification factor is greater than one, then the scheme is unstable [6]; if it is less than one, damping occurs [7]. The partial differential equation given by (2) is dissipative in nature due to the terms and .

The relative phase error (RPE) is a measure of the dispersive characteristics of a scheme. The relative phase error of a scheme approximating the 1D advection-diffusion equation is given by where is the Courant number, is phase angle, is the amplification factor of the numerical scheme approximating the 1D advection-diffusion equation, and and are the real and imaginary parts of , respectively [8].

We extend the work on the relative phase error in [8] for the case of the 2D advection-diffusion equation. The relative phase error for a numerical scheme approximating (2) is obtained on substituting by [9], where , is the time, and are wavenumbers in the directions of and , respectively, and is the dispersion relation. Thus, we get On substituting (4) into (2), we obtain

The exact phase velocity is and we obtain

The amplification factor can be written as , where and are the real and imaginary parts of , respectively. Also, we can express as , where is time step and is exponential growth rate [9]. Thus, we obtain

The numerical phase velocity is calculated as and we get The relative phase error is the ratio of the numerical phase velocity to the exact phase velocity [10]. It is calculated as The phase angles in the directions of and are given by and , respectively. Hence, the relative phase error is where and are the Courant numbers in the directions of and , respectively.

3. Time-splitting Procedures and Numerical Experiments

The domain we consider is . We divide the spatial interval along - and -directions into and nodes, respectively, such that and and also divide the time interval into grid points such that . Then the grid points are defined by For simplification we take . Let ; then the parameters and represent the spatial and temporal grid spacing, respectively. We denote the approximated value of at the grid point by .

3.1. Time-Splitting Procedures

Since one-dimensional schemes are easier to use than two-dimensional schemes, (2) is split into the following two one-dimensional equations: Each of (12) and (13) can be solved over half of a time step to be used for the complete 2D advection-diffusion equation, using the procedures developed for the 1D advection-diffusion equation.

Some work on time-splitting procedures can be found in [3, 11]. In this paper, we refer to [3] on how a 2D advection-diffusion equation is converted into two 1D advection-diffusion equations using the locally one-dimensional (LOD) time-splitting procedure. Solving (12) and (13) in each half time step is equivalent to solving the following equations over a full-time step: Then we can use the schemes used for solving the 1D advection-diffusion equation to solve (14) and (15).

3.2. Numerical Experiments

In [12], different explicit and implicit finite differences schemes are used to solve the 1D advection-diffusion equation Three values of the cell Reynolds number, namely, , and , are used for the numerical experiments. In [3], the time-splitting procedure is used to solve the two-dimensional transport equation in the region , , with and with the following boundary and initial conditions: for which the exact solution is We consider two time-splitting procedures LOD Lax-Wendroff and LOD to solve (2) with and , subject to the boundary conditions (18)–(22) at time .

We consider two values of ; and . Since and , we have . Thus for and 4, we have and , respectively. When , we have and and for we have and .

3.3. Quantification of Errors from Numerical Results

In this subsection, we describe how errors from numerical results can be quantified into dispersion and dissipation by a technique devised by Takacs [13].

The total mean square error in the 1D case [13] is calculated as where is the exact solution, is the numerical solution at a grid point , and is the number of grid points. The total mean square error can be expressed as

The right hand side of (25) can be rewritten as where and denote the variance of and , respectively, and and denote the mean values of and , respectively. Then we have Therefore, where .

The total mean square error can be expressed as where is the coefficient of correlation.

The term measures the dispersion error and the term measures the dissipation error.

We extend the work on quantification of errors in [14, 15] for the 2D case. The total mean square error for the 2D case is calculated as where and are the exact and numerical solutions at a grid point (), respectively, and Since we have

Hence, where . The dissipation error is and the dispersion error is .

The error rate with respect to norm for is calculated as

4. Construction of the LOD Lax-Wendroff Procedure

We use the following approximations in the first time step of the LOD procedure [3]: where is the spatial weighting factor in the direction of . The following approximations are used for the second step of the LOD procedure: where is the spatial weighting factor in the direction of . By using the following relationships the finite difference formula for the first half time step of the LOD Lax-Wendroff procedure is given by At the second half time step the finite difference is given by where We obtain a single expression for in a complete time step as follows: To find the modified equation of the scheme, we first find Taylor’s expansion of each term in (42) about . The Taylor series expansion of is given by

The Taylor series expansions for some grid points about are given as follows: Using (2), we convert the temporal derivatives , , and into spatial derivatives. We thus have Then on substituting the Taylor series expansions of each term of the difference scheme we get the following modified equation for the LOD Lax-Wendroff The scheme is second order accurate in space and the leading error terms are dispersive in nature (presence of odd-order derivatives and ). As the time and spatial increments go to zero, the modified equation (46) reduces to its original equation, that is, (2). Hence LOD Lax-Wendroff is consistent.

We now study the spectral analysis of the dispersive and dissipative properties of the scheme for the case and . To obtain the amplification factor we use the Von Neumann stability analysis by substituting by in (42), where . We thus have The real and imaginary parts of the amplification factor are given by respectively. The modulus of the amplification factor, AFM, is obtained as We find the region of stability using the approach of Hindmarsh et al. [16] and Sousa [17]. We consider the case when and and . When , (49) gives From the Von Neumann stability analysis, the scheme is stable if and only if . Thus, we get

On simplification, we get which reduces to

When and , we use the following approximations: We consider (49) and use the approximations in (54) to obtain

Thus LOD Lax-Wendroff is stable if . Therefore, we have

On combining (53) and (56), we obtain the region of stability for the LOD Lax-Wendroff procedure as We choose and [3]. For from (57), we have On solving for , we get

Therefore, for , the stability region for the LOD Lax-Wendroff procedure is . We choose such that is an integer as for our numerical experiments, . Then gives the number of time intervals. We choose , and .

For , we have and . When and , using (51) we get which gives Therefore, for , LOD Lax-Wendroff is stable if . We then choose some values of for the numerical experiments. for our numerical experiments.

3D plots of the modulus of the amplification factor versus phase angles in - and -directions for some different values of and are depicted in Figures 1 and 2. 2D plots of the modulus of the amplification factor versus , when , are shown in Figure 3. The scheme is in general less dissipative at as compared to . Out of the eight combinations of values of and , the scheme is least dissipative when and and it is most dissipative when and .

Figures 4 and 5 show the 3D plots of the relative phase error versus versus for some different values of and . Figure 6 shows the 2D plots of the relative phase error versus , for the case at and . For , we observe phase lag behaviour at and and phase lead behaviour at and . For , we have phase lag behaviour at , and . The scheme is least dispersive when ; and ; . In the following section, we consider the LOD scheme.

5. LOD Explicit Procedure

In this procedure the following approximations are used in the first half time step [3]: On substituting (62) into (14), we get the following finite difference equation for the first half time step: where

The following approximations are used in the second half time step: Then on substituting (65) into (15), we get the following difference equation for the second half time step: where

The complete LOD scheme is given by The Taylor series expansion of the terms on the right hand side of (68) about is given as follows: for and .

We obtain the following modified equation: The scheme is essentially dispersive as the leading error terms are dispersive in nature due to the presence of the odd-order derivative terms and . Also, the scheme is consistent and it is fourth order accurate in space.

The amplification factor of the LOD scheme is given by We consider and . We use the Von Neumann stability analysis and the approach of Hindmarsh et al. [16] to obtain the stability region. When , on simplification of (71), we get For stability, we must have We consider (71) and for and we use the following approximations: We thus have On neglecting higher order terms, we have Thus, the numerical method is stable if . Therefore, we have When , we have and . Using (73) and (77), we get respectively. On combining (78) and (79), we get the stability region when for the LOD procedure as When , we have and and the modulus of the amplification factor is given by From (77) and (81), we obtain the stability region for the LOD procedure when as We analyse the spectral analysis; for , we choose and for , we choose and 0.05.

3D plots of the modulus of the amplification factor versus phase angle in -direction versus phase angle in -direction at two values of , namely, 0.025 and 0.05, at some different values of are shown in Figures 7 and 8. 2D plots of the modulus of the amplification factor versus when are illustrated in Figure 9. The scheme is in general less dissipative at as compared to . Out of the eight combinations of and values, the scheme is least dissipative when and .

Figures 10 and 11 show the 3D plots of the relative phase error versus versus for some different values of and . Figure 12 shows the 2D plots of the relative phase error versus , for the case at and , respectively. We observe that the scheme is least dispersive when and and most dispersive when and . We have phase lead when and and there is phase lag for the other seven combinations of and . In the following section, we present the results of the numerical experiment described in Section 3.2.

6. Numerical Results

6.1. LOD Lax-Wendroff Procedure

The results of the numerical experiment at some values of using the LOD Lax-Wendroff procedure at at and are shown in Tables 1 and 2, respectively. We observe that the dispersion error is significantly greater than the dissipation error. Out of the five combinations of values of and , we observe that the dispersion error and total mean square error are both least when and and are both greatest when and .

6.2. LOD Procedure

The results of the numerical experiment using the LOD procedure at for and are shown in Tables 3 and 4, respectively. Out of the 12 combinations of values of and , we observe that the dissipation, dispersion error, and error rate are all least when and . Also, the dissipation and dispersion errors and error rate are greatest when and . We observe that at the total mean square error, dispersion error, and dissipation error are not much affected by the values of . Also, the dispersion error is greater than the at all values of and considered.

7. Optimization

In this section, we obtain the optimal value of at that minimizes the dispersion error for the two time-splitting procedures.

Since the partial differential equation we consider is slightly dissipative and also we observe from the numerical experiments carried out that the dissipative errors are much less than the dispersive errors, we choose to minimize the square of the dispersion error of the two splitting schemes.

7.1. Proposed Techniques of Optimization

Tam and Webb [18], Bogey and Bailly [19], and Hixon [20] among others have implemented techniques which enable coefficients to be determined in numerical schemes specifically designed for computational aeroacoustics. We now describe briefly how Tam and Webb [18] and Bogey and Bailly [19] define their measures and consequently their technique of optimization in computational aeroacoustics.

The dispersion-relation-preserving (DRP) scheme was designed so that the dispersion relation of the finite difference scheme is formally the same as that of the original partial differential equations. The integrated error is defined as where the quantities and represent the numerical and exact wavenumbers, respectively. The dispersion error and dissipation error are calculated as and , respectively.

Tam and Shen [21] set as 1.1 and optimize the coefficients in the numerical scheme such that the integrated error is minimized.

Bogey and Bailly minimize the relative difference between the exact wavenumber and the effective/numerical wavenumber and define their integrated errors as or

In computational fluid dynamics for a particular method under consideration, the dispersion error is calculated as

We have modified the measures used by Tam and Webb and Bogey and Bailly in a computational aeroacoustics framework to suit them in a computational fluid dynamics framework [22] such that the optimal parameter can be obtained. We have defined the following integrated errors integrated error from Tam and Webb (IETAM) integrated error from Bogey and Bailly (IEBOGEY) [22] as follows:

In [8], the integrated error for a scheme discretising the 1D advection-diffusion equation is obtained as The value of was fixed as 0.02 and the range of values of , was determined. Then the integrated error, which is a function of , was minimized and the optimal value of is determined using NLPSolve function in maple.

We extend the work on optimization of parameters in [8] which is for the 1D advection-diffusion equation for the case of the 2D advection-diffusion equation. We define the integrated error as We first obtain an expression for the RPE of the LOD Lax-Wendroff when discretizing the equation

Since and and we choose , we have and . Hence, the RPE is a function of , , and . Since we can have phase wrapping, we make use of Taylor’s expansion to obtain an approximation for the RPE up to the terms

The integrated error is obtained by using Simpson's method and it is a function of only. A plot of the integrated error versus is shown in Figure 13(a) for . Using the NLPSolve function in maple, this optimal value of is found to be 0.009593 and the minimum value of the integrated error is .

To validate our results, we perform the same numerical experiment described in Section 3.2 at and use the optimal value of or a value of close to this optimal value, in that case, , and compute the errors. The error rate, total mean square error, and dispersion error are , , . These three errors are all least as compared to when other values of are used as shown in Table 1.

We adopt the same procedure to compute the optimal value of for the LOD scheme when . We obtain an approximate expression for the RPE of the LOD when and . We use Simpson's rule to approximate the integral given by (90) which is a function of . A plot of the integrated error versus for is shown in Figure 13(b) and using NLPSolve function in maple the optimal value of is 0.013782 and also the minimum value of the integral is .

We perform the numerical experiment described in Section 3.2 with which is close to the optimal value of we have obtained with . The error rate, total mean square error, dispersion error, and dissipation error are , , , and , respectively. The total mean square error and dispersion error are both least when .

8. Conclusion

In this paper, two time-splitting procedures are used to solve a 2D advection-diffusion equation with constant coefficients when the advection velocity in both - and -directions is 0.8 and also when the coefficient of diffusivity in both - and -directions is 0.01. We perform a stability analysis and spectral analysis of the dispersion and dissipation properties of the two schemes at some values of and . Numerical experiments are carried out and various errors are computed. These errors are dependent on the values of and . It is observed that in general the dispersion error is more affected by the values of and for the LOD Lax-Wendroff scheme as compared to that of the LOD scheme at a given value of . We then use an optimization technique based on minimisation of the square of the dispersion error to find the optimal value of when is chosen as 0.025 and this is validated by numerical experiments.

Future extension of this work to consider other types of advection-diffusion equations when dissipation dominates and to find out which optimization techniques are suitable in these cases. Also, the work can be extended to 2D nonlinear convection-diffusion problems.

Nomenclature

:
: Spatial step
: Time step
: Reynolds number
RPE: Relative phase error per unit time step
AFM: Modulus of amplification factor
: Advection velocity in -direction
: Advection velocity in -direction
: Advection velocity in -direction
: Coefficient of diffusivity in -direction
: Coefficient of diffusivity in -direction
: Coefficient of diffusivity in -direction
: Wavenumber in -direction
: Wavenumber in -direction
: Phase angle in -direction
: Phase angle in -direction
:
:
:
:
:
:
DISP. ERROR: Dispersion error
DISS. ERROR: Dissipation error.

Acknowledgments

The work of Dr A. R. Appadu was funded through the Research Development Programme of the University of Pretoria and the period of funding is from January 2013 to January 2014. Mr Hagos is grateful to the African Institute for Mathematical Sciences (AIMS), South Africa, for a full bursary for the structured Master's program from August 2012 to July 2013. The authors are grateful to the two anonymous reviewers for their comments which were useful in clarifying and focusing the presentation.