A High Order Compact ADI Method for Solving 3D Unsteady Convection Diffusion Problems

In this paper, we develop a rational high order compact alternating direction implicit (RHOC ADI) method for solving the three dimensional (3D) unsteady convection diffusion equation. The present scheme, based on the idea of the fourth order rational compact finite difference operator for the spatial discretization and the Crank-Nicolson method for the time discretization, is fourth order accurate in space and second order accurate in time. The solution procedure consists of a number of tridiagonal matrix operations, which makes the computation to be cost-effective. It is shown by means of the discrete Fourier analysis that this method is unconditionally stable. Three test problems are given to demonstrate the performance of the present method. The numerical results show that the present RHOC ADI scheme has higher accuracy and better phase and amplitude error characteristics than the classical second order Douglas-Gunn ADI method [16] and some high order compact ADI methods including the Karaa’s HOC ADI method [26], Cao and Ge’s HOC ADI method [27], and our previous exponential HOC ADI method [28].


Introduction
In this paper, we consider the following 3D unsteady convection diffusion equation: is unknown function in a cubic domain Ω with initial condition 0 ( , , , 0) ( , , ), ( , , ) u x y z u x y z x y z = ∈ Ω , and Dirichlet boundary condition ( , , , ) ( , , , ), ( , , , ) (0, ] u x y z t g x y z t x y z t T = ∈ ∂Ω × where ∂Ω is the boundary of Ω , and 0 u , g and the source term S are given functions with sufficient smoothness. In Eq. (1), a , b and c are nonnegative diffusion coefficients and p , q and r are convective coefficients in the x-, y-and z-direction, respectively. The convection diffusion equation (1) is one of the most important models of mathematical physics, which may be considered as a simplified version of the Navier-Stokes equations and plays a very important role in computational fluids dynamics. It also can be seen in many applications to model the convection and diffusion of various physical quantities, such as mass, momentum, energy, vorticity, etc. [1][2][3]. Finite difference methods have been widely used to solve unsteady convection diffusion equations . Among them, traditional numerical discretization schemes, such as the classical explicit difference scheme (forward difference for time and central difference for space), the classical implicit difference scheme (backward difference for time and central difference for space), and explicit or implicit upwind difference schemes, are either first order or second order accurate in space, and get poor solutions for convection-dominated problems if the mesh is not sufficiently refined. In general, explicit ( , , , ) u x y z t difference schemes usually have a small stability region, which would lead the schemes diverge if the stability conditions are not satisfied. So, very small temporal step size must be chosen to keep computation converge. On the other hand, fully implicit difference schemes show good stability, and relatively big temporal step size can be used. But to get a solution of a large linear system arising from an implicit scheme at each time step, direct methods based on Gaussian elimination or conventional iterative methods, such as Gauss-Seidel, Jacobi, and SOR etc. may be too expensive to be used in practice, especially for higher dimensions.
To overcome these difficulties, there exist at least two alternative strategies. One is to use high order finite difference methods [4-14, 17-23, 25-28]. In the last few years, implicit high order compact difference schemes for solving 1D and 2D unsteady convection diffusion equations have been developed [5][6][7][8][9][10][12][13][14]. These schemes usually have third to fourth order accuracy in space and have a large stability region [5,8,10], and they are even unconditionally stable [7,12,13]. With relatively big temporal step size, they can obtain much more accurate solutions than the classical difference schemes. However, traditional iterative methods still be used to solve the sparse linear system arising from the implicit difference schemes at each time step [5,10]. Biconjugate gradient stabilized method (BiCGStab) is used in [13,14]. This method is more efficient than conventional iterative methods, such as Gauss-Seidel, Jacobi, and SOR etc. However, it is still expensive to be used at each time step, especially for higher dimensional problems. The other strategy is to develop efficient and cost-effective difference algorithms, such as alternating direction implicit (ADI) [11,15,16,18,[20][21][22][23][25][26][27][28] or locally one dimensional (LOD) procedures [29][30][31][32], which are based on reducing high dimensional problems in several space variables to collections of 1D problems and only requiring to solve tridiagonal matrices, are simple to implement and economical to use. It is well known that the classical ADI method developed by Peaceman and Rachford [15], and Douglas and Gunn [16] have been popular due to their computational cost-effectiveness. However, these two ADI schemes, which are second order accurate in space and often produce significant dissipation and phase error [11,21,22], are not ideally suited to deal with the spatial discretization of convection-dominated transport problems. And the Peaceman-Rachford ADI scheme for 3D convection diffusion problems is conditionally stable. To achieve higher spatial accuracy computational cost-effectiveness, recently, there has been a renewed interest in the development of HOC ADI methods for 2D and 3D unsteady convection diffusion equations. Karra and Zhang [20], You [21], Tian and Ge [22], Tian [11], developed, respectively, an HOC ADI method for solving 2D unsteady convection diffusion equations with constant coefficients. They are all second order accurate in time, fourth order accurate in space and unconditionally stable. Very recently, Sun and Li [23] presented a combined compact difference ADI (CCD ADI) method for solving 2D unsteady convection diffusion equations. This method is second order accurate in time, sixth order accurate in space and unconditionally stable. For 3D unsteady convection diffusion problems, Wang and Shen [25] proposed two splitting schemes, which are fourth order accurate for spatial diffusion term and second order accurate in time and a revised ADI scheme which obtains spatial fourth order accurate for convective term. Karaa [26] derived a so called delta formulation of HOC ADI (DHOC ADI) scheme, which is second order accurate in time and fourth order accurate in space. It is shown by using a discrete Fourier analysis that the method is unconditionally stable in the diffusion case. However, it is conditionally stable in the convection diffusion case. Recently, Cao and Ge [27] extended Karra and Zhang's HOC ADI scheme [20] to the 3D case. Ge et al. [28] extended Tian and Ge's EHOC ADI scheme [22] to the 3D case. These two generalized HOC ADI schemes are second order accurate in time, fourth order accurate in space and unconditionally stable. More recently, finite difference method is also employed by some authors to solve multidimensional unsteady convection diffusion reaction equations [35,36].
Among these HOC ADI schemes, Tian [11] proposed a rational HOC ADI (RHOC ADI) scheme for solving the 2D unsteady convection diffusion problems, which is temporally second order, spatially fourth order and unconditionally stable. The main advantage of the RHOC ADI scheme is that it is more accurate than the existing other HOC ADI schemes, such as Karra and Zhang [20], You [21], Tian and Ge [22]. So, in this article, we are aiming at developing an RHOC ADI to solve the 3D unsteady convection diffusion equation. This paper is organized as follows. Section 2 presents the RHOC ADI finite difference method for the 3D unsteady convection diffusion problems; in Section 3, the Fourier (or von Neumann) stability analysis is used to show that the proposed RHOC ADI method is unconditionally stable in the convection diffusion case; in Section 4, numerical experiments for three test problems are performed to validate the effectiveness of the present schemes; finally, Section 5 is devoted to some conclusion.

Rational High Order Compact ADI Scheme
For convenience of deduction, Eq. (1) is rewritten as: We introduce the basic idea of rational high order compact ADI method by using the 1D steady convection diffusion equation as following where a is the nonnegative constant conductivity, p is the constant convective velocity, and ( ) f x is a sufficiently smooth function of x. Using the technique outlined in [11], we can easily derive a three-point fourth order compact scheme for Eq. (5), which is formulated symbolically as 1 4 ( ) and x δ and 2 x δ are the second order central difference operators for the first and second derivatives respectively, x h is the mesh size in the x-direction, and the coefficients 1 α , 2 α and α are given as follows: Pe is the cell Reynolds number in the x-direction It is obvious that Eq. (6) with Eqs. (7) and (8) is a fourth order scheme for convection diffusion equation (5). This scheme is named as an RHOC difference scheme in [11]; i.e., the influencing coefficients of the difference scheme formulation are connected to the rational functions of the coefficients of the differential operators and mesh size.
For convenience, we define the fourth order finite difference operators about y-and z-direction as follows: and y δ , z δ , γ are given as follows: where y y Applying the rational fourth order compact difference operators to the 3D unsteady convection diffusion Eq. (4), it yields the following rational fourth order approximation Rearranging Eq. (13), we obtain In which ijk represents the spatial position of ( , , ) Combining Eq. (16) with Eq. (14), a fourth order difference approximation of Eq. (1) is obtained by ( ) Under the assumptions of the constant convective and diffusive coefficients, the difference operators x A , y A , z A , x L , y L and z L commute with each other, which yields that By using the Taylor expansions, and rearranging it, we have When applied to both sides of Eq. (19) with the difference operator x y z L L L and neglected It is clear that the approximation (20)  ω , 2 ω and 3 ω are nonnegative integers less than or equal to 2. Finally, we introduce two intermediate variables * u and ** u , and apply the D'Yakonov ADI-like scheme [33], then we can get the following RHOC ADI scheme for the 3D unsteady convection diffusion equation with a source term.

Stability Analysis
In the following, we use the Fourier (or von Neumann) method for linear stability analysis of the RHOC ADI scheme (21). Let the numerical solution be expressed by means of a Fourier series, whose typical term is In the same way, we can obtain the other two terms ( )

Numerical Experiments
In this section, three test problems possessing exact solutions are given to demonstrate the performance of the RHOC ADI method. We compare the numerical results of the present method with those of the Douglas-Gunn ADI scheme [16], the Karaa's DHOC ADI schemes [26], Cao and Ge's PHOC ADI scheme [27] and EHOC ADI scheme developed by Ge at al. [28]. We conduct our computations using double precision arithmetic on P4/3.4G dual-core personal computer with 4GB memory.

Problem 1
We first consider a 3D unsteady diffusion problem in the  The initial and Dirichlet boundary conditions are directly taken from this analytical solution. This test problem was used in [26][27][28]. For the 3D pure diffusion equation, the PHOC ADI scheme, the EHOC ADI scheme and the present RHOC ADI scheme are the same, because, actually, they all reduce to the standard fourth order Padé scheme. So, we just use the Douglas-Gunn ADI scheme, the DHOC ADI scheme, and the present RHOC ADI scheme to compute. We consider uniform grids ( with different mesh sizes. The compared quantities are the L ∞ norm errors, 2 L norm errors and convergence rate of the numerical solution with respect to the exact solution. The rate of convergence is defined by L norm errors with the spatial grid sizes h and / 2 h or temporal grid sizes t ∆ and / 2 t ∆ , respectively. In Table 1, we set and halve the spatial grid sizes h from 1/ 5 to 1/ 40 to demonstrate the spatial fourth order accuracy. We see that the present RHOC ADI and the DHOC ADI schemes are fourth order accurate in space, whereas the Douglas-Gunn ADI scheme is second order accurate in space. The present RHOC ADI scheme provides more accurate solution than the Douglas-Gunn ADI scheme or the DHOC ADI scheme. In Table 2, 0.01 h = , 0.2 T = and the change of temporal grid sizes t ∆ from 0.05 to 0.00625 are chosen for the demonstration of temporal second order accuracy. We observe that all three schemes have second order accuracy in time. However, the results of the present RHOC ADI scheme and the DHOC ADI scheme become more and more accurate with the reduction in time step than the ones of the Douglas-Gunn ADI scheme.

Problem 2
To further study the performance of the present RHOC ADI scheme, we apply it to a 3D unsteady convection diffusion problem which is defined in the unit cubic region , with an exact solution given, as in [26][27][28], by The Dirichlet boundary and initial conditions are set to satisfy this exact solution.
In Table 3, the L ∞ norm errors, the 2 L norm errors and the CPU time are shown by using the present RHOC ADI, the Douglas-Gunn ADI [16], the DHOC ADI [26], the PHOC ADI [27], and the EHOC ADI [28] schemes. We fix 0.01 a b c = = = and 0.025  Table 3 show that the present RHOC ADI scheme provides a higher accuracy solution than the Douglas-Gunn ADI scheme and the other HOC ADI schemes. We notice the DHOC ADI scheme is divergent when convective coefficients are larger than a certain number. There is nothing strange since the DHOC ADI scheme is conditionally stable for the convection diffusion problem [26]. In Figure 1 Figure 2. We also notice that the present RHOC ADI scheme can still give much better results than the other HOC ADI schemes. This fact can also be seen by observing the position of contour plots of solution surface. An analysis of Figure 1 and Figure 2 shows that the present RHOC ADI scheme produces a numerical solution in good agreement with the exact solution. However, noticeable phase differences are observed between the other HOC ADI schemes and the exact solution. We also notice that the RHOC ADI, the PHOC ADI, the EHOC ADI and the Douglas-Gunn ADI methods exhibit less CPU time than that of the DHOC ADI method from Table 3. The execution CPU time of the RHOC ADI method is much less than that of the DHOC ADI method. This clearly shows that the RHOC ADI method is the most effective in view of accuracy and time consumption.

Problem 3
We consider a 3D unsteady pure convection problem which is defined with the periodic boundary condition in the unit cubic The initial conditions is given by ( , , , 0) sin( ( )) u x y z x y z π = + + This problem is revised by the 2D pure convection problem considered in [11,34], which was used as a test problem of a sin-surface periodic flow. For this problem, the DHOC ADI, the PHOC ADI, and the EHOC ADI schemes become singular because 0 a b c = = = , so these three HOC ADI schemes cannot be used. We just use the present RHOC ADI and the Douglas-Gunn ADI schemes to compute the numerical solution of this problem. The compared quantities are the L ∞ and 2 L norm errors and convergence rate of the numerical solution with respect to the exact solution.  In Table 4, we set 2 t h ∆ = , 0.2 T = and halve the spatial grid sizes h from 0.2 to 0.025 to qualify the spatial fourth order accuracy. We see that the results of the present RHOC ADI scheme become more and more accurate with the reduction in spatial step sizes, while the ones of the Douglas-Gunn ADI scheme are almost invariable. In Table 5, we fix 0.05 h = and change t ∆ from 0.08 to 0.01 to show the temporal accuracy order. We find that the present RHOC ADI scheme is second order accurate in time, whereas the Douglas-Gunn ADI scheme gives very poor solution.

Conclusion
In this article, a rational high order compact ADI method for solving the 3D unsteady convection diffusion equation is established. By using the fourth order compact method for spatial derivatives and using the second order Crank-Nicolson method for temporal derivative, the present RHOC ADI scheme is fourth order accurate in space, second order accurate in time and only involves three-point stencil for each 1D operator. So, in each ADI solution step, it can be solved by simple tridiagonal Gaussian decomposition and may be used by application of the 1D tridiagonal Thomas algorithm with a considerable amount of saving in computing time. It is shown by a discrete Fourier analysis that the RHOC ADI scheme is unconditionally stable. Numerical experiments for three test problems are performed to demonstrate its performance and to show its superiority over the classical Douglas-Gunn ADI scheme and the other existing HOC ADI schemes.