Numerical Solution of Advection Diffusion Equation Using Semi-Discretization Scheme

Khandoker Nasrin Ismet Ara¹, Md. Masudur Rahaman², Md. Sabbir Alam¹

¹Department of Mathematics, Jahangirnagar University, Savar, Bangladesh.
²Department of Mathematics, Patuakhali Science & Technology University, Dumki Upazila, Bangladesh.

DOI: 10.4236/am.2021.1212079                

Abstract

Numerical diffusion and oscillatory behavior characteristics are averted applying numerical solutions of advection-diffusion equation are themselves immensely sophisticated. In this paper, two numerical methods have been used to solve the advection diffusion equation. We use an explicit finite difference scheme for the advection diffusion equation and semi-discretization on the spatial variable for advection-diffusion equation yields a system of ordinary differential equations solved by Euler’s method. Numerical assessment has been executed with specified initial and boundary conditions, for which the exact solution is known. We compare the solutions of the advection diffusion equation as well as error analysis for both schemes.

Keywords

Advection Diffusion Equation, Finite Difference Scheme, Semi-Discretization, Rate of Convergence, Error Analysis

1. Introduction

The advection diffusion equation (ADE) is the model that can be used for simulation natural processes. Two categories of the advection-diffusion equation: advection is first due to the movement of materials from one region to another; the second category is called diffusion which is due to the movement of materials from higher concentration to low concentration. This mathematical model has a wide range of applications in natural science and engineering. These applications include where simulation techniques are useful for transport of air, river water, adsorption of pollutants in soil, food processing, modeling of the biological system, finance, electromagnetism, fluid mechanics structural dynamics, quantum physical process, etc. The analytical and numerical solutions along with an initial and two boundary conditions help to comprehend pollutant concentration distribution behavior through an open medium like rivers, air, lakes, and porous medium. Various works have been appeared to solve and use this equation in their simulation using finite difference methods [1] [2] [3]. Numerical Solution of the 1D advection-diffusion Equation solved using standard and nonstandard Finite Difference Schemes [4]. The significant application of the linear advection diffusion equation lies in fluid dynamics, heat transfer, and mass transfer [5]. Researchers examine numerical solution of Advection Diffusion Equation using operator splitting method [6]. These methods have been implemented by a characteristic method with cubic spline interpolation (MOC-CS) and Crank-Nicolson (CN) finite difference scheme. Obtained results were compared with analytical solutions. It is seen that the implemented method has lower error than other methods also produces accurate results even when the time steps are great. The linear advection diffusion equation (ADE) is a model which describes the contaminant transport due to the combined effect of advection and diffusion in a porous media [7]. In this study, the advection diffusion equation is solved by explicit finite difference schemes and investigates a different approach, the semi-discretization method: a spatial variable which yields a system of ODE with the temporal independent variable. We solve this system of ODE’s by Euler’s method and develop an algorithm of the Euler method for the system of ODE’s and implement it for the computation of the concentration $u\left(x,t\right)$.

2. Advection Diffusion Equation

We consider the following partial differential equation, which has both an adventive and diffusive terms together.

$u_t\left(x,t\right)+c\left(x,t\right)u_x\left(x,t\right)=D{u}_{xx}\left(x,t\right)$ (1)

with initial condition:

$u\left(x,t_0\right)=f\left(x\right);a<x<b$.

And boundary conditions:

$u\left(a,t\right)=u_a\left(t\right);t_0\ll T$

$u\left(b,t\right)=u_b\left(t\right);t_0<t<T$

where $u_a,u_b$ are concentration values and $u\left(x,t\right)$ is the unknown solution being investigated which indicates concentration, $c\left(x,t\right)$ is the velocity of the medium in the x direction, $D\left(x,t\right)$ is the diffusion coefficient. To introduce numerical scheme for Equation (1), an advection diffusion problem whose general solution [8] [9] [10] is (Figure 1)

$u\left(x,t\right)=\frac{1}{\sqrt{4\pi Dt}}{\text{e}}^{-\frac{{\left(x-ct\right)}^2}{4Dt}}$ (2)

3. Numerical Methodology

In mathematics, our goal is to approximate the solution of the differential equations.

Figure 1. The general solution of the advection diffusion equation.

This gives a large algebraic system of equations to be solved in replace of the differential equation, which can be easily solved [11] [12] on a computer by Matlab code.

3.1. Computational Grid

We consider some simple space discretization on a uniform grid. We divide the spatial interval $\left[0,L\right]$ into $M+1$ equal sub-interval such that

$x_1<x_2<x_3<\cdots <L$ with $x_m=\left(m-1\right)\Delta x,m=1,2,3,\cdots ,M+1$ and $\Delta x=\frac{L}{M}$.

Approximations $u\left(t\right)\approx u{\left(x_m,t\right)}_0$ are found by replacing the spatial derivatives by difference quotients. we also divide the time interval $\left[0,T\right]$ into $N+1$ equal subinterval such that $t_1<t_2<t_3<\cdots <T$ with $t_n=\left(n-1\right)\Delta x,n=1,2,3,\cdots ,N+1$,

and $\Delta t=\frac{T}{N}$. For purpose of the notation $\Delta x=h$ and $\Delta t=k$.

This gives a finite difference discretization in space. Setting

$u\left(t\right)={\left(u_1\left(t\right),\cdots ,u_m\left(t\right)\right)}^{\text{T}}$.

Therefore, we get a system of ordinary difference equations (ODEs) of (1.1)

$u^{\prime }\left(t\right)=F\left(t,u\left(t\right)\right),t>0,u\left(0\right)=u_0$ (3)

with a given initial value $u\left(0\right)$.

3.2. Discretization of Explicit Finite Difference Scheme (ECDS)

To approximate the solution to Equation (1) using the Explicit Centered Difference Scheme, we use the following approximations

$u_t\left(x_m,t^n\right)\cong \frac{u_m^{n+1}-u_m^n}{\Delta t}$ (a)

$u_x\left(x_m,t^n\right)\cong \frac{u_{m+1}^n-u_{m-1}^n}{2\Delta x}$ (b)

$u_{xx}\left(x_m,t^n\right)\cong \frac{u_{m+1}^{n+1}-2u_m^n+u_{m-1}^n}{{\left(\Delta x\right)}^2}$ (c)

where $\Delta x$ is the spatial step, $\Delta t$ is the time step, m and n is spatial and temporal node respectively. Substituting Equation (a), (b), (c) in Equation (1) and solving

for unknown $u_m^{n+1}$. We obtain $u_m^{n+1}=\left(\frac{\alpha }{2}+\gamma \right)u_{m-1}^n+\left(1-2\gamma \right)u_m^n+\left(\gamma -\frac{\alpha }{2}\right)u_{m+1}^n$ where $\frac{c\Delta t}{\Delta x}=\alpha$ and $\frac{D\Delta t}{\Delta x^2}=\gamma$. Stability condition $\frac{c\Delta t}{\Delta x}\le 1$ & $\frac{D\Delta t}{\Delta x^2}\le \frac{1}{2}$.

3.3. Spatial Discretization Technique for ADE and Solved by Euler Method

We consider the following figure for ADE (Advection Diffusion Equation).

We draw vertical grid line as shown in the picture. These lines are parallel to the t-axis and cross the x-axis in $x=x_m$, $m=1,\cdots ,M+1$. $x_m=m\times \Delta x$,

$\Delta x=\frac{1}{M+1}$.

In semi-discretization method (SDM), we assume that the PDE system with its boundary conditions has spatially discretized, and thus we focus on ODE system $u^{\prime }\left(t\right)=F\left(t,u\left(t\right)\right)$, representing semi-discrete advection-diffusion problems.

This equation is same as Equation (2). Now Consider ADE (1) with $c>0$ and $D>0$.

We introduce the function of one variable: $u_m\left(t\right)\approx u\left(t,x_m\right)$, $m=1,\cdots ,M+1$.

Now approximate the first derivative and second derivative ${\partial }_x$ and ${\partial }_{xx}$ respectively as:

$u_x\left(x_m,t\right)\cong \frac{u_{m+1}\left(t\right)-u_{m-1}\left(t\right)}{2\Delta x}$ (d)

$u_{xx}\left(x_m,t\right)\cong \frac{u_{m+1}\left(t\right)-2u_m\left(t\right)+u_{m-1}\left(t\right)}{{\left(\Delta x\right)}^2}$ (e)

This gives the semi-discrete form of (1)

$\begin{array}{l}\frac{\text{d}{u}_m\left(t\right)}{\text{d}t}=-\alpha \frac{\left(u_{m+1}-u_{m-1}\right)}{2\Delta x}+\gamma \frac{\left(u_{m-1}+2u_m+u_{m+1}\right)}{{\left(\Delta x\right)}^2}\\ u_1\left(t\right)=f_0\left(t\right)\\ u_{M+1}\left(t\right)=f_1\left(t\right)\\ m=2,\cdots ,M;\alpha =\frac{c}{2\Delta x}\&\gamma =\frac{D}{{\left(\Delta x\right)}^2}\end{array}$ (4)

Setting $u_m\left(0\right)=f\left(x_m\right),m=2,\cdots ,M$

We get, $\stackrel{\dot{}}{u}\left(t\right):={\left[u_2\left(t\right),\cdots ,u_m\left(t\right)\right]}^{\text{T}}$

$\begin{array}{c}\frac{\text{d}}{\text{d}t}\left(\begin{array}{c}{u}_2\\ ⋮\\ u_M\end{array}\right)=-\alpha \left(\begin{array}{cccccc}0& 1& 0& & & \\ -1& 0& & & & \\ 0& & & & & \\ & & & 0& 1& 0\\ & & & -1& 0& 1\end{array}\right)\left(\begin{array}{c}{u}_2\\ ⋮\\ u_M\end{array}\right)\\ +\gamma \left(\begin{array}{cccccc}-2& 1& 0& & & \\ 1& -2& & & & \\ 0& & & & & \\ & & & -2& 1& 0\\ & & & 1& -2& 1\end{array}\right)\left(\begin{array}{c}{u}_2\\ ⋮\\ u_M\end{array}\right)+\left(\begin{array}{c}-\alpha u_1\\ ⋮\\ -\alpha u_{M+1}\end{array}\right)+\left(\begin{array}{c}\gamma u_1\\ ⋮\\ \gamma u_{M+1}\end{array}\right)\end{array}$

So, we obtain a linear system of ordinary differential equations (ODE’s) of the type.

$\begin{array}{l}\stackrel{\dot{}}{u}\left(t\right)=\frac{\text{d}u}{\text{d}t}=\underset{f\left(t\right)}{\underset{︸}{Au+b\left(t\right)}},u\left(0\right)=u_0={\left[f\left(x_1\right),\cdots ,f\left(x_M\right)\right]}^{\text{T}}\\ \therefore \stackrel{\dot{}}{u}\left(t\right)=F\left(t,u\right),u\left(0\right)=u_0,\text{where}u\text{and}F\left(t,u\right)\text{}\text{are}\text{vector}\end{array}$

Now, we get with time step $\Delta t$ for the numerical solution

$\begin{array}{l}\frac{u^{n+1}-u^n}{\Delta t}=F\left(t_n,u^n\right)\\ u^{n+1}=u^n+\Delta tF\left(t_n,u^n\right)\end{array}$

which is the semi discretization using Euler method of $F\left(t,u\right)=A\stackrel{\dot{}}{u}+b\left(t\right)$.

4. Numerical Result’s and Discussion

4.1. Algorithm for the Semi Discretization of ADE Solved by Euler Method

To approximate the solution to the partial differential equation

$\begin{array}{l}\frac{\partial u}{\partial t}\left(x,t\right)+c\frac{\partial u}{\partial x}\left(x,t\right)-D\frac{{\partial }^{2}u}{\partial x^2}\left(x,t\right)=0,a<x<b\text{and}{t}_0<t<T\\ \text{Subject to the boundary condition}\\ u\left(a,t\right)=u_a\left(t\right),t_0<t<T\\ u\left(b,t\right)=u_b\left(t\right),t_0<t<T\\ \text{And Initial condition,}\\ u\left(x,t_0\right)=f\left(x\right);a<x<b\end{array}$

Input: dt, dx, constant co-efficient $C,D$, $t_0$, the left and right end point $t_f$ of $\left(0,T\right)$ ; $x_d$, the right end point of $\left(0,b\right)$.

Output: approximation $u_{mm}$ to $u\left(x_m,t^n\right)$ for each $mm=2,\cdots ,m$, $n=1,\cdots ,N$.

Step 1: set $\begin{array}{l}nx=\frac{xd-x_0}{dx};nt=\frac{tf-f_0}{dt};\\ \alpha =-\frac{c}{2dx},\lambda =\frac{D}{d{x}^2}\end{array}$

Step 2: for $m=2,\cdots ,nx$ $\begin{array}{l}\text{Set}{u}_{m,1}=f\left(x\right)\left(\text{Initial value}\right)\\ un=u_{m,1};\end{array}$

Step 3: $\begin{array}{l}\text{Set},u_{1,n}=u_a\left(\text{left boundary condition}\right)\\ u_{nx+1,n}=u_b;\left(\text{right boundary condition}\right)\end{array}$

Step 4-Step 6: (solve a tri-diagonal linear system)

Step 4: Construct matrix A and B

Set, $T1=\left(\alpha \ast A\right)+\left(\beta \ast B\right)$

$\begin{array}{l}b\left(1\right)=u\left(1,1\right)\\ b\left(nx-1\right)=u\left(n,1\right)\end{array}\}$ Boundary

$unew=un+dt\left(T1\ast un\right)+b$ (Unknown solution for first time step)

Step 5: For $n=1,\cdots ,nt+1$

Set $u_{mm,n}=unew;un=unew;$

$unew=un+dt\left(T1\ast un\right)+b$ ;

Output $\left(x,u_{m,n}\right)$,

Step 6: Stop (the procedure complete)

4.2. Case Study

Numerical implementation of Euler method [10] [11]: our solving equation (1): $u_t+c{u}_x=D{u}_{xx}$ initial condition: $u_m^0=f\left(x_m\right)=u\left(x,t=0\right)$

Boundary condition: $\begin{array}{l}u\left(x=0,t\right)=u_0^n=f_1\left(t_n\right)\\ u\left(x=b,t\right)=u_b^n=f_b(\; t\; n\; )\end{array}$

We will solve numerically for the concentration u using matrix system of equation.

Suppose we use 4 grid points $x_1,x_2,x_3,x_4=x_{m+1}$ i.e. $m=3$ in this example.

We let, ${\stackrel{\to }{u}}^n=\left(\begin{array}{c}{u}_2^n\\ u_3^n\end{array}\right)$, Solution for concentration vector ${\stackrel{\to }{u}}^n$ at time $t_n$.

The boundary condition gives

$u_1^n=u\left(x=0,t_n\right)$ and $u_{m+1}^n=u_4^n=u\left(x=b,t_b\right)$.

We can rewrite general nth term in Equation (4) to required Euler method of advection diffusion equation.

${\stackrel{\to }{u}}^{n+1}=\left(\begin{array}{c}{u}_2^n\\ u_3^n\end{array}\right)-\alpha \underset{A}{\underset{︸}{\left(\begin{array}{cc}0& 1\\ -1& 0\end{array}\right)}}\left(\begin{array}{c}{u}_2^n\\ u_3^n\end{array}\right)+\lambda \underset{B}{\underset{︸}{\left(\begin{array}{cc}-2& 1\\ 1& -2\end{array}\right)}}\left(\begin{array}{c}{u}_2^n\\ u_3^n\end{array}\right)-\alpha \left(\begin{array}{c}{u}_1^n\\ u_4^n\end{array}\right)+\lambda \left(\begin{array}{c}{u}_1^n\\ u_4^n\end{array}\right)$

$\begin{array}{l}{\stackrel{\to }{u}}^{n+1}=\left(\begin{array}{c}{u}_2^n\\ u_3^n\end{array}\right)+\Delta t\left\{\underset{P}{\underset{︸}{-\alpha A\left(\begin{array}{c}{u}_2^n\\ u_3^n\end{array}\right)+\lambda B\left(\begin{array}{c}{u}_2^n\\ u_3^n\end{array}\right)-\alpha \left(\begin{array}{c}{u}_1^n\\ u_4^n\end{array}\right)+\lambda \left(\begin{array}{c}{u}_1^n\\ u_4^n\end{array}\right)}}\right\}\\ {\stackrel{\to }{u}}^{n+1}=\left(\begin{array}{c}{u}_2^n\\ u_3^n\end{array}\right)+\Delta t\stackrel{\to }{P};\text{}\text{}\alpha =\frac{c}{2\Delta x};\&\lambda =\frac{D}{{\left(\Delta x\right)}^2}\end{array}$

4.3. Semi Discretization Scheme for ADE Solved Using Euler Method

For this case the time step is increased to $\Delta t=0.008$, $\Delta x=0.1$, and the parameter $c=0.2$ and $D=0.02$, $\frac{c\Delta t}{\Delta x}=\alpha$ known as advection equation number and $\frac{D\Delta t}{{\left(\Delta x\right)}^2}=\gamma$ known as diffusion term. Regarding this application, (Figure 2)

$\alpha =\frac{0.2\times 0.008}{0.1}=0.0160\le 1,$

$\gamma =\frac{0.02\times 0.008}{{0.1}^2}=0.0160\le 0.5$

Summary of elapsed time for different temporal grid point present below:

Figure 2. The concentration distribution for $\Delta x=0.1$ at different time of Euler and ECDS.

5. Error Analysis and Convergence

5.1. Error Analysis

The comparison of relative error for two finite difference schemes as a function of time is shown in Figure 3. Two different curves show relative error for Explicit (red) and Euler (blue). We found the relative error for ECDS remains below 0.0009 and the relative error for Euler remains below 0.0011.

5.2. Convergence

Figure 4 illustrates the convergence of relative error by the scheme FTCS techniques. In this figure, numerical computation of ADE is presented by using explicit finite difference methods by FTCS and compared with an exact solution of

Figure 3. Plot shows Relative error of ADE for both schemes.

Figure 4. Plot shows Convergence of relative error for ECDS.

the ADE. A good agreement between the numerical solutions and the analytical solutions is obtained and the error becomes clear when using large size step for the time.

Figure 5 presents the convergence of semi discretization by Euler for ADE with respect to time. It can be seen that, the relative error rate of the convergence curves for the Euler scheme with respect to time. This figure shows a very good rate of convergence.

5.3. Problem Discussion

Figure 6 describes the concentration distribution profile as a function of distance. Different curves represent concentration distribution profiles at different times starting from $t=10\sec$ to $t=50\sec$. The plot marked by “red curved”

Figure 5. Plot shows error rate of convergence for Euler scheme.

Figure 6. Plot shows concentration distribution at different time.

represents the concentration is high for 10 seconds and the green curve shows that concentration is decreased for 20 seconds. The curve marked by “blue” represents concentration is high for 30 seconds and. The plot marked by “yellow” represents concentration is high for 50 seconds we can see that when time is increased concentration profile is decreasing. As can be seen in Figure 6, five curves show concentration distribution profiles for different times, advection, and diffusion parameter as a function of distance.

Similarly, we observe the above figure with the concentration distribution profile for different distances is decreased in a position with respect to time in Figure 7.

Figure 8 represents the concentration profile for different velocities and different diffusion coefficients regarding distances. Different curves show concentration profiles for different applied velocities and diffusion rates. If we vary

Figure 7. Plot shows concentration distribution at different position.

Figure 8. Plot demonstrated varying advection & diffusion rate a time $t=20\sec$.

Figure 9. The plot demonstrated varying diffusion rate a time $t=20\sec$.

both velocity and diffusion coefficient at t = 20 secs, the solutions appeared which can be seen in Figure 8.

Similarly, concentration profile by solving numerically at a different diffusion rate can be seen in Figure 9 where the maximum concentration profile is shown at time $t=20\sec$.

6. Conclusion

We present numerical solutions with exact solutions for the advection-diffusion equation with an initial condition and two boundary conditions by using ECDS and semi discretize method. A numerical study of the ADE has been presented graphically for two different schemes. The numerical solution of ADE by semi-discretization scheme shows a good agreement with the exact solution as well as for ECDS. Though, ECDS seems a more efficient scheme in terms of elapse timing; we compute the relative errors for two different schemes, both schemes show a very good rate of convergence. In comparison to Euler’s method, ECDS shows less error, which is obvious. Semi discretize methods have to deal with a large number of systems of ordinary differential equations in comparison with ECDS. In our next work, we would like to upgrade this work with higher-order accuracy.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Mojtabi, M.D.A. (2015) One-Dimensional Linear Advection–Diffusion Equation: Analytical and Finite Element Solutions. Computers and Fluids, Elsevier, 189-195.
[2]	Ahmed, S. (2012) A Numerical Algorithm for Solving Advection-Diffusion Equation with Constant and Variable Coefficients. The Open Numerical Methods Journal, 4, 1-7.
[3]	Dehghan, M. (2005) On the Numerical Solution of the One Dimensional Convection-Diffusion Equation. Hindawi Publishing Corporation, 61-74.
[4]	Appadu, A.R. (2013) Numerical Solution of the 1D Advection-Diffusion Equation Using Standard and Nonstandard Finite Difference Schemes. Hindawi Publishing Corporation Journal of Applied Mathematics, 2013, Article ID: 734374.
[5]	Bahar, E. and Gurarslan, G. (2017) Numerical Solution of Advection-Diffusion Equation Using Operator Splitting Method. International Journal of Engineering & Applied Sciences, 9, 76-88.
[6]	Azad, T.M.A.K. and Andallah, L.S. (2016) Stability Analysis of Finite Difference Schemes for an Advection Diffusion Equation. Bangladesh Journal of Scientific Research, 29, 143-151.
[7]	Chen-Charpentier, B.M. and Kojouharov, H.V. (2013) Unconditionally Positivity Presentation Scheme for Advection Diffussion Reaction Equation. Mathematical and Computer Modelling, 57, 1277-1285.
[8]	Willem, H. and Amsterdam, C.W.I. (2000) Numerical Solution of Advection-Diffusion-Reaction Equations. Lecture Notes, Thomas Stieltjes Institute, 2000.
[9]	Olsen-Kettle, D.L. (2011) In Numerical Solution of Partial Differential Equations. https://www.researchgate.net/publication/296672444_Numerical_solution_of_partial_differential_equations_and_code
[10]	Socolofsky, S.A. and Jirka, G.H. (2002) Environmental Fluid Mechanics Part I: Mass Transfer and Diffusion, Germany: Universit at Karlsruhe, 76128-Karlsruhe, Institut fur Hydromechanik. 2nd Edition, Germany, 2002.
[11]	Mathews, J.H. and Kurtis, K. (2004) Fink, Book: Numerical Methods Using Matlab. 4th Edition, Prentice Hall, Hoboken.
[12]	Burden, R.L. and Faires, J.D. (2010) Numerical Analysis. 9th Edition, 20 Channel Center Street Boston, 2010.