# Numerical Solution of Advection-Diffusion Equation Using a Sixth-Order Compact Finite Difference Method

Authors: Gurhan Gurarslan; Halil Karahan; Devrim Alkaya; Murat Sari; Mutlu Yasar
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

This study aims to produce numerical solutions of one-dimensional advection-diffusion equation using a sixth-order compact difference scheme in space and a fourth-order Runge-Kutta scheme in time. The suggested scheme here has been seen to be very accurate and a relatively flexible solution approach in solving the contaminant transport equation forPe≤5. For the solution of the present equation, the combined technique has been used instead of conventional solution techniques. The accuracy and validity of the numerical model are verified through the presented results and the literature. The computed results showed that the use of the current method in the simulation is very applicable for the solution of the advection-diffusion equation. The present technique is seen to be a very reliable alternative to existing techniques for these kinds of applications.


---

## Body

## 1. Introduction

Problems of environmental pollution can always be reduced to the solution of a mathematical model of advection diffusion. The unknown quantity in these cases is the concentration,C, a scalar physical quantity, which represents the mass of a pollutant or the salinity or temperature of the water [1]. Advection-diffusion equation (ADE) illustrates many quantities such as mass, heat, energy, velocity, and vorticity [2]. The ADE has been used as a model equation in many engineering problems such as dispersion of tracers in porous media [3], pollutant transport in rivers and streams [4], the dispersion of dissolved material in estuaries and coastal seas [5, 6], contaminant dispersion in shallow lakes [7], long-range transport of pollutants in the atmosphere [8], thermal pollution in river systems [9], and flow in porous media [10].  The advection-diffusion transport equation in one-dimensional case without source terms is as follows:
(1)∂C∂t+U∂C∂x-D∂2C∂x2=0,0≤x≤L,0<t≤T,
with initial condition
(2)C(x,t=0)=C0(x)
and boundary conditions
(3)C(x=0,t)=f0(t),C(x=L,t)=fL(t)or-D∂C(L,t)∂x=g(t),
where t is time, x is space coordinate, D is diffusion coefficient, C(x,t) is concentration, U(x,t) is velocity of water flow, and L is length of the channel, respectively. Here C0, f0, fL, and g are prescribed functions, whilst C is the unknown function. Notice that D>0 and U>0 are considered to be positive constants quantifying the diffusion and advection processes, respectively.It is known that the use of the standard finite difference and finite element method is not effective and often leads to unreal results. For that reason, several alternative methods are proposed in the literature for solving the ADE with high accuracy [11]. These include method of characteristic with Galerkin method (MOCG) [11], finite difference method [12–14], high-order finite element techniques [15], high-order finite difference methods [16–24], Green-element method [25], cubic B-spline [26], cubic B-spline differential quadrature method (CBSDQM) [27], method of characteristics integrated with splines (MOCS) [28–30], Galerkin method with cubic B-splines (CBSG) [31], Taylor-Collocation (TC) and Taylor-Galerkin (TG) methods [32], B-spline finite element method [33], Least squares finite element method (FEMLSF and FEMQSF) [34], Lattice Boltzmann method [35], Taylor-Galerkin B-spline finite element method [36], and meshless method [37, 38].Utility of higher-order numerical methods in solving many problems accurately is required. Lately, a noticeable interest in the development and application of CD methods for solving the Navier-Stokes [39–41] and other partial differential equations [42–46] has been renovated. Narrower stencils are required in the CD schemes, and by a comparison to classical difference schemes, they have less truncation error.  In the current paper, accurate solutions of the ADE are obtained by using a sixth-order compact difference (CD6) [47], a fourth-order Runge-Kutta (RK4) schemes in space and time, respectively.


## 2. The Compact Finite Difference Method

CD methods are very popular in the fluid dynamics community because of their high accuracy and advantages associated with stencils [48]. These methods are efficient for higher accuracies without any increase in a stencil, while traditional high-order finite difference methods use larger stencil sizes that make boundary treatment hard. It can also be noted that the CD schemes have been demonstrated to be more precise and computationally economic. Use of smaller stencil sizes in the CD methods is useful when dealing with nonperiodic boundary conditions.A uniform one-dimensional mesh is considered, consisting ofN points: x1,x2,…,xi-1,xi,xi+1,…,xN with the mesh size h=xi+1-xi. The first-order derivatives of the unknown function can be given at interior nodes as follows [47]:
(4)αCi-1′+Ci′+αCi+1′=bCi+2-Ci-24h+aCi+1-Ci-12h
leading to an α-family of fourth-order tridiagonal schemes with
(5)a=23(α+2),b=13(4α-1).
A sixth-order tridiagonal scheme is obtained by α=1/3,
(6)Ci-1′+3Ci′+Ci+1′=112h(Ci+2+28Ci+1-28Ci-1+Ci-2).
Approximation formulae for the derivatives of nonperiodic problems can be derived with the consideration of one-sided schemes for the boundary nodes. Interested readers are referred to the work of Lele [47] for details of the derivations for the first- and second-order derivatives.The third-order formula at boundary point 1 is as follows:(7)C1′+2C2′=1h(-52C1+2C2+12C3).
The fourth-order formula at boundary points 2 and N-1 is as follows:
(8)Ci-1′+4Ci′+Ci+1′=3h(Ci+1-Ci-1).
The third-order formula at boundary point N is as follows:
(9)2CN-1′+CN′=1h(52CN-2CN-1-12CN-2).
Use of the above formulae leads to following matrix equation:
(10)BC′=AC,
where C=(C1,…,Cn)T. Consideration of the first-order operator twice will give us the second-order derivative terms; that is,
(11)BC′′=AC′.
The RK4 scheme is considered to obtain the temporal integration in the present study. Utility of the CD6 technique to (1) gives rise to the following differential equation in time:
(12)dCidt=LCi,
where L indicates a spatial linear differential operator. The spatial and temporal terms are approximated by the CD6 and the RK4 schemes, respectively. The semidiscrete equation (12) is solved using the RK4 scheme, through the following operations:
(13)C(1)=Cn+12ΔtL(Cn),C(2)=Cn+12ΔtL(C(1)),C(3)=Cn+ΔtL(C(2)),Cn+1=Cn+16Δt[L(Cn)+2L(C(1))+2L(C(2))+L(C(3))].
To obtain the approximate solution of (1) with the boundary and initial conditions using the CD6-RK4, the domain [0,L] is first discretized such that 0=x1<x2<⋯<xN=L where N is the number of grid points.


## 3. Numerical Illustrations

Let us consider the advection-diffusion equation with the initial and boundary conditions. The numerical results are compared with the exact solutions. The differences between the computed solutions and the exact solutions are shown in Tables1–5. Three examples for which the exact solutions are known are used to test the method described for solving the advection-diffusion equation. The technique is applied to solve the ADE with C0(x), f0(t), fL(t), and g(t) prescribed. To test the performance of the proposed method, L2 and L∞ error norms are used as follows:
(14)L∞=maxi|ciexact-cinumerical|,L2=∑i=1N|ciexact-cinumerical|2.Table 1
Peak concentration values att=9600 s for various Cr numbers (Δt=50 s).
Cr
0.25
0.50
0.75
1.00
MOCS [28]
9.677
9.756
9.805
10.000
MOCG [11]
9.816
9.836
9.934
10.000
CBSG [31]
9.986
9.986
9.993
9.986
FEMLSF [34]
9.647
9.864
9.918
9.943
FEMQSF [34]
9.926
9.932
9.949
9.961
TC [32]
9.940
9.984
9.993
9.986
TG [32]
9.989
9.991
9.996
9.991
CD6
9.999
10.000
10.000
10.000
Exact
10.000
10.000
10.000
10.000Table 2
Error norms for various Cr values att=9600 s.
Cr
h
CSDQM [27]
FEMLSF [34]
FEMQSF [34]
CD6
L
2
L
∞
L
2
L
∞
L
2
L
∞
L
2
L
∞
0.1250
200
34.734
1.156
32.874
1.350
12.555
0.518
0.8511
0.4293
0.2500
100
2.685
0.136
10.596
0.494
7.951
0.376
0.0218
0.0100
0.5000
50
0.170
0.008
7.984
0.380
7.908
0.373
0.0024
0.0008
1.0000
25
0.023
0.001
7.881
0.377
7.908
0.379
0.0029
0.0007Table 3
Comparison between numerical solutions and the exact solution.
x
MOCS
MOCG
CBSG
FEMLSF
FEMQSF
TC
TG
CD6
Exact
[28]
[11]
[31]
[34]
[34]
[32]
[32]
Δ
t
=
10 s
Δ
t
=
1 s
0
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
18
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
1.000
19
1.000
0.999
1.000
1.000
1.000
0.999
0.999
0.999
0.999
0.999
20
1.000
0.998
0.999
0.999
1.000
0.999
0.998
0.998
0.998
0.998
21
1.000
0.996
0.996
0.997
0.999
0.999
0.996
0.996
0.996
0.996
22
1.000
0.990
0.991
0.993
0.996
0.998
0.991
0.992
0.991
0.991
23
1.000
0.978
0.981
0.985
0.989
0.994
0.980
0.982
0.982
0.982
24
1.000
0.957
0.961
0.970
0.974
0.987
0.960
0.965
0.964
0.964
25
1.000
0.922
0.927
0.943
0.946
0.972
0.926
0.936
0.935
0.934
26
0.996
0.870
0.874
0.902
0.900
0.945
0.874
0.891
0.889
0.889
27
1.013
0.799
0.800
0.842
0.832
0.902
0.800
0.827
0.824
0.823
28
1.047
0.708
0.706
0.763
0.743
0.838
0.705
0.743
0.739
0.738
29
0.897
0.602
0.596
0.666
0.638
0.755
0.595
0.641
0.637
0.636
30
0.457
0.488
0.479
0.556
0.524
0.653
0.479
0.528
0.523
0.523
31
0.067
0.375
0.366
0.442
0.411
0.541
0.366
0.413
0.408
0.408
32
−0.036
0.272
0.265
0.332
0.306
0.427
0.264
0.306
0.301
0.301
33
−0.010
0.185
0.181
0.235
0.218
0.320
0.181
0.212
0.208
0.208
34
0.002
0.118
0.118
0.156
0.147
0.227
0.117
0.138
0.135
0.135
35
0.000
0.070
0.072
0.096
0.095
0.152
0.072
0.084
0.082
0.082
36
0.000
0.038
0.042
0.055
0.058
0.096
0.041
0.048
0.047
0.046
37
0.000
0.020
0.023
0.030
0.034
0.057
0.023
0.025
0.025
0.024
38
0.000
0.009
0.012
0.015
0.019
0.032
0.012
0.012
0.012
0.012
39
0.000
0.004
0.006
0.007
0.010
0.017
0.006
0.006
0.005
0.005
40
0.000
0.002
0.003
0.003
0.005
0.008
0.002
0.002
0.002
0.002
41
0.000
0.001
0.001
0.001
0.003
0.004
0.001
0.001
0.001
0.001
42
0.000
0.000
0.001
0.000
0.001
0.001
0.000
0.000
0.000
0.000Table 4
A comparison of the peak errors of different solution techniques forPe=4 with h=0.025.
Cr
Δ
t
CN [12]
FD3 [17]
CD6
0.016
0.0005
1.54
E
-
03
9.13
E
-
04
5.64
E
-
09
0.032
0.001
1.51
E
-
03
9.31
E
-
04
5.64
E
-
09
0.064
0.002
1.47
E
-
03
9.73
E
-
04
5.78
E
-
09
0.08
0.0025
1.44
E
-
03
9.96
E
-
04
5.99
E
-
09
0.16
0.005
1.34
E
-
03
1.13
E
-
03
1.10
E
-
08
0.32
0.01
1.20
E
-
03
1.47
E
-
03
8.18
E
-
08
0.64
0.02
1.12
E
-
03
2.26
E
-
03
9.02
E
-
07
0.80
0.025
NA
NA
1.79
E
-
06Table 5
A comparison of analytical and CD6 solutions for various values ofx with h=0.025, Δt=0.005, Cr=0.16, and Pe=4.
x
Exact
CD6
Absolute error
3.50
0.0000000
0.0000000
3.78
E
-
13
4.00
0.0000159
0.0000159
1.02
E
-
09
4.50
0.0201770
0.0201769
1.94
E
-
08
5.00
0.2182179
0.2182179
1.10
E
-
08
5.50
0.0201770
0.0201770
2.51
E
-
08
6.00
0.0000159
0.0000160
9.45
E
-
10
6.50
0.0000000
0.0000000
4.40
E
-
13Example 1.
Here, pure advection equation is considered in an infinitely long channel of constant cross-section and bottom slope, and velocity is taken to beU=0.5 m/s. Concentration is accepted to be the Gaussian distribution of ρ=264 m, and initial peak location is x0=2000 m. The initial distribution is transported downstream in a long channel without change in shape by the time t=9600 s. Exact solution of this problem is as follows [11]:
(15)C(x,t)=10exp[-(x-x0-Ut)22ρ2].
At the boundaries the following conditions are used:
(16)C(0,t)=0,-D(∂c∂x)(9000,t)=0.
Initial conditions can be taken from exact solution. The initial Gaussian pulse at t=0, the concentration distribution obtained using the CD6 solution, and  concentration distribution obtained using exact solution at t=9600 s are compared in Figure 1. Both quantitative and qualitative agreements between the exact and the CD6 solutions are excellent (see Tables 1 and 2, Figure 1).Figure 1
Comparison of the exact solution and the numerical solution obtained with CD6 scheme forh=25 and Δt=50 s.As seen from Table1, the CD6 method has given the closest result to the exact peak concentration value. Note that, since an explicit time integration scheme RK4 is used in this study, the CD6 scheme cannot produce any result for Cr>1. To show the accuracy of the obtained results, L2 and L∞ error norms have been calculated using the CD6 scheme and exhibited in Table 2. As seen in the corresponding table, the CD6 solution is better than its rivals.Example 2.
Flow velocity and diffusion coefficient are taken to beU=0.01 m/s and D=0.002 m2/s in this experiment. Let the length of the channel be L=100 m and be divided into 100 uniform elements. The Pe number is accepted to be 5. The Cr numbers are selected as 0.01, 0.1, and 0.6 for the present work.  Exact solution of the current problem is [11]
(17)C(x,t)=12erfc(x-Ut4Dt)+12exp(UxD)erfc(x+Ut4Dt).
At the boundaries the following conditions are used:
(18)C(0,t)=1,-D(∂c∂x)(L,t)=0.
Initial conditions can be taken from exact solution. Comparison between numerical solutions and the exact solution is given in Table 3. In the calculation of the exact results given by Szymkiewicz [11], there has erroneously been a mistake. Therefore the exact results have been recalculated in MATLAB. As seen in Table 3, the solutions produced by other researchers [11, 28, 31, 32, 34] for Cr=0.6 do not converge enough. This case proves that the selected time step (Δt=60 s) is greater than it needs to be. Note also that the CD6 scheme gives stable results for Cr=0.6, and the computed results are nearly as good as in FEMQSF [34]. Since this problem cannot be solved accurate enough for Δt=60 s, the calculations have been repeated for the cases Δt=10 s (Cr=0.1) and Δt=1 s (Cr=0.01), and the corresponding results have been given in Table 3. The results produced by the CD6 scheme for Δt=1 s are the same as with the exact solution, while the results of the CD6 scheme for Δt=10 s are seen to be acceptable level. Comparison of the exact solution and the numerical solution obtained with CD6 scheme for h=1 m and Δt=10 s is shown in Figure 2. As can be seen in Figure 2, there is an excellent agreement between CD6 and exact solutions.Figure 2
Comparison of the exact solution and the numerical solution obtained with CD6 scheme forh=1 m and Δt=10 s.Example 3.
Consider the quantitiesU=0.8 m/s and D=0.005 m2/s in (1). The following exact solution for this example can be found in [49]:
(19)C(x,t)=14t+1exp[-(x-1-Ut)2D(4t+1)].
At the boundaries the following conditions are used:
(20)C(0,t)=14t+1exp[-(-1-Ut)2D(4t+1)],C(9,t)=14t+1exp[-(8-Ut)2D(4t+1)].
Initial conditions can be taken from exact solution. The distribution of the Gaussian pulse at t=5 s is computed using the exact solution and compared with the concentration distribution obtained using the CD6 solution as shown in Figure 3. As can be seen in Table 4, the CD6 results in Example 3 are far more accurate comparison to Crank-Nicolson (CN) scheme [12] and third-order finite difference (FD3) scheme [17]. In Table 5, a comparison of analytical and CD6 solutions is carried out for various values of x with h=0.025, Δt=0.005, Cr=0.16, and Pe=4.Figure 3
Comparison of analytical and CD6 solutions for transport of one-dimensional Gaussian pulse.


## 4. Conclusions

This paper deals with the advection-diffusion equation using the CD6 scheme in space and the RK4 in time. The combined method worked very well to give very reliable and accurate solutions to these processes. The CD6 scheme provides an efficient and alternative way for modeling the advection-diffusion processes. The performance of the method for the considered problems was tested by computingL2 and L∞ error norms. The method gives convergent approximations for the advection-diffusion problems for Pe≤5. Note that numerical solution cannot obtained while Pe>5 and Cr>1. To overcome these disadvantages, upwind compact schemes and implicit time integration need to be used. For further research, special attention can be paid on the use of compact difference schemes in computational hydraulic problems such as sediment transport in stream and lakes, contaminant transport in groundwater, and flood routing in river and modeling of shallow water waves.



---

*Source: 672936-2013-04-16.xml*