- About this Journal ·
- Abstracting and Indexing ·
- Aims and Scope ·
- Annual Issues ·
- Article Processing Charges ·
- Articles in Press ·
- Author Guidelines ·
- Bibliographic Information ·
- Citations to this Journal ·
- Contact Information ·
- Editorial Board ·
- Editorial Workflow ·
- Free eTOC Alerts ·
- Publication Ethics ·
- Reviewers Acknowledgment ·
- Submit a Manuscript ·
- Subscription Information ·
- Table of Contents

Mathematical Problems in Engineering

Volume 2013 (2013), Article ID 634657, 20 pages

http://dx.doi.org/10.1155/2013/634657

## Time-Splitting Procedures for the Numerical Solution of the 2D Advection-Diffusion Equation

^{1}Department of Mathematics and Applied Mathematics, University of Pretoria, Pretoria 0002, South Africa^{2}African Institute for Mathematical Sciences (AIMS), 6 Road Melrose Road, Muizenberg 7945, Cape Town, South Africa

Received 5 June 2013; Accepted 12 August 2013

Academic Editor: Waqar Khan

Copyright © 2013 A. R. Appadu and H. H. Gidey. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### 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