Computational Solutions of Two Dimensional Convection Diffusion Equation Using Crank-Nicolson and Time Efficient ADI

To develop an efficient numerical scheme for two-dimensional convection diffusion equation using Crank-Nicholson and ADI, time-dependent nonlinear system is discussed. These schemes are of second order accurate in apace and time solved at each time level. The procedure was combined with Iterative methods to solve non-linear systems. Efficiency and accuracy are studied in term of 2 L , L∞ norms confirmed by numerical results by choosing two test examples. Numerical results show that proposed alternating direction implicit scheme was very efficient and reliable for solving two dimensional nonlinear convection diffusion equation. The proposed methods can be implemented for solving non-linear problems arising in engineering and physics.


Introduction
In this paper we have extended our previous approach associated to two dimension Convection-diffusion equation.The great Physicist Johannes Martinus Burgers discovered Burgers equation, which is non-linear parabolic partial differential equation (PDE) and widely used as a model in many engineering problems, which explains such as physical flow phenomena in fluid dynamics, turbulence, boundary layer behavior, shock wave formation, and mass transport [1].Two dimensional convection-diffusion equation is given by the following equation.
( ) , , , , u f f g g are given sufficiently smooth functions and ( ) , , u x y t may represent heat, diffusion, etc. Re is the Reynolds number.
This equation established the interaction between the non-linear convection processes and the diffusive viscous processes [2].As Burgers equation is probably one of the simplest non-linear PDE for which it is possible to obtain an exact solution [3].Also depending on the magnitude of the various terms in Burgers equation, it behaves as an elliptic, parabolic or hyperbolic PDE, consequently, it is one of the principle model equations used to test the accuracy of new numerical methods or computational algorithms [4].It is widely known that non-linear PDEs do not have precise analytic solutions [5].The first attempt to solve the Burgers equation analytically was done by Batman [6], who derived the steady state solution of this equation as a test solution to one dimension, which was used to model turbulence nature of the phenomena [7] [8].The two dimensional non-linear Burgers equations are a special form of in compressible Naiver-Stokes equations without the pressure term and the continuity equation [9] [10].Due to its wide range of applicability, several researchers, both scientists and engineers, have been interested in studying the properties of the Burgers equation using various numerical techniques [11].They have successfully used it to develop new computational algorithms and to test the existing ones [12].Vineet and Tamsir [10] used two different test problems to analyses the accuracy of the Crank Nicholson(CN) scheme [10].From literature review, it came to know that Newton's method is also applicable to reform the Jacobean matrix to get the linear algebraic sparse matrix.Solution of such algebraic system of equations can be found by Gauss elimination with partial pivoting technique [13] [14].Bahadir also used same technique to test the accuracy of scheme, using fully implicit finite difference scheme [15] [16].The terminology of the Burgers equation explains that with viscous term the Burgers Equation is parabolic while without viscous term it is hyperbolic.In the later case it possesses discontinuous solutions due to the non-linear term and even if smooth initial condition is considered the solution may be discontinuous after finite time [8].It also governs the phenomenon of shock waves [12].
Many different researchers used Burgers equation to develop new algorithms and to test various existing algorithms [4].For exact solution of such non-linear problem, researchers used Hopf-Cole transformation to linearize the Burgers equations into parabolic partial differential equation [17].Some of the researchers also tried to tackle the non-linear Burgers equation directly (without Hopf-Cole transformation), by applying Crank-Nicholson finite difference method to the linearized Burgers equation by Hopf-Cole transformation which is unconditionally stable and is second order convergent, in both space and time with no restriction on mesh size [10].In another result due to Kutluay et al. a direct approach via least square quadratic B-spline finite element method is discussed [18].Recently Pandey et al. discussed Douglas finite difference scheme on linearized Burgers equation which is fourth order convergent in space and second order convergent in time [1] [18].

Problem 1
From literature review, [19] we found that earlier work done by Mittal,Jain and Holla in 2012 [20] on convection-diffusion equation in one dimension.We extended our work to enhance our knowledge towards two dimensional Convection diffusion equation.Two test problems were taken to understand the numerical solution with finite difference schemes.By setting some parameters with arbitrary constants in bounded domain ( )

{ }
, : 0.5 0.5, 0.5 0.5 . Exact solution of the above two dimensional equation is where ( ) and R is a parameter, known as Reynolds number.
Boundary conditions and initial conditions can be taken from exact solution of u(x,y,t) [21].

Problem 2
In this problem the rectangular domain of two dimensional nonlinear convectiondiffusion Equation ( 1) is given as ( ) where ( ) and R is a parameter, known as Reynolds number.
Boundary conditions and initial conditions can be taken from exact solution of ( ) , , u x y t .where Ω is a rectangular domain in 2  R .The main objective of the paper is to find efficient solution of unknown ( ) , u x t .Two test problems were described to understand the numerical solution by taking two finite difference schemes.Also Convection diffusion equation has been extensively studied to describe various kinds of phenomena which can be seen from equation [21].

Second Order Implicit Scheme
We apply Crank-Nicholson implicit finite difference scheme to equation [21], by integrating Equation (1) in the compact way: , when substitute these terms in to Equation ( 1), the Crank-Nicholson Scheme is given by where , 8 2 The scheme shows that the accuracy is of ( ) O k h + .A Jacobian matrix is now Penta-diagonal, but unfortunately due to large number of iterations it extends from the diagonal at least n entries away in every direction,but another methods which can be used to handle such problems (discussed later), because of the large bandwidth, increasing grid points the calculation become more difficult.To overcome this difficulty another method solution is needed.Newton method is used for solving nonlinear task (discussed later).The Crank-Nicholson is computationally inefficient.

Computationally Efficient Implicit Scheme
In search of a time efficient alternate, we analyzed that the Crank-Nicholson scheme for the two dimensional equation, and find out that scheme is not time efficient [8] [11] [12].To get high time efficiency, the common name of Alternating Direction Implicit (ADI) method, can be used [22].American Journal of Computational Mathematics In this approach, the finite difference equations are written in terms of quantities at two levels However, two different finite difference approximations were used alternately, one to advance the calculations from the plane n to a plane n*, and the second to advance the calculations from (n*)-plane to the (n + 1).Same parameters were used in this method as described above.The derivation of ADI scheme, we have following steps; 2 2 O k h + , newton's iterative method is used to solve tridiagonal system.
The family of linear system in x-direction as: where The family of linear system in y-direction as: ( ) , similarly for the reaction term in y-direction: , , finally the scheme makes tridiagonal family of linear system.Iterative methods was carried out to solved this system.The trick used in constructing the ADI scheme is to split time step into two, and apply two different stencils in each half time step, therefore to increment time by one time step in grid point, we first compute both of these stencils are chosen such that the resulting linear system is tridiagonal [5] [22].To obtain the numerical solution, we need to solve a non-linear tridiagonal system at each time step.We have done this by using Newton's iterative method Algorithm 1 To construct Newton iterative method for the two dimensional Convectiondiffusion equation.The non-linear system in equations [23] and [24], can be written in the form: where , , , G G G −  were system of nonlinear equations obtained from the system in [23] and [24].The system of equations, is solved by Newton's iterative method using the following steps 1) Specify ( ) 0 u as an initial approximation.
• Solve the linear system Jacobian matrix, which is computed analytically and is the correction vector.In the iteration method solution at the previous time step is taken as the initial guess.Iteration at each time step is stopped when with Tol is a very small prescribed value.The linear system obtained from Newton's iterative method, is solved by Court's method.

Algorithm 2
Clearly, the system is tridiagonal and can be solved with Thomas algorithm.
The dimension of J is l m × .In general a tridiagonal system can be written as, , with 0 above system can be written as in a matrix-vector form,

Ju S =
where J is a coefficient matrix (Jacobean Matrix), which is known, comes from Newton's iterative method.Right hand side is column vector which is known.Our main goal is to find the resultant vector u .Now we have , , , , , , , , technique is explained in the following steps,

J LU =
where By equating both sides of the Ju S = , we get the elements of the matrices L and U .The computational tricks for the implementation of Thomas algorithm are shown in results, taken from a specific examples.

Error Norms
The accuracy and consistency of the schemes is measured in terms of error norms specially 2 L and L ∞ which are defined as: where ( )

Results and Discussion
Numerical computations were performed using the uniform grid.In Table 1 & Table 2 Crank-Nicholson and ADI results was performed, compared with the analytical results at grids size 20 × 20 in problem 1. Moreover by fixing some parameter such as time step k = 0.0001 time level t = 0.5 and Reynolds number Re = 100,500 with different mesh points.The obtained results was compared with the existing results in literature.By describing results same parameters, notations was used as other researchers were used in their studies.For convergence 2 , L L ∞ norm were used for the unknown ( ) , , u x y t Khater et al.  4 more refine results was obtained with Reynolds no Re = 50.In Figure 1 analytical and numerical results were compared at grid size 20 × 20 k = 0.0001 time = 0.5 and time level = 0.0001 corresponding results in Table 1. Figure 2 and Figure 3 changing time = 0.1 and grid size 30 × 30 with same time level and Reynolds no as in Figure 1.In both Figure 1  boundary and initial conditions were taken from the exact solution showed stable results time step and increasing grid size (refine mesh size).In Table 5 analytical and numerical results of Crank-Nicholson and ADI were compared at gird size 30 × 30, Re = 200, time = 1, time step = 0.001.In problem 2, Table 6 by reducing grid size 20 × 20, increasing time step k = 0.0001 and time = 3 with fixed Reynolds number attained stable results.
In Figure 4           Reynolds no Re value.Solution profiles at t = 0.5 and t = 1 have been presented in Figure 6 and Figure 7 for grids sizes 20 × 20 and 30 × 30 at k = 0.001, Re = 50 and 100 respectively.Figure 8 Crank-Nicholson results obtained using same parameters as in Figure 7 with grid size 25 × 25.Solution presented in Figure 9 at grid point 20 × 20, time = 0.5, step size = 0.0001 with large Reynolds no Re = 200 error increased.Solution reported in Table 11 by changing Reynolds no Re and grid sizes at time = 0.1 and k = 0.001, 0.0001 2 , L L ∞ norm increased.Root mean square value also increased by changing Reynolds no Re.Among all the interior grid points found both 2 , L L ∞ norm of the numerical solution.          .Moreover, got similar results with small grid size.From graphical illustrations, obtained numerical results give steady state solution and the scheme is stable.In Table 9 & Table 10 we observed that results becomes better by changing grid sizes from low to high and by reducing step size.
The obtained results gives excellent agreement with the solutions available in the literature.When the number of grid points get larger than several hundred, the memory and storage of the Crank-Nicholson starts to become a serious issue and it is better to solve this method using different approaches that take advantage of the special form.This method is computationally inefficient.
Thomas algorithm avoids having to store having the whole matrix J (Jacobean) in its memory and solve the system much more expediently.ADI methods reduces to the CN scheme and this method solve the system very efficiently.The order of truncation error: ( ) O k h + .The implementation of ADI computationally is in a time efficient manner.Alternating direction implicit method is fastest when it works and it works well for simple,ideal problems and give efficient results.

Conclusion
In this research work, finite difference methods has been discussed for solving two dimensional convection-diffusion equation.Two test problems were considered, explained the efficiency, accuracy and stability of the schemes.The numerical results showed that Alternating Direction Implicit method is easy to implement and excellent in time efficient manner.The accuracy and stability of these methods were compared to the other numerical methods, shows good agreement with the exact solution.Both ADI and Crank-Nicholson are unconditionally stable and highly accurate.For convergence 2 L and L ∞ norms were treated towards zero when grid size was increased.Numerical results showed that both methods are good but ADI method is consistent and timeefficient.The approach used in this paper may be useful to solve higher dimensional partial differential equations appearing in various applications of science and engineering.
tridiagonal matrix and the array 1 1 , b d depends on l and m The reaction term is x-direction and

Table 7 &
Table 8 error norm reduced when changing step

Table 5 .
Comparison of Analytical and Exact solution at different time level with fixedReynolds number Re = 200 at t = 1, k = 0.001 and grid size = 30 × 30 for unknown

Table 6 .
Comparison of Analytical and Exact solution at different time level with fixed

Table 8 .
Calculating errors using different parameters for unknown values u(x, t) at different grid size and time step k at t = 0.25,

Table 9 .
Calculating Errors using different parameters for unknown values ( )

Table 10 .
Calculating Errors using different parameters for unknown values ( ) and grid size with fixed Reynolds no Re = 1, t = 0.05 and 0.25.Good results obtained when compared the values of this exact solution with those of the approximation gained in Tables1-3.Furthermore in Table9error reduced when changing grid size at fixed Reynolds no Re = 100,k = 0.005, time = 0 .In our next computations in Table10changing Reynolds no Re and grid size at time = 0.1 error increased.It means that error increase by choosing high size

Table 11 .
Calculating errors using different parameters for unknown values ( ) , u x t at different Reynold's number, grid sizes at t = 0.1 Problem 1.