American Journal of Computational Mathematics, 2017, 7, 208-227
http://www.scirp.org/journal/ajcm
ISSN Online: 2161-1211
ISSN Print: 2161-1203
Computational Solutions of Two Dimensional
Convection Diffusion Equation Using
Crank-Nicolson and Time Efficient ADI
Muhammad Saqib, Shahid Hasnain, Daoud Suleiman Mashat
Department of Mathematics, Numerical Analysis, King Abdulaziz University, Jeddah, Saudi Arabia
How to cite this paper: Saqib, M., Has-
nain, S. and Mashat, D.S. (2017) Computa-
tional Solutions of Two Dimensional Con-
vection Diffusion Equation Using Crank-
Nicolson and Time Efficient ADI. Ameri-
can Journal of Computational Mathemat-
ics, 7, 208-227.
https://doi.org/10.4236/ajcm.2017.73019
Received: February 23, 2017
Accepted: July 31, 2017
Published: August 3, 2017
Copyright © 2017 by authors and
Scientific Research Publishing Inc.
This work is licensed under the Creative
Commons Attribution International
License (CC BY 4.0).
http://creativecommons.org/licenses/by/4.0/
Open Access
Abstract
To develop an efficient numerical scheme for two-dimensional convection
diffusion equation using Crank-Nicholson and ADI, time-dependent nonli-
near 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
2L , L∞ norms confirmed by numerical results by choosing two test
term of
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.
Keywords
Crank-Nicholson, Taylor’s Series, Newton’s Iterative Method, Alternating
Direction Implicit (ADI)
1. 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 dif-
ferential 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.
DOI: 10.4236/ajcm.2017.73019 Aug. 3, 2017
u
t
+
uu
x
+
uu
y
−
(
u
xx
+
u
yy
)
=
0
1
R
(1)
208
American Journal of Computational Mathematics
x y t
,
,
where (
(
0,
with initial conditions
∈Ω×
)
T
]
(
u x y
,
,0
)
=
u x y
,
0
(
)
,
(
x y
,
)
∈Ω
The Dirichlet boundary conditions are given by
)
(
u b y t
,
(
)
u x d t
,
f
1
g x y t
,
1
)
x y t
,
,
)
,
(
(
=
=
,
,
,
=
f
2
=
g
2
(
(
,
(
)
u a y t
,
)
(
u x c t
,
]
∈Ω×
T
,
,
{
(
)
(
where (
x y t
0,
,
,
]
2R , (
domain in
smooth functions and
Reynolds number.
)
,
a
:
0,T is the time interval.
Ω =
x y
,
≤ ≤
is a rectangular
x b c
,
≤ ≤
f g g are given sufficiently
u f
,
,
2
0
1
may represent heat, diffusion, etc. Re is the
y
,
(
u x y t
,
)
,
,
1
M. Saqib et al.
x y t
,
,
)
)
x y t
,
,
}
d
2
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, conse-
quently, 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
209
American Journal of Computational Mathematics
DOI: 10.4236/ajcm.2017.73019
M. Saqib et al.
equations into parabolic partial differential equation [17]. Some of the resear-
chers 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].
1.1. 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
. Exact
x
solution of the above two dimensional equation is
y
+ −
2
0.5 tanh
(
u x y t
,
0.5, 0.5
}
0.5
: 0.5
−
,x y ∈Ω ,
where (
t > 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].
≤ ≤
≤ ≤
)
t R
(2)
Ω =
x y
,
{
(
(
x
−
,
)
=
−
y
)
0
)
1.2. Problem 2
In this problem the rectangular domain of two dimensional nonlinear convection-
)
diffusion Equation (1) is given as
. Exact
x y
1
,
solution of the above two dimensional equation is
≤ ≤
≤ ≤
Ω =
1,0
: 0
x
y
(
(
u x y t
,
,
)
=
(
x
1
+
1
y
+ −
2
)
t R
(3)
,
)
)
0
,x y ∈Ω ,
. where Ω is a rectangular domain in
where (
t > 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
2R . The main objective of the
,
(
),
u x t
paper is to find efficient solution of unknown
. 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].
2. Numerical Methods
Numerical solution of the two dimensional non-linear equation in a finite
210
American Journal of Computational Mathematics
DOI: 10.4236/ajcm.2017.73019
M. Saqib et al.
)
(
)
(
=
−
=
hy
and
b a L
−
d c M
domain Ω . The first step is to choose integers L and M to define step sizes
hx
in x and y directions respectively.
xh and the interval [c, d]
Partition the interval [a, b] into L equal parts of width
yh . Place a grid by drawing vertical and horizontal
into M equal parts of width
lines through the points with coordinates (
for each
l
0,1,2,
=
also the lines
for each
x
x=
are grid lines, and their intersections are the mesh points of
l
the grid. For each mesh point in the interior of the grid, (
, for
l
=
, we apply different algorithms to
approximate the numerical solution to the problem in equation [21] also we
assume
where t is the time.
, where
0,1, 2,
)
x y
,l
m
m
=
and
L
−
=
−
c mk
and
and
a lh
nk n
,
x y
,l
= +
= +
lx
M
1,2,
L
y
0,1,
y=
1, 2,
NT
my
M
m
1
=
=
1
)
,
m
,
m
,
,
,
nt
2.1. Second Order Implicit Scheme
We apply Crank-Nicholson implicit finite difference scheme to equation [21], by
integrating Equation (1) in the compact way:
n
1
+
u
l m
,
n
u
l m
,
,
+
k
n
1
+
u
l m
,
u
=
n
1
+
u
l m
1,
+
−
n
1
+
u
l m
1,
−
n
u
l m
,
+
2
n
u
l m
1,
−
,
u
t
=
u
x
=
−
n
u
+
l m
1,
+
h
4
ˆ
u
2
+
l m
,
2
h
ˆ
u
l m
1,
−
ˆ
u
l m
1,
+
−
2
δ
x
ˆ
u
=
,
2
δ
y
ˆ
u
=
=
u
y
ˆ
u
l m
,
n
1
+
u
l m
,
1
+
−
n
1
+
u
l m
,
1
−
n
u
+
l m
,
h
4
ˆ
u
l m
1
,
−
−
1
+
+
ˆ
u
2
l m
,
2
h
−
n
u
l m
,
1
−
1
+
+
n
u
l m
,
when substitute these terms in to Equation (1), the Crank-Nicholson Scheme is
given by
(
R u
1
n
1
+
u
l m
1,
−
k
h
8
){
n
1
+
u
l m
1,
+
n
1
+
u
+
l m
,
1
−
n
1
+
u
l m
1
,
+
k
Rh
n
u
l m
,
n
u
4
l m
,
n
u
l m
,
0
n
1
+
u
4
l m
,
n
1
+
u
l m
,
1
+
n
1
+
u
l m
,
1
−
n
u
l m
1,
+
n
u
l m
1,
+
n
u
l m
1,
−
n
1
+
u
l m
1,
−
n
u
l m
1,
−
n
1
+
l m
1,
+
(4)
n
u
l m
,
n
u
l m
,
n
u
l m
,
1
+
n
1
+
l m
,
R
2
R
1
}
+
−
+
−
+
−
+
−
+
+
−
+
+
+
+
−
=
=
=
2
1
+
1
−
1
−
,
2
−
n
1
+
u
l m
,
R u
2
+
where
(
O k
)
2
2
h+
The scheme shows that the accuracy is of
. 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.
2.2. 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].
211
American Journal of Computational Mathematics
DOI: 10.4236/ajcm.2017.73019
M. Saqib et al.
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;
=
−
*
u
l m
,
n
u
l m
,
Sweep in x-direction
(
*
P u
l m
1,
1
+
(
n
P u
+
l m
,
2
(
P u
3
(
2
P
δ
x
1
n
u
l m
,
*
u
l m
,
n
l m
,
+
−
=
−
*
u
2
+
l m
,
n
u
2
−
l m
,
−
1
+
(
1
+
n
u
l m
,
(
*
u
l m
,
+
*
u
l m
1,
−
n
u
+
l m
1,
−
n
u
l m
,
1
−
+
n
u
l m
,
)
)
)
)
+
n
n
u
u
2
−
+
l m
l m
,
1,
+
(
)
n
n
u
P u
+
l m
l m
1,
3
,
+
(
)
n
u
l m
1,
−
n
u
−
l m
1,
−
)
)
(5)
(
2
δ
y
(
P
2
n
u
l m
,
)
)
+
(
kF u
n
)
(6)
where
n
)
=
−
*
u
l m
,
(
P u
3
n
1
+
u
l m
,
(
F u
Sweep in y-direction
(
n
1
+
P u
l m
1
,
1
+
(
*
P u
+
l m
1,
2
+
(
(
*
*
P u
u
l m
l m
,
3
,
(
(
2
P
δ
y
1
n
1
+
u
l m
,
*
u
l m
,
=
−
+
−
=
where
*
)
=
(
F u
(
P u
3
k
Re h
2
defines similarly to
(
2
O k
P
1
=
*
,l mu
where
The method has accuracy
solve tridiagonal system.
)
)
+
+
)
)
k
Reh
n
δ
l m x
,
(
n
u
l m
,
)
)
+
(
P u
3
n
δ
l m y
,
n
u
l m
,
)
+
n
1
+
u
2
l m
,
*
u
2
−
l m
,
−
1
+
+
n
1
+
u
l m
,
1
−
*
u
+
l m
1,
−
*
u
l m
,
1
−
*
u
l m
,
)
+
*
u
2
−
+
l m
1
,
+
(
*
*
u
P u
l m
l m
1,
3
,
+
*
u
l m
,
1
−
*
u
−
l m
1,
−
(
)
)
)
(7)
n
1
+
u
l m
,
+
*
u
l m
,
)
)
(
2
δ
x
(
P
2
*
u
l m
,
)
)
+
(
kF u
*
)
(8)
*
δ
l m x
,
(
*
u
l m
,
(
P u
3
*
u
l m
,
)
,
P
2
=
2
,
P
3
2
*
δ
l m y
,
k
h
2
=
1
n
l mu + . This method is unconditionally stable.
,
, newton’s iterative method is used to
h+
2
)
The family of linear system in x-direction as:
(
n
F ku
x m
,
a u
1
b
m
1
u
*
x m
,
=
+
where
m
=
1,
−
M
,
1
The family of linear system in y-direction as:
(
*
F ku
y m
,
n
1
+ =
y m
,
d
m
1
c u
1
+
)
)
,
1
1,
L
−
l
=
1,a c develops tridiagonal matrix and the array
1
where
where
m
The reaction term is x-direction
=
)
similarly for the reaction term in y-direction:
(
n
kF u
1,
(
F ku
n
x m
,
)
m
,
kF u −
L
1
n
1,
(
,
(9)
(10)
,b d depends on l and
1
1
)
DOI: 10.4236/ajcm.2017.73019
212
American Journal of Computational Mathematics
M. Saqib et al.
(
F ku
*
y l
,
)
(
*
kF u
1,
=
M
1
−
)
,
,
(
*
kF u
1,1
)
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] [7] [8] [9] [11] [17] [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 Convection-
diffusion equation. The non-linear system in equations [23] and [24], can be
written in the form:
) 0
(
G S =
(11)
where
S
≈
n
1
+
u
=
S
1,
m
,
S
2,
m
,
,
S
(
2
L
−
)
2 ,
m
n
1
+
u
m
1,
,
u
n
1
+
m
2,
,
n
1
+
u
m
3,
,
,
u
n
1
+
L m
1,
−
T
T
1
+
n
mu
:
G −
(
L
2
=
T
=
,
,
,
m
m
G
G G
1,
2,
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
, where
G G
1,
2,
G −
(
L
)
2 ,
)
2 ,
m
m
m
m
,
,
,
2
1) Specify
k =
2) For
( )0u
0,1,2,
as an initial approximation.
until convergence achieve.
)
(
)
R u
= −
u
∆
)
k
k
k
(
)
(
)
(
(
A u
(
k
)
)
(
,
u
)k
(
)
• Solve the linear system
)
(
k
1k
+ =
u
u
• Specify
+ ∆
(
)
is (
m m×
A u
where
Jacobian matrix, which is computed analytically and
)ku∆
(
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
(
)k
R u
with Tol is a very small prescribed value. The linear system
obtained from Newton’s iterative method, is solved by Court’s method.
Convergence done with iterations along less CPU time [5] [14].
Tol
≤
∞
(
)
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,
c
l
above system can be written as in a matrix-vector form,
, with
a x
l
l
b x
l
l
c x
l
l
a
1
+
+
=
=
S
1
+
1
−
l
=
0
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
213
American Journal of Computational Mathematics
DOI: 10.4236/ajcm.2017.73019
M. Saqib et al.
J
0
0
0
0
0
c
2
c
1
b
2
b
0
1
a
0
2
=
a
n
1
−
c
n
1
−
b
n
1
−
b
n
a
n
]
]
n
t
u
s
n
t
u
=
S
=
2
,
[
u u u
,
1
3
[
s s
,
1
2
s
3
,
,
,
,
,
technique is explained in the following steps,
LU=
J
where
and
L
2
0
0
0
0
0
0
0
0
0
µ
1
0
0
β µ
2
0
α β µ
3
3
=
0
0
n
n
1
−
α β µ
n
α β µ
n
1
−
1
−
n
n
3
U
δ λ
1
1
1
0
0
δ λ
2
2
1
0
0
δ λ
3
0
1
0
0
0
0
=
0
...
0
λ
2
0
δ
n
1
−
1
δ
n
−
1
0
1
0
0
3
2
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.
3. Error Norms
The accuracy and consistency of the schemes is measured in terms of error
norms specially
2L and L∞ which are defined as:
)2
(
U
∑
u
i
−
L
i
j
j
,
,
RMSError
=
i
,
j
1
=
L L
×
(
U
i j
.
−
u
i
,
j
)
L
∞
=
L
∑
max
i L j
1
≤ ≤
1
=
(12)
(13)
DOI: 10.4236/ajcm.2017.73019
(
U
ρ=
i j
,
L
2
−
u
i j
,
t
) (
U
i j
,
−
u
i j
,
)
(14)
214
American Journal of Computational Mathematics
M. Saqib et al.
,
)
)
(
)
λ
i j
,
i j
,
−
u
i j
,
i j
,
=
max
(
U
ρ
2,L L∞ norm were used for the unknown
,
. In this method
)
respectively.
denote the numerical and exact solutions at
and λ is
)
(
(
U x y t
u x y t
,
,
where
and
the grid point (
)
x y t
,
,
m n
l
an eigen value of (
u−
U
4. 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
Khater et al.
(2008), Mittal and Jiwari (2012), Kutluay and Yagmurlu (2012) and R.C. Mittal
and Amit Tripathi were considered this problem in (2016) [21] [23] [24].
Approximate results in problem 1 Table 3, comparison of analytical and
numerical results of CN and ADI at fixed Reynolds no Re = 500, time level =
0.0005, time = 1, grid size 30 × 30 with different typical mesh points at (-0.3,-0.3),
(−0.45, −0.45), (−0.35, −0.35), (0.25, −0.25). The obtained results makes a very
good agreement with the exact solution. To attain more refine and better results
by changing time level = 0.001, grid size = 25 × 25 with the same mesh points in
Table 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 and Figure 2 almost same behavior
corresponds to results in Tables 1-4. Similar patterns were depicted by Kutluay
and Yagmurlu (2012), Mittal (2016) [21] [24]. The obtained solutions were
better than those obtained in earlier studies (Khater et al., 2008; Kutluay and
Yagmurlu, 2012, Mittal (2016).
(
u x y t
,
,
)
In problem 2, considering Equation (1) over the domain [0.5, 0.5] × [0.5, 0.5],
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 showed ADI results by reducing Reynolds no Re = 10 with gird
size 20 × 20 k = 0.0001 and time = 0.5 and in Figure 5 increasing grid size 101 ×
101 at Reynolds 50, time level = 0.0001, time = 0.5, obtained good accuracy
2yδ are more complex in
corresponds to the exact results.Matrices
Crank-Nicholson case depend on the order where the grid points are arranged in
to the array u. In Table 7 & Table 8 error norm reduced when changing step
2xδ and
215
American Journal of Computational Mathematics
DOI: 10.4236/ajcm.2017.73019