A mass-conservative higher-order ADI method for solving unsteady convection–diffusion equations

In the paper, a high-order alternating direction implicit (ADI) algorithm is presented to solve problems of unsteady convection and diffusion. The method is fourth- and second-order accurate in space and time, respectively. The resulting matrix at each ADI computation can be obtained by repeatedly solving a penta-diagonal system which produces a computationally cost-effective solver. We prove that the proposed scheme is mass-conserved and unconditionally stable by means of discrete Fourier analysis. Numerical experiments are performed to validate the mass conservation and illustrate that the proposed scheme is accurate and reliable for convection-dominated problems.


Introduction
At the present time, the population explosion, the high standard of living, and the industrial expansion have led to the rise of pollutants into the aquatic environment. In eliminating or at least reducing this pollution, the pollutant transport processes induced by the advection-diffusion equation need to be well studied. Moreover, the processes to be committed should be adapted to the nature of these processes [1]. The following is the twodimensional (2D) mathematical expression of the unsteady advection-diffusion equation without the source term: u tau xxbu yy + pu x + qu y = 0. ( Here, u is an unknown function. The constants p and q are nonnegative speeds of convection, and the constants a and b are positive diffusivity, respectively, in the x-and ydirection. Engaging a conjunction of convection and diffusion, these processes are encountered in science and engineering phenomena very frequently. Also, these problems appear in a number of applications, for instance, in ground water pollutants and air flow, the transportation of oil reservoir, semi-conductor modeling, and so on. By convection, we refer to an actual behavior with some characteristic carried by the flow's ordered motion. On the other hand, diffusion means a physical process with the property transported by random motion of the fluid's molecules. The behavior of fluid subjected to mass, vorticity, and heat transfer can be investigated using the mathematical model formulated by the conservation laws of momentum, mass, and energy. Difficulty in the development of a numerical method is treating the first derivative term that contains large coefficients (convection dominated) in a convection-diffusion equation. It follows that a numerical method that simultaneously solves both hyperbolic term (advection) and parabolic term (dispersion) is needed. However, there is no numerical method that can completely solve this problem [2]. Therefore, every effort has been put to develop the efficient and stable numerical techniques. To solve the convection-diffusion equation, a huge body of research on proficient numerical methods, such as the finite element method (see, e.g., [3][4][5]), the finite difference method (see, e.g., [6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21][22][23]), and others, has been created. We now concentrate on the research about the finite difference method. According to ease of implementation as a conventional approach, the finite difference method is always utilized to solve partial differential equations. High-order compact (HOC) finite difference methods happen to be one of the most commonly used numerical methods to solve the unsteady convection-diffusion equation, that may extremely enhance the approximate precision (see, e.g., [6][7][8][9][10][11][12]).
Since the objective of the alternating direction implicit (ADI) method is decreasing multi-dimensional problems to a series of one-dimensional (1D) problems, and it is only required to solve the tridiagonal systems, the ADI method is greatly effective in solving parabolic and hyperbolic initial boundary value problems. Presented by Peaceman and Rachford [13], the ADI scheme has been considered to be one of the most approved means to solve the higher-dimensional problems due to its unconditional stability and huge efficiency. Nevertheless, the Peaceman-Rachford ADI scheme appears to be second-order accurate in space, and it can generate significant dissipation and phase errors. For gaining resolutions with higher order, many efforts are attempted to use HOC schemes for spatial estimations of Eq. (1). To solve unsteady convection-diffusion equations, a highorder compact ADI scheme was presented by Karaa and Zhang [14]. On that account, great accuracy and great computational efficiency are coincidentally achieved. For 2D unsteady convection-diffusion equations, a Padé scheme-based ADI (PDE-ADI) model was presented by You [15]. It has finer phase and amplitude characteristics. Tian and Ge presented [16] an exponential high-order compact alternating direction implicit (EHOC-ADI) method used to solve 2D unsteady convection-diffusion equations. It highly accomplishes solving convection-dominated equations with high Reynolds numbers. Subsequently, the logical HOC scheme with the ADI (RHOC-ADI) method to solve unsteady convection-diffusion equations was investigated by Tian [17]. On this account, its great effectiveness on solution accuracy and computational efficiency was demonstrated. All of those methods are second-order accurate in time and fourth-order accurate in space (see more [18][19][20]). Additional developments on a series of ADI methods have recently been completed. Some researchers are using the higher-order difference schemes with ADI methods [21][22][23] for achieving sixth-order accuracy in space.
In this paper, we consider the 2D unsteady advection-diffusion equation with the homogeneous boundary conditions and the initial condition where φ(x, y) is a known smooth function. The solution is supposed to have the following asymptotic values: |u(x, y)| → 0 as |x|, |y| → +∞. For that reason, if x L 0, y L 0, x R 0, and y R 0, the initial-boundary value problem is in agreement with the Cauchy problem of Eq. (1). From the past, many conservative schemes have been presented by researchers for examining solutions of various nonlinear partial differential equations. The conservativeness is obviously one of the necessary characteristics of fluid or propagating waves. Confirmed by schemes of the energy or mass conservation, the accomplishment for a numerical approximation in the long-time behavior is one example. A mass-conservative finite difference scheme that conserves the solution for the convection-diffusion equation is therefore required. In the same way, the significant feature which the convection-diffusion problems occupy is the balance of mass. That is, if there exists no source, the mass must be conserved. Some of conservative methods shown in [24][25][26][27][28] preserve mass balance identity well. At this point, work on the family of ADI methods is not necessary.
According to the benefit of ADI methods and conservative discretization simultaneously, the improvement of the solution accuracy in the spatial domain based on the HOC-ADI method, which is a compact sixth-order scheme that inherits mass preserving properties, will be focused on. The novel HOC scheme for the 1D steady convection-diffusion equation will be first implemented. Then, a scheme for solving the 2D unsteady advectiondiffusion equation based on the ADI technique and an exponential difference operator are introduced. We note that no additional computational cost is needed for the compact scheme. Thus, the compact scheme has better effectiveness than the noncompact scheme does.
The paper is structured as follows: Numerical methods consisting of spatial and temporal discretized methods of the 2D convection-diffusion equations are explained in Sect. 2 in detail. Solvability and the conservative property are proved in Sect. 3. In Sect. 4, the analysis of the linearity of Fourier stability for the proposed ADI scheme is presented.
Various numerical examples and detailed numerical results in 1D and 2D cases of unsteady advection-diffusion equation are provided to justify the theory analysis and show the performance of the present ADI scheme. The behavior of numerical solutions is also discussed when the Reynolds number increases. Eventually, remarkable conclusions are summarized to finalize the research paper.

Description of the numerical methods for 2D convection-diffusion equations 2.1 Notations
In this section, we provide some notations which are used in this paper. For a positive integer N , denote τ = T N , t n = nk, n = 0, 1, 2, . . . , N . Let Ω = Ω × (0, T] be the solution domain which is covered by a uniform grid (x i , y j , t n ) as defined by where h x = (x Rx L )/M x and h y = (y Ry L )/M y are the uniform step sizes in the x and y spatial directions, respectively. Denote u n i,j = u(x i , y j , t n ) and the solution space The following notations are introduced for the simplicity for u n ∈ Z 0 h : and u n ∞ = max i,j |u n ij |. For convenience, we let

Double compact high-order ADI finite difference method (DHOC-ADI)
To develop a new HOC-ADI method, we begin with the elementary 1D steady convection diffusion equation where a is the positive constant conductivity, p is the constant convective velocity, and f is a sufficiently smooth function of x. The technique outline is to combine the HOC-ADI proposed by [14] and Richardson's extrapolation for solving the steady convection diffusion equation. By using the Taylor series expansion, we have First, we differentiate both sides of Eq. (4) with respect to x and obtain Differentiating both sides of Eq. (6), we have Substituting Eqs. (6) and (7) into Eq. (5) and rearranging it, we get Apply the fourth-order difference operator for the second-order derivative, and it yields To this end, again, we differentiate both sides of Eq. (7) with respect to x and obtain Also, differentiating both sides of Eq. (10), we have Similarly, substituting Eqs. (10) and (11) into Eq. (9) and rearranging it, we get By using the Richardson extrapolations, we get Combine the above equation and Eq. (12) using the second-order difference operator for the third-and fourth-order derivative; it yields The truncation error analysis shows that Eq. (13) is a sixth-order scheme for the convection-diffusion equation (4). Scheme (13) may be named double compact high-order (DHOC) finite difference scheme. In fact, Eq. (13) can be symbolically formulated as where Notice that the operator L -1 x has symbolic meaning only. By observation, to derive the higher-order compact scheme for numerical solution of transport problems involving convective and diffusive processes, the symbolic operator approximate technique has been extensively used (for more details, see [14,16,17]. An analogous symbolic sixth-order compact approximation operator can also be applied for the y variable. Therefore, we define the difference operator in the y-direction as follows: This yields when the sixth-order compact difference operators L -1 x A x and L -1 y A y are applied to the unsteady 2D convection-diffusion equation (1). For 2D unsteady convection-diffusion problem, the term f ij can be replaced by -∂u/∂t to obtain where O(h 6 ) refers to O(h 6 x + h 6 y ). To complete the discrete scheme, we apply the semidiscrete approximation technique in [14,16,17]. Using the forward Taylor series development, we get By applying the commutativity of the difference operators and Taylor expansions with Eq. (17), we obtain which yields It is worth to notice that the DHOC-ADI scheme (20)- (21) is second-order accurate in time and sixth-order in space. The solution to the resulting DHOC-ADI schemes can be computed by solving the one-dimensional penta-diagonal matrix. The boundary condition for u * can be obtained by Eq. (21). Therefore, by considering the boundary conditions Eq. (2), it is logical to assume that u * has the same asymptotic value as boundary as Eq. (2).

Solvability
The following lemmas and theorems are well-known and elementary results which are essential for proving the solvability of the new DHOC-ADI scheme (20)- (21).

Lemma 1 For any two mesh functions
Moreover,

Lemma 2 For any mesh function u
Theorem 4 There exist u n , u * ∈ Z 0 h satisfying the difference scheme (20)- (21).
Proof It follows from the original problem that u 0 satisfies scheme (20)- (21). By using the mathematical induction, for n ≤ N -1, assume that there exist u 0 , u 1 , . . . , u n , with corresponding u * in each step, satisfying Eqs. (20)- (21). Next, we prove that there exist u n+1 , u * satisfying scheme (20)- (21). By considering the homogenous form of Eq. (20) taking the inner product of Eq. (22) with u * , we obtain By using the boundary conditions, Lemma 1, and Lemma 2, we see that By using Lemma 3, we have which yields where Eqs. (24) and (25) are applied in Eq. (22). This implies that Hence, it uniquely admits a zero solution satisfying scheme (22). Therefore, u * is uniquely solvable. Repeating the same arguments to the homogenous form of Eq. (21), we can easily obtain the solvability and the uniqueness of u n+1 . This completes the proof.

Mass-preserving properties
Theorem 5 Suppose u 0 (x, y) ∈ Z h 0 (Ω), then the finite difference scheme (20)-(21) is conservative for discrete mass in the sense of Proof Multiplying Eq. (21) with h x h y and summing up for i from 1 to M x -1 and j from 1 to M y -1, we obtain According to the boundary conditions, we see that Hence, Eq. (26) gives Next, we again multiply Eq. (20) with h x h y , summing up for i from 1 to M x -1 and j from 1 to M y -1, we obtain Again, by using the boundary conditions, we obtain Consider the right term of Eq. (28), we see that It is not difficult to see that Hence, Eq. (28) gives where Eqs. (29) and (30) are used. Consider the right term of Eq. (31) and use the boundary conditions From Eqs. (27) and (31), we finally arrive at as desired.

Stability analysis
In this section, we study the stability of the proposed scheme (20)-(21) via the von Neumann method for linear stability analysis. Following the method, we assume that the numerical solution can be expressed by virtue of a Fourier series, whose typical term is where I = √ -1, η n is the amplitude at time level n, and θ x = k x θ x and θ y = k y h y are phase angles with the wave number k x and k y in the x-and y-directions, respectively. Substituting the discrete Fourier mode into Eq. (19), the amplification factor G(θ x , θ y ) = η n+1 /η n is found to be where g x (θ x ) is given by with The other term g y (θ y ) can be defined in a similar way by replacing x, a, p, and h x with y, b, q, and h y , respectively, in the above expression. To obtain the stability, it is sufficient to show that |g x (θ x )| ≤ 1 (and also |g y (θ y )| ≤ 1), which is equivalent to γ 1 γ 2 + γ 3 γ 4 ≥ 0. Simple calculations show that which yields Consider Hence γ 1 γ 2 + γ 3 γ 4 ≥ 0, it follows that |g x (θ x )| ≤ 1 for all θ x ∈ [-π, π]. In a similar way, we also have an inequality for g y (θ y ). We then conclude that the present scheme is unconditionally stable.

Numerical experiment
For verifying the effectiveness and correctness of our theoretical results in the previous sections, several numerical experiments are performed in this section and compared with some existing schemes in the recent literature [14,16,17,[30][31][32]. The comparison of numerical solutions with the exact solutions can prove the accuracy of the method, and then the rate of convergence can be measured by using the ratio Rate = log 2 Error 1 Error 2 , where Error 1 and Error 2 are the norm errors relevant with the grid sizes h and h/2, respectively.

1D convection-diffusion equation
Consider the 1D convection-diffusion equation whose analytic solution is of the form where the boundary condition is from the above analytic solution.
In the simulations, we take the computational domain in the interval [-12, 28] with p = 0.8 and a = 0.1. Adapted into the 1D case, the comparison of numerical results using the implementation of the present ADI scheme and the similar approach of ADI, such as the Karaa and Zhang ADI (HOC-ADI) scheme [14], the Tian and Ge ADI (EHOC-ADI) scheme [16], and the rational higher-order compact ADI (RHOC-ADI) [17], is given in Tables 1 and 2. In Table 1, to neglect the error domination in terms of temporal step size for verifying only the spatial sixth-order accuracy, we choose the value of τ to be sufficiently small as of 0.0001. The final time T is set to be 5. The rate of convergence is approximately 6 for the proposed DHOC-ADI which is in good agreement as to the theoretical prediction. However, for the HOC-ADI, EHOC-ADI, and RHOC-ADI, the rate of convergence is approximately 4 due to the limitation of those schemes. When h reduces to 0.1, both L 2 and L ∞ error norms obtained by the proposed ADI scheme provide more accuracy  than others up to three digits. Especially, the most accurate simulation obtained from our method is with the choices for which the L ∞ error norm is less than 4.80 × 10 -9 , while the most accurate simulation of other methods can be reduced to only 2.43 × 10 -6 . In conclusion, we see that the proposed DHOC-ADI scheme in this paper provided more accurate results in comparison with others.
To verify the rate of convergence in time, one can have the fixed value of the spatial step size h = 0.01. It is clearly seen that all schemes provide the second-order rate of convergence which is close to the theoretical prediction, as recorded in Table 2. However, no significant difference between results obtained by our method and the ones obtained by other methods can be observed in these cases.
Next, we consider the absolute errors on the interval [0, 1] at different locations and make a comparison with those mentioned ADI schemes and other recorded results proposed in [30][31][32] including the cubic B-spline method [30], the weighted finite difference method [31], and the exponential B-spline method [32]. The results are reported in Ta-   Table 3 that the numerical solutions obtained with the methods in [30][31][32] are more accurate than those obtained with HOC-ADI schemes in the case of h = 0.01 and τ = 0.001. However, when the coarse grid size h = 0.1 and the smaller temporal step size τ = 0.0001 are used, the great numerical resolution is obtained by the present ADI method, while other ADI methods give very poor results. That is, the present scheme can give accurate solutions comparable with solutions by the methods in [30][31][32] when the step size τ is small enough even when using the coarse grid size h. We also provide the distribution of the Gaussian pulse in Fig. 1, which is in excellent agreement with the exact solutions, especially the sub-plot (top) in the sub-interval [0, 1]. The corresponding absolute error distribution is presented in the sub-plot (bottom).

2D convection-diffusion equation
Consider the 2D convection-diffusion equation The analytic solution is given, as in [33], by where the boundary and initial conditions can be taken directly from Eq. (38).  In our experiments, we set the computational domain to be (x, y, t) ∈ [-1, 3] × [-1, 3] × (0, T]. First, we employed the present scheme in Eqs. (20)-(21) for a = b = 0.01 and p = q = 0.8 in our simulations. To explore the efficiency of the proposed DHOC-ADI scheme, we provide the comparison of numerical results using the present DHOC-ADI scheme and the original HOC-ADI scheme, which are given in Tables 4 and 5. Besides, to show the accuracy of the scheme proposed here, a temporal step size is chosen to be sufficiently small as τ = 2.5 × 10 -6 so that it does not have impact on the accuracy in space while the fine grid h x = h y = 0.01 is chosen to verify the accuracy in time. The simulations run up to T = 1.25. Numerical results in terms of errors and rates of convergence are listed in Tables 4 and 5. In Table 4, the rate of convergence of the DHOC-ADI scheme in space tends to 6 as reported, when h x = h y = 0.2, 0.1, 0.02, and 0.025, while the rate of convergence in space of the HOC-ADI scheme seems to be 4, the same as the theoretical prediction. In Table 5, the errors and rates of convergence in time are reported, where τ = 0.1, 0.05, 0.025, and 0.0125 are used. From the table, we may see that the convergence rate obtained by the present DHOC-ADI scheme and the HOC-ADI scheme is 2, which is expected to confirm the theoretical prediction. In terms of the grid point number, the computational performance of the new DHOC-ADI scheme is clearly better than that of the HOC-ADI schemes. That is, the results of our DHOC-ADI scheme show improvement over the HOC-ADI scheme. Especially, our DHOC-ADI scheme can reduce errors from the HOC-ADI scheme up to 97% in the case h x = h y = 0.025 and τ = 2.5 × 10 -6 . However, in terms of the temporal grid point number, the L 2 and L ∞ error norms are comparable and no significant improvement is noticed.  (5,9) [34] 0.025 0.0125 1.02 × 10 -3 2.25 × 10 -2 Kalita et al. (9,9) [34] 0.025 0.0125 5.24 × 10 -5 1.19 × 10 -3 Kalita et al. [35] 0.025 0.0125 7.77 × 10 -5 1.69 × 10 -3 Noye and Tan [33] 0.025 0.0125 1.43 × 10 -5 4.84 × 10 -4 We also listed the comparison on errors versus the recorded results [33][34][35] using the average error and the L ∞ norm error with the same parameters a = b = 0.01 and p = q = 0.8 at final time T = 1.25 in Table 6. One can find that the computational efficiency of the DHOD-ADI scheme satisfies in the case h x = h y = 0.025 with τ = 0.0125. Moreover, when the coarse grid h x = h y = 0.05 and τ = 0.00625 are used, the error obtained by present ADI scheme is much better than the that by other ADI schemes and is considered as a comparable accurate solution with the recorded ones. In addition, the most accurate simulation as obtained from our scheme is with the choice h x = h y = 0.025 and τ = 0.0001 for which the average error and L ∞ norm error is less than 4.98 × 10 -6 , which is clearly better than others. Figure 2 presents the numerical solutions and its obtained contours based on the present DHOC-ADI scheme within the sub-region [1,2]  The absolute error distributions and the contour plots are shown in Fig. 3. It can be observed that the maximum error is taken around the peak amplitude of the Gaussian pulse in both cases. However, using coarse grid sizes h x = h y = 0.05 and τ = 0.00625, one can observe that errors occur a lot at the bottom left of the pulse base.
As observed in Table 3, in 1D case, the numerical resolution can be made by choosing a small time step; even coarse grid spacing is used. Also, as reported in Table 6, one can see that the most accurate solution is obtained when small temporal step size is used. Figure 4 shows errors obtained by the present ADI scheme and mentioned ADI schemes as a function of time step. From the simulations, the errors can be improved when the  time step is decreased as we expected. The performance of the DHOC-ADI scheme is clearly better than that of others. However, the error slightly increases as the time step is decreased less than τ = 0.0625, presumably because of the increase in round-off error. Next, the proposed DHOC-ADI scheme is applied to verify the conservation of the numerical model. It results from the present method in case h x = h y = 0.01, τ = 0.025 and h x = h y = 0.05, τ = 2.5 × 10 -6 that the values of Q n at any time t ∈ [0, 1.25] coincide with Theorem 5 as presented in Table 7. The results show that the mobility constant slightly differs from the exact value by less than 5.7 × 10 -15 when h x = h y = 0.01 and τ = 0.025 is used while less than 5.3 × 10 -12 when h x = h y = 0.05 and τ = 2.5 × 10 -6 .
To explore the superiority of the present ADI scheme, the numerical simulation at a high Reynolds number is the major indicator as studied before in [15,17]. Here, three cases of Reynolds numbers are presented: Pe = 100, 1000, and 10,000, corresponding to the convective velocities p = q = 100, p = q = 1000, and p = q = 10,000, respectively. The value of viscosity coefficients is set to be a = b = 0.01. As previously observed, the numerical resolution can be achieved when the temporal step size is sufficiently small as shown in Fig. 4. We then choose the temporal step size of τ = 2.5 × 10 -5 , 2.5 × 10 -6 , and 2.5 × 10 -7 for Pe = 100, Pe = 1000, and Pe = 10,000, respectively. To compare the accuracy of the numerical solution, the spatial step size h x = h y = 0.01 is used.
The numerical solutions against the exact solution in the sub-region 1.2 ≤ x, y ≤ 1.8 for each test carried out are plotted in Figs. 5-7, where the solid line indicates the exact solution and the dash-dot line indicates the numerical solutions. Results of the cases present the moving center of the initial Gaussian pulse from (0.5, 0.5) to (1.5, 1.5). The numerical solutions obtained by the presented DHOC-ADI and the RHOC-ADI schemes produce the solution that is in good agreement with the exact solution for all cases (see Figs. 5-7(a), (d)), while small pulse distortion can be observed by the EHOC-ADI scheme (see Figs. 5-7(c)). However, the HOD-ADI scheme produces a noticeable dissipated solution, which can be observed from Figs. 5-7(b) with the same agreement as [15,17]. Thus, it is worth noting that at the high Reynolds numbers the present scheme produces better results than the HOC-ADI scheme does. Table 8 reports the average error, the L 2 norms error, and the L ∞ norms error using Pe = 100, 1000, and 10,000. The results from the present DHOC-ADI scheme are compared with those from the HOC-ADI scheme, the EHOC-ADI scheme, and the RHOC-ADI scheme. It is clear that our scheme can reduce errors from other ADI schemes. Especially, the present DHOC-ADI scheme obtains over 90% less error than the HOC-ADI scheme does and up to approximately 95% when Pe = 10,000. Again, in terms of errors, this clearly shows that the DHOC-ADI scheme provides more incredible resolution than the standard HOC-ADI scheme. Finally, the computations using the DHOC-ADI scheme

Concluding remarks
In this paper, with the help of the Richardson extrapolation and the technique used in [14], we proposed a new high-order alternating direction implicit (DHOC-ADI) method for solving 2D unsteady convection-diffusion equations which inherit the mass-preserving property. The method is second-order in time and sixth-order in space and the unconditional stability was proved through a discrete Fourier analysis. Numerical studies were carried out in both 1D and 2D cases to demonstrate the accuracy and superiority over some existing schemes. In the 2D case, the DHOC-ADI scheme gave high accurate solutions and preserved the mass exactly. Additionally, it is obvious that moving of the Gaussian pulse by the proposed scheme can be smoothed out at high Reynolds numbers, which is superior to the classical HOC-ADI scheme. Furthermore, the multi-dimensional convectiondiffusion equation can be easily extended by a similar approach.