You are currently viewing a beta version of our website. If you spot anything unusual, kindly let us know.

Preprint

Article

Convergence Analysis of Numerical Solutions of Advection-Diffusion-Reaction Equations Using a Finite Difference Method an Application to Air Pollution Problems

Altmetrics

Downloads

219

Views

Comments

,Sutthisak Phongthanapanich

,Sutthisak Phongthanapanich

This version is not peer-reviewed

Abstract

In this paper, we analyze the convergence of a finite difference method with the implicit forward time central space (FTCS) scheme for solving the three-dimensional advection-diffusion-reaction equation (ADRE). It is proved that the method is unconditionally convergent. We apply the scheme to obtain numerical solutions for the transport of pollutants in street tunnel problems with various reaction coefficients and various rates of change of concentrations of sources or sinks of pollution.

Keywords:

Subject: Computer Science and Mathematics - Computational Mathematics

1. Introduction

Advection-diffusion-reaction equations (ADRE) have been widely used as mathematical models in many areas of science and engineering, for example, in areas such as water pollution, air pollution, molecular diffusion, and chemical engineering. There are now many numerical methods that have been developed to solve linear and nonlinear ADREs. Some examples are as follows. In 2012, Diego et al [1] developed a model consisting of a decoupled system of advection-diffusion-reaction equations, along with the Navier-Stokes equations of incompressible flow, and solved the model by using the finite element method. In 2012, Savovic and Djordjevich [2] proposed an explicit finite difference method to solve a one-dimensional advection-diffusion equation with variable coefficients in semi-infinite media for three dispersion problems. In 2013, Appadu and Gidey [3] proposed two time-splitting procedures that they used to solve a 2D advection-diffusion equation with constant coefficients. In 2013, Bause and Schwegler [4] proposed a higher order finite element approximation for systems of coupled convection-dominated transport equations. In 2015, Mojtabi and Deville [5] proposed a finite element method to solve analytically using separation of variables and numerical solution of a time-dependent one-dimensional linear advection-diffusion equation with Dirichlet homogeneous boundary conditions and an initial sine function. In 2016, Gharehbaghi [6] proposed an explicit and an implicit differential quadrature method to compute a numerical solution of a one-dimensional time-dependent advection-diffusion equation with variable coefficients in a semi-infinite domain. In 2017, Gyrya and Lipnikov [7] proposed an arbitrary order mimetic finite difference discretization for the diffusion equation with non-symmetric tensorial diffusion coefficient in a mixed formulation on general polygonal meshes. In 2017, Bahar and Gurarslan [8] studied the effects of Lie-Trotter and Strang splitting methods on the solution of a one-dimensional advection-diffusion equation. In 2017, Oyjindal and Pochai [9] proposed models for numerical simulations of air pollution measurements near an industrial zone. The numerical experiments consisted of different conditions including normal and controlled emissions. The models were then solved using an explicit finite difference technique and the solutions were compared with the measurements for the controlled and uncontrolled point sources at the monitoring points. In 2018, Suebyat and Pochai [10] developed a three-dimensional air pollution measurement model for a heavy traffic area under a Bangkok sky train platform and used finite difference techniques to approximate the solutions. In 2018, Al-Jawary et al. [11] proposed a semi-analytical technique for finding the exact solutions for different types of Burger’s equations and systems of equations in 1D, 2D, and 3D. Also, the method was applied to solve the diffusion and advection-diffusion equations. In 2018, Kusuma et al. [12] proposed a finite difference (FTCS) method to solve a model of pollution distribution in a street tunnel using two dimensional advection and three dimensional diffusion in a rectangular box domain. In 2018, Lou et al. [13] developed reconstructed Discontinuous Galerkin methods for solving linear advection-diffusion equations on hybrid unstructured grids based on a first-order hyperbolic system formulation. In 2018, Bhatt et al, [14] developed a Krylov subspace approximation-based locally extrapolated exponential time differencing method and studied its accuracy and efficiency for solving three-dimensional nonlinear advection-diffusion-reaction systems. In 2020, Pananu et al. [15] analyzed the convergence of the finite difference method with the implicit forward time central space (FTCS) scheme for the two-dimensional advection-diffusion-reaction equation (ADRE) and applied the scheme to a pollutant dispersion with removal mechanism model in a reservoir. In 2020, Heng and Guodong [16] improved the element-free Galerkin method and used it for solving 3D advection-diffusion problems. In 2020, Cruz-Quintero and Jurado [17] proposed a backstepping design for the boundary control of a reaction-advection-diffusion equation with constant coefficients and Neumann boundary conditions. In 2021, Para et al. [18] proposed the characteristic finite volume method for solving a convection-diffusion problem on two-dimensional triangular grids and compared the accuracy of four piecewise linear reconstruction techniques on structured triangular grids, namely, Frink, Holmes-Connell, Green-Gauss, and least squares methods. In 2021, Hong et al. [19] discussed the numerical solution of 3D unsteady advection-diffusion equations using a meshless numerical scheme with a space-time backward substitution method. In 2021, Irfan and Hidayat [20] proposed a meshless finite difference method with B-splines for the numerical solution of coupled advection-diffusion-reaction problems. In 2021, Shahid et al [21] studied an epidemic type model with advection and diffusion terms for the transmission dynamics of a computer virus model with fixed population density. In 2021, García and Jurado [22] designed an adaptive boundary control for a parabolic type reaction-advection-diffusion PDE under the assumption of unknown parameters for both advection and reaction terms and Robin and Neumann boundary conditions. In 2022, Para et al. [23] proposed an explicit characteristic-based finite volume method (FVM) for the numerical solution of advection-diffusion-reaction equations (ADRE) and applied it to solve some 1D and 2D water pollution problems which can be modeled in terms of ADREs. Then, the FVM results were compared with numerical results obtained using a finite difference method (FDM) with an implicit forward time central space (FTCS) scheme.

In this article, we consider the application of a finite difference method based on the implicit forward time central space (FTCS) scheme to the numerical solution of an advection-diffusion-reaction equation of the following form:

\begin{matrix} \frac{\partial ϕ}{\partial t} + \nabla \cdot (v ϕ - ε \nabla ϕ) + κ ϕ = q, \end{matrix}

(1)

where

ϕ

is a scalar quantity,

v = v (x)

is a given advection velocity vector,

ε \geq 0

is a diffusion coefficient,

κ

is a reaction coefficient,

q = q (x, t)

is a prescribed source term,

x

is a position vector and time

t \in [0, T]

, where

T > 0

. The present paper is organized as follows. The convergence theory is briefly provided in Section 2. The consistency and stability analysis and numerical results of the method are given in Section 3. Using the method, numerical results and their graphs for the air pollution problems written in terms of the ADRE in Eq. (1) are shown in Section 4. The conclusion of our work is presented in Section 5.

2. Convergence Theory

We consider the following initial time, space boundary value problem (IBVP) for the advection-diffusion-reaction equation

L ϕ (x, t) = \frac{\partial ϕ (x, t)}{\partial t} + \nabla \cdot (v (x) ϕ (x, t) - ε \nabla ϕ (x, t)) + κ ϕ (x, t) = q (x, t),

(2)

where

(x, t) \in Ω \times (0, T]

, and with the initial time condition

ϕ (x, 0) = ψ (x) for x \in Ω,

(3)

and either Dirichlet or Neumann boundary space conditions

\begin{matrix} Dirichlet & ϕ (x, t) = ξ_{D} (x, t) \\ Neumann & ε \frac{\partial ϕ (x, t)}{\partial n} = ξ_{N} (x, t) for (x, t) \in \partial Ω \times [0, T] . \end{matrix}

(4)

In Eqs. (2)–(4), T is a finite time,

\partial Ω

is the boundary of the spatial domain

Ω

ξ_{D}

and

ξ_{N}

are given functions,

v (x)

is an advection vector,

ε > 0

is a diffusion coefficient,

κ

is a reaction constant and

q (x, t)

is a source term. In addition, L is a linear partial differential operator from a space of continuous functions F to a space of continuous functions H, where the dependent variable

ϕ (x, t) \in F

, and where

L ϕ (x, t) \in H

and

q (x, t) \in H

The first step in developing the finite difference method is to discretize

Ω

\partial Ω

and

[0, T]

. We assume that there are

N_{1}

step sizes of length

Δ x = \frac{L_{1}}{N_{1}}

in the x direction,

N_{2}

step sizes of length

Δ y = \frac{L_{2}}{N_{2}}

in the y direction and

N_{3}

step sizes of length

Δ z = \frac{L_{3}}{N_{3}}

in the z direction. For simplicity, we will assume that

Δ x = Δ y = Δ z = h

and therefore

L_{1} = N_{1} h

L_{2} = N_{2} h

and

L_{3} = N_{3} h

. Then the discretized version of

Ω

is the spatial grid

Ω_{h} = {(x_{i}, y_{j}, z_{k}) | x_{i} = i h, y_{j} = j h, z_{k} = k h, i \in J_{1}, j \in J_{2}, k \in J_{3}},

(5)

where

J_{1} = {0, 1, 2, . . ., N_{1}}

J_{2} = {0, 1, 2, . . ., N_{2}}

J_{3} = {0, 1, 2, . . ., N_{3}}

. The boundary of the spatial grid is

\begin{matrix} \partial Ω_{h} & = & {x_{0} = 0; y_{j} = j h, j \in J_{2}, z_{k} = k h, k \in J_{3}} \\ \cup {x_{N_{1}} = L_{1}; y_{j} = j h, j \in J_{2}, z_{k} = k h, k \in J_{3}} \\ \cup {x_{i} = i h, i \in J_{1}; y_{0} = 0; z_{k} = k h, k \in J_{3}} \\ \cup {x_{i} = i h, i \in J_{1}; y_{N_{2}} = L_{2}; z_{k} = k h, k \in J_{3}} \\ \cup {x_{i} = i h, i \in J_{1}, y_{j} = j h, j \in J_{2}; z_{0} = 0} \\ \cup {x_{i} = i h, i \in J_{1}, y_{j} = j h, j \in J_{2}; z_{N_{3}} = L_{3}} . \end{matrix}

(6)

Similarly, the interval

[0, T]

is discretized with N step sizes of size

Δ t = \frac{T}{N}

. The set of discretized times is then

T_{K} = {t^{n} | t^{n} = n Δ t, n \in K}, K = {0, 1, 2, \dots, N} .

(7)

Finally, we define

F_{h}

and

H_{h}

as discretized function spaces with domain

Ω_{h} \times T_{K}

. We also define a projection

P_{h} (ϕ) \in F_{h}

of the exact solution

ϕ (x, t)

of (2) onto

Ω_{h} \times T_{K}

and a projection

P_{h} (q) \in H_{h}

of the source term

q (x, t)

onto

Ω_{h} \times T_{K}

We can then write a finite difference scheme (FDS) corresponding to the continuous partial differential equation (2) in the form

L_{h} ϕ^{h} (x, t) = q^{h} (x, t), (x, t) \in Ω_{h} \times T_{K},

(8)

where

L_{h} : F_{h} \to H_{h}

is a finite difference operator, corresponding to the partial differential operator L, which acts from the discrete function space

F_{h}

to the discrete function space

H_{h}

. The initial condition in (3) can then be rewritten as the discretized initial condition

ϕ^{h} (x, 0) = ψ^{h} (x), x \in Ω_{h},

(9)

and, for the Dirichlet case, the boundary conditions in (4) can be rewritten as the discretized boundary conditions

ϕ^{h} (x, t) = ξ^{h} (x, t) for (x, t) \in \partial Ω_{h} \times T_{K} .

(10)

Lax’s equivalence theorem [24]

An initial time, space boundary value problem (IBVP) is called a properly posed IBVP if it has a unique solution. In this paper, all advection-diffusion-reaction equations that we study will be properly posed IBVP.

Then Lax’s equivalence theorem states that necessary and sufficient conditions for the solution of a finite difference approximation to converge to the unique solution of the IBVP are the following conditions of convergence, consistency and stability.

Using the definition given above of the projection

P_{h}

of a function in the continuous space H onto the discrete space

H_{h}

, we define the truncation error (T.E.) between the finite difference system Eq. (8) and the IBVP Eq. (2) as follows

T . E . = P_{h} (L ϕ (x, t)) - L_{h} ϕ^{h} (x, t), (x, t) \in Ω_{h} \times T_{K} .

(11)

Introducing the norm in the set of discrete functions as the infinity norm

\begin{matrix} \begin{matrix} ∥ ϕ^{h} ∥_{F_{h}} & = \max_{\begin{matrix} i, j . k \in J \\ n \in K \end{matrix}} | ϕ^{h} (x_{i}, y_{j}, z_{k}, t^{n}) |, \end{matrix} \end{matrix}

(12)

we have the definition of convergence for numerical methods as follows

Definition 1.

Convergence

A solution

ϕ^{h}

of the FDS in Eq. (8)-(10) converges to the solution ϕ of the IBVP in Eq. (2)-(4) if

∥ P_{h} (ϕ) - ϕ^{h} ∥_{F_{h}} \to 0 a s h \to 0 .

(13)

Definition 2.

Convergence with order m

The FDS in Eq. (8) converges with order m if

∥ P_{h} (ϕ) - ϕ^{h} ∥_{F_{h}} \leq C h^{m},

(14)

where C is a positive constant that does not depend on h.

1) Consistency : A finite difference representation of a PDE is said to be consistent [24] if the difference between the PDE and its difference representation vanishes as the mesh is refined, i.e.,

Definition 3.

The FDS (8) is a consistent approximation of the IBVP (2)-(4) with order m if

{∥ T . E . ∥}_{H_{h}} \leq C h^{m},

(15)

where C is a positive constant that does not depend on h.

2) Stability : The finite difference scheme defined by (8) with linear operator

L_{h}

will be called stable, if there exists

h_{0} > 0

such that for arbitrary

h < h_{0}

and for

∥ f^{h} ∥_{H_{h}} = \max_{\begin{matrix} i, j, k \in J \end{matrix}} | ψ (x_{i}, y_{j}, z_{k}, t^{n}) | + \max_{\begin{matrix} i, j, k \in J \\ n \in K - {0} \end{matrix}} | q (x_{i}, y_{j}, z_{k}, t^{n}) |,

(16)

the solution

ϕ^{h}

of the FDS,

L_{h} ϕ^{h} = q^{h}

in Eq. (8), exists and is unique and satisfies the inequality

∥ ϕ^{h} ∥_{F_{h}} \leq C {∥ f^{h} ∥}_{H_{h}},

(17)

where the positive constant C does not depend on h.

The inequality (17) has to be true for any initial conditions

ψ (x_{i}, y_{j}, z_{,} 0)

and source terms

q (x_{i}, y_{j}, z_{k}, t^{n})

of Eq. (2). In particular, if

q (x_{i}, y_{j}, z_{k}, t^{n}) = 0

, then the condition (17) reduces to the necessary condition for stability of the homogeneous version of (8).

For the special case of Fourier or Von Neumann Analysis [24], a solution of the FDS (8) for

q (x_{i}, y_{j}, z_{k}, t_{n}) = 0

can be written in the form

ϕ^{h} (x_{i}, y_{j}, z_{k}, t^{n}) = λ^{n} e^{I (α i + β j + η k)}, (i, j, k) \in J, n \in K, I = \sqrt{- 1},

(18)

where

α

β

and

η

are wave numbers and

e^{I (α i + β j + η k)}

are eigenfunctions corresponding to eigenvalues

λ

L_{h}

. The necessary condition for stability of FDS (8) for

q (x_{i}, y_{j}, z_{k}, t^{n}) = 0

will then hold for all

α, β, η \in ℜ

if the following inequality holds:

| λ | \leq 1 .

(19)

3. Consistency and Stability Analysis

Consider the following initial-boundary value problem

\begin{matrix} \frac{\partial ϕ}{\partial t} + u \frac{\partial ϕ}{\partial x} + v \frac{\partial ϕ}{\partial y} + w \frac{\partial ϕ}{\partial z} - ε (\frac{\partial^{2} ϕ}{\partial x^{2}} + \frac{\partial^{2} ϕ}{\partial y^{2}} + \frac{\partial^{2} ϕ}{\partial z^{2}}) + κ ϕ = q (x, y, z, t), \end{matrix}

\begin{matrix} (x, y, z) \in (0, L_{1}) \times (0, L_{2}) \times (0, L_{3}), t \in (0, T], \\ subject to & ϕ (x, y, z, 0) = ψ (x, y, z), \\ ϕ (0, y, z, t) = α_{1}, ϕ (L_{1}, y, z, t) = α_{2}, \\ ϕ (x, 0, z, t) = β_{1}, ϕ (x, L_{2}, z, t) = β_{2}, \end{matrix}

(20)

\begin{matrix} ϕ (x, y, 0, t) = η_{1}, ϕ (x, y, L_{3}, t) = η_{2}, \end{matrix}

(21)

in which the above partial differential equation is the component form of (1). Here

ϕ = ϕ (x, y, z, t)

is the dependent variable and

v = (u, v, w)

is the advection velocity where

u, v

and w are the velocities in the x-, y- and z-directions, respectively.

In this section, we will analyze the convergence of the finite difference method with the implicit FTCS scheme applied to (20). We discretize the domain of the problem:

(x_{i}, y_{j}, z_{k}, t^{n}) \in Ω_{h} \times T_{K}

, where the spatial grid

Ω_{h}

is defined in (5) and the time grid

T_{K} = [0, T]

has points

t^{n} = n Δ t

with

n \in K = {0, 1, 2, . . ., N} .

Following the method in Section 2, we use Lax’s equivalence theorem to prove the convergence of the FTCS by proving its consistency and stability. We first show the consistency of the FTCS. The implicit FTCS scheme of Eq. (20) is as follows.

\begin{matrix} L_{h} ϕ^{h} = & \frac{ϕ_{i, j, k}^{n + 1} - ϕ_{i, j, k}^{n}}{Δ t} + u \frac{ϕ_{i + 1, j, k}^{n + 1} - ϕ_{i - 1, j, k}^{n + 1}}{2 Δ x} + v \frac{ϕ_{i, j + 1, k}^{n + 1} - ϕ_{i, j - 1, k}^{n + 1}}{2 Δ y} + w \frac{ϕ_{i, j, k + 1}^{n + 1} - ϕ_{i, j, k - 1}^{n + 1}}{2 Δ z} \\ - ε [\frac{ϕ_{i + 1, j, k}^{n + 1} - 2 ϕ_{i, j, k}^{n + 1} + ϕ_{i - 1, j, k}^{n + 1}}{{(Δ x)}^{2}} + \frac{ϕ_{i, j + 1, k}^{n + 1} - 2 ϕ_{i, j, k}^{n + 1} + ϕ_{i, j - 1, k}^{n + 1}}{{(Δ y)}^{2}} \end{matrix}

(22)

\begin{matrix} + \frac{ϕ_{i, j, k + 1}^{n + 1} - 2 ϕ_{i, j, k}^{n + 1} + ϕ_{i, j, k - 1}^{n + 1}}{{(Δ z)}^{2}}] + κ ϕ_{i, j, k}^{n + 1} = q_{i, j, k}^{n}, \\ subject to & ϕ (x_{i}, y_{j}, z_{k}, 0) = ψ (x_{i}, y_{j}, z_{k}), \\ ϕ (0, y_{j}, z_{k}, t^{n}) = α_{1}, ϕ (L_{1}, y_{j}, z_{k}, t^{n}) = α_{2}, \\ ϕ (x_{i}, 0, z_{k}, t^{n}) = β_{1}, ϕ (x_{i}, L_{2}, z_{k}, t^{n}) = β_{2}, \end{matrix}

(23)

\begin{matrix} ϕ (x_{i}, y_{j}, 0, t^{n}) = η_{1}, ϕ (x_{i}, y_{j}, L_{3}, t^{n}) = η_{2}, \end{matrix}

where

ϕ_{i, j, k}^{n} = ϕ^{h} (x_{i}, y_{j}, z_{k}, t^{n})

and

q_{i, j, k}^{n} = q (x_{i}, y_{j}, z_{k}, t^{n})

The truncation error for (22) obtained using the method is

\begin{matrix} T . E . & = & - \frac{Δ t}{2!} \frac{\partial^{2} ϕ}{\partial t^{2}} - u \frac{{(Δ x)}^{2}}{3!} \frac{\partial^{3} ϕ}{\partial x^{3}} - v \frac{{(Δ y)}^{2}}{3!} \frac{\partial^{3} ϕ}{\partial y^{3}} - w \frac{{(Δ z)}^{2}}{3!} \frac{\partial^{3} ϕ}{\partial z^{3}} \\ + ε (\frac{2 {(Δ x)}^{2}}{4!} \frac{\partial^{4} ϕ}{\partial x^{4}} + \frac{2 {(Δ y)}^{2}}{4!} \frac{\partial^{4} ϕ}{\partial y^{4}} + \frac{2 {(Δ z)}^{2}}{4!} \frac{\partial^{4} ϕ}{\partial z^{4}}) \\ + O ({(Δ t)}^{2}, {(Δ x)}^{4}, {(Δ y)}^{4}, {(Δ z)}^{4}) . \end{matrix}

(24)

Therefore, from (24), we have

\begin{matrix} T . E . = O ((Δ t), {(Δ x)}^{2}, {(Δ y)}^{2}, {(Δ z)}^{2}) . \end{matrix}

Taking

(Δ t, Δ x, Δ y, Δ z) \to 0

, then we have

\begin{matrix} lim_{(Δ t, Δ x, Δ y, Δ z) \to 0} ∥ T . E . ∥ = lim_{(Δ t, Δ x, Δ y, Δ z) \to 0} ∥ O ((Δ t), {(Δ x)}^{2}, {(Δ y)}^{2}, {(Δ z)}^{2}) ∥ = 0 . \end{matrix}

We will now prove the stability of the finite difference method with the implicit FTCS scheme for numerical solution of (20). In order to illustrate stability of the method, we initially show that condition (19) holds for the homogeneous equation of (22) and then that condition (17) is satisfied for the nonhomogeneous equation (22). For convenience, we set

τ = Δ t, h = Δ x = Δ y = Δ z

Case I: Homogeneous equation

Setting

q_{i, j, k}^{n} = 0

in equation (22), we obtain the discretized homogeneous version of (22). Assuming a discretized solution of the resulting equation as

ϕ_{i, j, k}^{n} = λ^{n} e^{I (α i + β j + η k)}

I = \sqrt{- 1}

[24], we obtain the following equation for the eigenvalue

λ

\begin{matrix} \frac{λ^{n + 1} e^{I (α i + β j + η k)} - λ^{n} e^{I (α i + β j + η k)}}{τ} + u \frac{λ^{n + 1} e^{I (α (i + 1) + β j + η k)} - λ^{n + 1} e^{I (α (i - 1) + β j + η k)}}{2 h} \\ + v \frac{λ^{n + 1} e^{I (α i + β (j + 1) + η k)} - λ^{n + 1} e^{I (α i + β (j - 1) + η k)}}{2 h} \\ + w \frac{λ^{n + 1} e^{I (α i + β j + η (k + 1))} - λ^{n + 1} e^{I (α i + β j + η (k - 1))}}{2 h} \end{matrix}

\begin{matrix} - ε [\frac{λ^{n + 1} e^{I (α (i + 1) + β j + η k)} - 2 λ^{n + 1} e^{I (α i + β j + η k)} + λ^{n + 1} e^{I (α (i - 1) + β j + η k)}}{h^{2}} \\ + \frac{λ^{n + 1} e^{I (α i + β (j + 1) + η k)} - 2 λ^{n + 1} e^{I (α i + β j + η k)} + λ^{n + 1} e^{I (α i + β (j - 1) + η k)}}{h^{2}} \\ + \frac{λ^{n + 1} e^{I (α i + β j + η (k + 1))} - 2 λ^{n + 1} e^{I (α i + β j + η k)} + λ^{n + 1} e^{I (α i + β j + η (k - 1))}}{h^{2}}] \\ + κ λ^{n + 1} e^{I (α i + β j + η k)} = 0 . \end{matrix}

(25)

Then, we can solve (25) for the eigenvalue

λ

as follows.

\begin{matrix} λ^{n} e^{I (α i + β j + η k)} [\frac{λ - 1}{τ}] - λ^{n + 1} e^{I (α i + β j + η k)} [- u \frac{e^{I α} - e^{- I α}}{2 h} - v \frac{e^{I β} - e^{- I β}}{2 h} \\ - w \frac{e^{I η} - e^{- I η}}{2 h} + ε \frac{e^{I α} - 2 + e^{- I α}}{h^{2}} + ε \frac{e^{I β} - 2 + e^{- I β}}{h^{2}} + ε \frac{e^{I η} - 2 + e^{- I η}}{h^{2}} - κ] = 0, \\ λ = \frac{1}{1 + \frac{τ}{h} [\frac{4 ε}{h} ({sin}^{2} \frac{α}{2} + {sin}^{2} \frac{β}{2} + {sin}^{2} \frac{η}{2}) + I (u sin α + v sin β + w sin η) + κ h]} . \end{matrix}

Letting

γ = \frac{4 ε τ}{h^{2}}

and

μ = \frac{τ}{h}

, we then obtain,

\begin{matrix} λ = \frac{1}{1 + γ ({sin}^{2} \frac{α}{2} + {sin}^{2} \frac{β}{2} + {sin}^{2} \frac{η}{2}) + μ I (u sin α + v sin β + w sin η) + κ τ} . \end{matrix}

Consequently, the magnitude of

λ

\begin{matrix} | λ | & = & \frac{1}{\sqrt{{[1 + γ ({sin}^{2} \frac{α}{2} + {sin}^{2} \frac{β}{2} + {sin}^{2} \frac{η}{2}) + κ τ]}^{2} + {[μ (u sin α + v sin β + w sin η)]}^{2}}}, \\ \leq & 1 . \end{matrix}

Therefore, condition (19) for Von Neumann stability, which is independent of

τ

and h, has been proved.

Case II: Nonhomogeneous equation

We now prove that condition (17) is satisfied for the nonhomogeneous equation (22) with

q_{i, j, k}^{n} \neq 0

by finding a solution

ϕ_{i, j, k}^{n + 1}

for the linear system of difference equations (22) and (23) if

ϕ_{i, j, k}^{n}

is assumed known. Substituting

Δ x = Δ y = Δ z = h

and

Δ t = τ

in (22) and (23), we obtain

\begin{matrix} L_{h} ϕ^{h} & = & \frac{ϕ_{i, j, k}^{n + 1} - ϕ_{i, j, k}^{n}}{τ} + u \frac{ϕ_{i + 1, j, k}^{n + 1} - ϕ_{i - 1, j, k}^{n + 1}}{2 h} + v \frac{ϕ_{i, j + 1, k}^{n + 1} - ϕ_{i, j - 1, k}^{n + 1}}{2 h} + w \frac{ϕ_{i, j, k + 1}^{n + 1} - ϕ_{i, j, k - 1}^{n + 1}}{2 h} \\ - ε [\frac{ϕ_{i + 1, j, k}^{n + 1} - 2 ϕ_{i, j, k}^{n + 1} + ϕ_{i - 1, j, k}^{n + 1}}{h^{2}} + \frac{ϕ_{i, j + 1, k}^{n + 1} - 2 ϕ_{i, j, k}^{n + 1} + ϕ_{i, j - 1, k}^{n + 1}}{h^{2}} \end{matrix}

\begin{matrix} + \frac{ϕ_{i, j, k + 1}^{n + 1} - 2 ϕ_{i, j, k}^{n + 1} + ϕ_{i, j, k - 1}^{n + 1}}{h^{2}} + \frac{ϕ_{i, j, k + 1}^{n + 1} - 2 ϕ_{i, j, k}^{n + 1} + ϕ_{i, j, k - 1}^{n + 1}}{h^{2}}] + κ ϕ_{i, j, k}^{n + 1} = q_{i, j, k}^{n} \\ subject to ϕ_{0, j, k}^{n} = α_{1}, ϕ_{N_{1}, j, k}^{n} = α_{2}; j = 0, 1, . . ., N_{2}; k = 0, 1, . . ., N_{3} \\ ϕ_{i, 0, k}^{n} = β_{1}, ϕ_{i, N_{2}, k}^{n} = β_{2}; i = 0, 1, . . ., N_{1}; k = 0, 1, . . ., N_{3} \end{matrix}

(26)

\begin{matrix} ϕ_{i, j, 0}^{n} = η_{1}, ϕ_{i, j, N_{3}}^{n} = η_{2}; i = 0, 1, . . ., N_{1}; j = 0, 1, . . ., N_{2} \end{matrix}

(27)

After rearranging, we can rewrite (26) in the form

a ϕ_{i - 1, j, k}^{n + 1} + b ϕ_{i, j - 1, k}^{n + 1} + c ϕ_{i, j, k - 1}^{n + 1} + d ϕ_{i, j, k}^{n + 1} + e ϕ_{i + 1, j, k}^{n + 1} + f ϕ_{i, j + 1, k}^{n + 1} + g ϕ_{i, j, k + 1}^{n + 1} = ϕ_{i, j, k}^{n},

(28)

\begin{matrix} where & a = (\frac{u τ}{2 h} + \frac{ε τ}{h^{2}}), b = (\frac{v τ}{2 h} + \frac{ε τ}{h^{2}}), c = (\frac{w τ}{2 h} + \frac{ε τ}{h^{2}}), \\ d = (- 1 - \frac{6 ε τ}{h^{2}} - κ τ), e = (- \frac{u τ}{2 h} + \frac{ε τ}{h^{2}}), f = (- \frac{v τ}{2 h} + \frac{ε τ}{h^{2}}), \\ g = (- \frac{w τ}{2 h} + \frac{ε τ}{h^{2}}), ϕ_{i, j, k}^{n} = - ϕ_{i, j, k}^{n} - τ q_{i, j, k}^{n}, \end{matrix}

and where the step sizes

τ

and h can be chosen so that the coefficients

a, b, c, d, e, f, g

satisfy the conditions

| d | > | a | + | b | + | c | + | e | + | f | + | g | + δ, δ > 0 .

(29)

Lemma 1.

If the coefficients of Eq.(28) satisfy the conditions Eq.(29) then the solution of Eq.(28) exists and is unique and satisfies the inequality

| ϕ_{i, j, k}^{n + 1} | \leq \max \{| α_{1} |, | α_{2} |, | β_{1} |, | β_{2} |, | η_{1} |, | η_{2} |, \frac{1}{δ} \max_{p, q, r} | ϕ_{p, q, r}^{n} |\} .

(30)

Proof

First, we prove inequality (30). Assume that

| ϕ_{p, q, r}^{n + 1} | = \max {| ϕ_{i, j, k}^{n + 1} |, i = 0, 1, . . ., N_{1}, j = 0, 1, . . ., N_{2}, k = 0, 1, . . ., N_{3}}

Let

0 < p < N_{1}, 0 < q < N_{2}, 0 < r < N_{3}

. Then

\begin{matrix} | d | | ϕ_{p, q, r}^{n + 1} | & = & | - a ϕ_{p - 1, q, r}^{n + 1} - b ϕ_{p, q - 1, r}^{n + 1} - c ϕ_{p, q, r - 1}^{n + 1} - e ϕ_{p + 1, q, r}^{n + 1} - f ϕ_{p, q + 1, r}^{n + 1} - g ϕ_{p, q, r + 1}^{n + 1} + ϕ_{p, q, r}^{n} | \\ \leq & | a ϕ_{p - 1, q, r}^{n + 1} | + | b ϕ_{p, q - 1, r}^{n + 1} | + | c ϕ_{p, q, r - 1}^{n + 1} | + | e ϕ_{p + 1, q, r}^{n + 1} | + | f ϕ_{p, q + 1, r}^{n + 1} | + | g ϕ_{p, q, r + 1}^{n + 1} | \\ + | ϕ_{p, q, r}^{n} | \\ = & | a | | ϕ_{p - 1, q, r}^{n + 1} | + | b | | ϕ_{p, q - 1, r}^{n + 1} | + | c | | ϕ_{p, q, r - 1}^{n + 1} | + | e | | ϕ_{p + 1, q, r}^{n + 1} | + | f | | ϕ_{p, q + 1, r}^{n + 1} | \\ + | g | | ϕ_{p, q, r + 1}^{n + 1} | + | ϕ_{p, q, r}^{n} | \\ \leq & (| a | + | b | + | c | + | e | + | f | + | g |) | ϕ_{p, q, r}^{n + 1} | + | ϕ_{p, q, r}^{n} |, \\ | ϕ_{p, q, r}^{n + 1} | & \leq & \frac{| ϕ_{p, q, r}^{n} |}{| d | - | a | - | b | - | c | - | e | - | f | - | g |} \\ \leq & \frac{| ϕ_{p, q, r}^{n} |}{δ} . \end{matrix}

This completes the proof.

Now we will prove the stability of the implicit FTCS scheme by showing that inequality (17) is satisfied.

Proof

Let

i^{*}

j^{*}

and

k^{*}

be the smallest non-negative integers such that

| ϕ_{i^{*}, j^{*}, k^{*}}^{n + 1} | = \max_{\begin{matrix} (i, j, k) \in J \end{matrix}} | ϕ_{i, j, k}^{n + 1} |,

(31)

where

J = {(i, j, k) | i = 0, 1, 2, . . ., N_{1}, j = 0, 1, 2, . . ., N_{2}, k = 0, 1, 2, . . ., N_{3}} .

Now let

\bar{J}

be the discretized boundary set of J, i.e.,

\begin{matrix} \bar{J} & = & {(i, j, k) \in J | i = 0, N_{1}} \cup {(i, j, k) \in J | i \neq 0, i \neq N_{1}, j = 0, N_{2}} \\ \cup {(i, j, k) \in J | i \neq 0, i \neq N_{1}, j \neq 0, j \neq N_{2}, k = 0, N_{3}} . \end{matrix}

(32)

It is obvious that if

(i^{*}, j^{*}, k^{*}) \in \bar{J},

then we have

\max_{(i, j, k) \in J} | ϕ_{i, j, k}^{n + 1} | \leq \max {| α_{1} |, | α_{2} |, | β_{1} |, | β_{2} |, | η_{1} |, | η_{2} |} .

(33)

(i^{*}, j^{*}, k^{*}) \in J = J - \bar{J},

then replacing

(i, j, k)

in (26) with

(i^{*}, j^{*}, k^{*})

we obtain

\begin{matrix} - ϕ_{i^{*}, j^{*}, k^{*}}^{n} - τ q_{i^{*}, j^{*}, k^{*}}^{n} = (\frac{u τ}{2 h} + \frac{ε τ}{h^{2}}) ϕ_{i^{*} - 1, j^{*}, k^{*}}^{n + 1} + (\frac{v τ}{2 h} + \frac{ε τ}{h^{2}}) ϕ_{i^{*}, j^{*} - 1, k^{*}}^{n + 1} + (\frac{w τ}{2 h} + \frac{ε τ}{h^{2}}) ϕ_{i^{*}, j^{*}, k^{*} - 1}^{n + 1} \\ - (1 + \frac{6 ε τ}{h^{2}} + κ τ) ϕ_{i^{*}, j^{*}, k^{*}}^{n + 1} - (\frac{u τ}{2 h} - \frac{ε τ}{h^{2}}) ϕ_{i^{*} + 1, j^{*}, k^{*}}^{n + 1} \\ - (\frac{v τ}{2 h} - \frac{ε τ}{h^{2}}) ϕ_{i^{*}, j^{*} + 1, k^{*}}^{n + 1} - (\frac{w τ}{2 h} - \frac{ε τ}{h^{2}}) ϕ_{i^{*}, j^{*}, k^{*} + 1}^{n + 1} . \end{matrix}

(34)

Next we make the following hypothesis:

\begin{matrix} (H 1) : \begin{matrix} u [ϕ_{i^{*} - 1, j^{*}, k^{*}}^{n + 1} - ϕ_{i^{*} + 1, j^{*}, k^{*}}^{n + 1}] \leq 0, \\ v [ϕ_{i^{*}, j^{*} - 1, k^{*}}^{n + 1} - ϕ_{i^{*}, j^{*} + 1, k^{*}}^{n + 1}] \leq 0, \\ w [ϕ_{i^{*}, j^{*}, k^{*} - 1}^{n + 1} - ϕ_{i^{*}, j^{*}, k^{*} + 1}^{n + 1}] \leq 0 . \end{matrix} \end{matrix}

(35)

Without loss of generality we assume that

ϕ_{i^{*}, j^{*}, k^{*}}^{n + 1} > 0

and further assume that the hypothesis (H1) holds. We can estimate the right side of (34) as

\begin{matrix} \frac{ε τ}{h^{2}} [ϕ_{i^{*} - 1, j^{*}, k^{*}}^{n + 1} - ϕ_{i^{*}, j^{*}, k^{*}}^{n + 1}] + \frac{ε τ}{h^{2}} [ϕ_{i^{*}, j^{*} - 1, k^{*}}^{n + 1} - ϕ_{i^{*}, j^{*}, k^{*}}^{n + 1}] + \frac{ε τ}{h^{2}} [ϕ_{i^{*}, j^{*}, k^{*} - 1}^{n + 1} - ϕ_{i^{*}, j^{*}, k^{*}}^{n + 1}] \\ + \frac{ε τ}{h^{2}} [ϕ_{i^{*} + 1, j^{*}, k^{*}}^{n + 1} - ϕ_{i^{*}, j^{*}, k^{*}}^{n + 1}] + \frac{ε τ}{h^{2}} [ϕ_{i^{*}, j^{*} + 1}^{n + 1} - ϕ_{i^{*}, j^{*}, k^{*}}^{n + 1}] + \frac{ε τ}{h^{2}} [ϕ_{i^{*}, j^{*}, k^{*} + 1}^{n + 1} - ϕ_{i^{*}, j^{*}, k^{*}}^{n + 1}] \\ + \frac{u τ}{2 h} [ϕ_{i^{*} - 1, j^{*}, k^{*}}^{n + 1} - ϕ_{i^{*} + 1, j^{*}, k^{*}}^{n + 1}] + \frac{v τ}{2 h} [ϕ_{i^{*}, j^{*} - 1, k^{*}}^{n + 1} - ϕ_{i^{*}, j^{*} + 1, k^{*}}^{n + 1}] \\ + \frac{w τ}{2 h} [ϕ_{i^{*}, j^{*}, k^{*} - 1}^{n + 1} - ϕ_{i^{*}, j^{*}, k^{*} + 1}^{n + 1}] - ϕ_{i^{*}, j^{*}, k^{*}}^{n + 1} - κ τ ϕ_{i^{*}, j^{*}, k^{*}}^{n + 1} \\ \leq - ϕ_{i^{*}, j^{*}, k^{*}}^{n + 1} . \end{matrix}

(36)

Then, substituting (36) into the right hand side of (34), we obtain

ϕ_{i^{*}, j^{*}, k^{*}}^{n + 1} \leq ϕ_{i^{*}, j^{*}, k^{*}}^{n} + τ q_{i^{*}, j^{*}, k^{*}}^{n} .

(37)

Therefore,

\begin{matrix} \max_{\begin{matrix} (i, j, k) \in J \end{matrix}} | ϕ_{i, j, k}^{n + 1} | \leq \max_{\begin{matrix} (i, j, k) \in J \end{matrix}} | ϕ_{i, j, k}^{n} + τ q_{i, j, k}^{n} | \leq \max_{\begin{matrix} (i, j, k) \in J \end{matrix}} | ϕ_{i, j, k}^{n} | + τ \max_{\begin{matrix} (i, j, k) \in J \\ n \in K \end{matrix}} | q_{i, j, k}^{n} |, \end{matrix}

(38)

where

K = {0, 1, 2, . . ., N} .

We now establish the following inequality representing the maximum principle as

\begin{matrix} \max_{\begin{matrix} (i, j, k) \in J \end{matrix}} | ϕ_{i, j, k}^{n + 1} | \leq \max \{| α_{1} |, | α_{2} |, | β_{1} |, | β_{2} |, | η_{1} |, | η_{2} |, \max_{\begin{matrix} (i, j, k) \in J \end{matrix}} | ϕ_{i, j, k}^{n} | + τ \max_{\begin{matrix} (i, j, k) \in J \\ n \in K \end{matrix}} | q_{i, j, k}^{n} |\} . \end{matrix}

(39)

We will now prove that for the given initial and boundary conditions in (27), the maximum solution

\max_{(i, j, k) \in J} | ϕ_{i, j, k}^{n} |

is bounded for all

n \leq N

for finite

t^{N} = T

At time

t = t^{0} = 0

, the maximum

ϕ (x_{i}, y_{j}, z_{k}, 0)

\max_{(i, j, k) \in J} ψ (x_{i}, y_{j}, z_{k}, 0)

Then, from (39), the maximum at time

t = t^{1} = τ

\max_{\begin{matrix} (i, j, k) \in J \end{matrix}} | ϕ_{i, j, k}^{1} | \leq \max \{| α_{1} |, | α_{2} |, | β_{1} |, | β_{2} |, | η_{1} |, | η_{2} |, \max_{\begin{matrix} (i, j, k) \in J \end{matrix}} | ψ_{i, j, k} | + τ \max_{\begin{matrix} (i, j, k) \in J \\ n \in K \end{matrix}} | q_{i, j, k}^{0} |\} .

(40)

Similarly, the maximum at time

t = t^{2} = 2 τ

\begin{matrix} \max_{\begin{matrix} (i, j, k) \in J \end{matrix}} | ϕ_{i, j, k}^{2} | & \leq & \max \{| α_{1} |, | α_{2} |, | β_{1} |, | β_{2} |, | η_{1} |, | η_{2} |, \max_{\begin{matrix} (i, j, k) \in J \end{matrix}} | ϕ_{i, j, k}^{1} | + τ \max_{\begin{matrix} (i, j, k) \in J \\ n \in K \end{matrix}} | q_{i, j, k}^{1} |\} \\ = & \max {| α_{1} |, | α_{2} |, | β_{1} |, | β_{2} |, | η_{1} |, | η_{2} |, \\ \max_{\begin{matrix} (i, j, k) \in J \end{matrix}} | ψ_{i, j, k} | + τ (\max_{\begin{matrix} (i, j, k) \in J \\ n \in K \end{matrix}} | q_{i, j, k}^{0} | + \max_{\begin{matrix} (i, j, k) \in J \\ n \in K \end{matrix}} | q_{i, j, k}^{1} |)}, \\ \leq & \max {| α_{1} |, | α_{2} |, | β_{1} |, | β_{2} |, | η_{1} |, | η_{2} |, \\ \max_{\begin{matrix} (i, j, k) \in J \end{matrix}} | ψ_{i, j, k} | + 2 τ (\max_{\begin{matrix} (i, j, k) \in J \\ n \in K \end{matrix}} | q_{i, j, k}^{0} | + \max_{\begin{matrix} (i, j, k) \in J \\ n \in K \end{matrix}} | q_{i, j, k}^{1} |)} . \end{matrix}

(41)

Then, by iteration we find that the maximum at time

t = t^{n}

\begin{matrix} \max_{\begin{matrix} (i, j, k) \in J \end{matrix}} | ϕ_{i, j, k}^{n} | & \leq & \max {| α_{1} |, | α_{2} |, | β_{1} |, | β_{2} |, | η_{1} |, | η_{2} |, \\ \max_{\begin{matrix} (i, j, k) \in J \end{matrix}} | ψ_{i, j, k} | + n τ \max (\max_{\begin{matrix} (i, j, k) \in J \\ n \in K \end{matrix}} | q_{i, j, k}^{0} |, \dots, \max_{\begin{matrix} (i, j, k) \in J \\ n \in K \end{matrix}} | q_{i, j, k}^{n - 1} |)}, \end{matrix}

(42)

and the maximum at time

t = t^{N} = T

\begin{matrix} \max_{\begin{matrix} (i, j, k) \in J \end{matrix}} | ϕ_{i, j, k}^{N} | & \leq & \max {| α_{1} |, | α_{2} |, | β_{1} |, | β_{2} |, | η_{1} |, | η_{2} |, \\ \max_{\begin{matrix} (i, j, k) \in J \end{matrix}} | ψ_{i, j, k} | + N τ \max (\max_{\begin{matrix} (i, j, k) \in J \\ n \in K \end{matrix}} | q_{i, j, k}^{0} |, \dots, \max_{\begin{matrix} (i, j, k) \in J \\ n \in K \end{matrix}} | q_{i, j, k}^{N - 1} |)}, \end{matrix}

(43)

Then, as in Eq. (16), for

n \in K

, we can define

\begin{matrix} {∥ f ∥}_{H_{h}} & = & \max_{n \in K} {| α_{1} |, | α_{2} |, | β_{1} |, | β_{2} |, | η_{1} |, | η_{2} |, \\ \max_{\begin{matrix} (i, j, k) \in J \end{matrix}} | ψ_{i, j} | + n τ \max {\max_{\begin{matrix} (i, j, k) \in J \end{matrix}} | q_{i, j, k}^{0} |, \dots, \max_{\begin{matrix} (i, j, k) \in J \end{matrix}} | q_{i, j, k}^{n - 1} |}} . \end{matrix}

(44)

The above inequality is true for all

n \in K

, and hence we have that

∥ ϕ^{h} ∥_{Ω_{h}} \leq C {∥ f^{h} ∥}_{H_{h}},

(45)

where C = 1 and ∥ ϕ^{h} ∥_{Ω_{h}} = \max_{\begin{matrix} i, j, k \in J \\ n \in K \end{matrix}} | ϕ_{i, j, k}^{n} | .

(46)

Hence, inequality (17) is satisfied and the solution is stable. This completes the proof.

4. Numerical Results

In this section, we use the implicit forward time central space (FTCS) finite difference method to solve two cases (homogeneous and nonhomogeneous) of three-dimensional ADRE problems of air pollution. In three dimensions, the advection-diffusion-reaction equation can be written as follows:

\begin{matrix} \frac{\partial C}{\partial t} + u \frac{\partial C}{\partial x} + v \frac{\partial C}{\partial y} + w \frac{\partial C}{\partial z} = D_{x} \frac{\partial^{2} C}{\partial x^{2}} + D_{y} \frac{\partial^{2} C}{\partial y^{2}} + D_{z} \frac{\partial^{2} C}{\partial z^{2}} - R C (x, y, z, t) + Q (x, y, z, t), \end{matrix}

(47)

where

C = C (x, y, z, t)

[kg/m

^{3}

] is the air pollutant concentration at

(x, y, z)

[m] and time t s, u, v, w are wind velocity components [m/s] in

x -

y -

and

z -

directions, respectively,

D_{x}

and

D_{y}

are constant diffusion coefficients in the horizontal direction [m

^{2}

/s],

D_{z}

is a constant diffusion coefficient in the

z -

direction (vertical) [m

^{2}

/s], R is a reaction coefficient, and

Q (x, y, z, t)

is the rate of change of concentrations of sources or sinks of air pollutants [kg/m

^{3} \cdot

s].

In this section, we use an implicit forward time central space (FTCS) scheme to solve the problems. Using the results from section 2, we can write the finite diference equation as

\begin{matrix} \frac{C_{i, j, k}^{n + 1} - C_{i, j, k}^{n}}{Δ t} + u \frac{C_{i + 1, j, k}^{n + 1} - C_{i - 1, j, k}^{n + 1}}{2 Δ x} + v \frac{C_{i, j + 1, k}^{n + 1} - C_{i, j - 1, k}^{n + 1}}{2 Δ y} + w \frac{C_{i, j, k + 1}^{n + 1} - C_{i, j, k - 1}^{n + 1}}{2 Δ z} \\ - D_{x} [\frac{C_{i + 1, j, k}^{n + 1} - 2 C_{i, j, k}^{n + 1} + C_{i - 1, j, k}^{n + 1}}{{(Δ x)}^{2}}] - D_{y} [\frac{C_{i, j + 1, k}^{n + 1} - 2 C_{i, j, k}^{n + 1} + C_{i, j - 1, k}^{n + 1}}{{(Δ y)}^{2}}] \\ - D_{z} [\frac{C_{i, j, k + 1}^{n + 1} - 2 C_{i, j, k}^{n + 1} + C_{i, j, k - 1}^{n + 1}}{{(Δ z)}^{2}}] + R C_{i, j, k}^{n + 1} = Q_{i, j, k}^{n} . \end{matrix}

(48)

We assume that the initial condition is

C (x, y, z, 0) = 0, 0 \leq x \leq 1; 0 \leq y \leq 1; 0 \leq z \leq 1,

(49)

and that the boundary conditions are as shown in Table 1.

After rearrangement, we can rewrite (48) in the form

\begin{matrix} C_{i, j, k}^{n + 1} + S_{1} C_{i - 1, j, k}^{n + 1} + S_{2} C_{i, j - 1, k}^{n + 1} + S_{3} C_{i, j, k - 1}^{n + 1} + S_{4} C_{i + 1, j, k}^{n + 1} + S_{5} C_{i, j + 1, k}^{n + 1} + S_{6} C_{i, j, k + 1}^{n + 1} \\ = S_{7} C_{i, j, k}^{n} + S_{8} Q_{i, j, k}^{n}, \end{matrix}

(50)

where

\begin{matrix} S_{1} = (- \frac{u Δ t}{2 Δ x} - \frac{D_{x} Δ t}{{(Δ x)}^{2}}) / (1 + \frac{2 D_{x} Δ t}{{(Δ x)}^{2}} + \frac{2 D_{y} Δ t}{{(Δ y)}^{2}} + \frac{2 D_{z} Δ t}{{(Δ z)}^{2}} + R Δ t) \\ S_{2} = (- \frac{v Δ t}{2 Δ y} - \frac{D_{y} Δ t}{{(Δ y)}^{2}}) / (1 + \frac{2 D_{x} Δ t}{{(Δ x)}^{2}} + \frac{2 D_{y} Δ t}{{(Δ y)}^{2}} + \frac{2 D_{z} Δ t}{{(Δ z)}^{2}} + R Δ t) \\ S_{3} = (- \frac{w Δ t}{2 Δ z} - \frac{D_{z} Δ t}{{(Δ z)}^{2}}) / (1 + \frac{2 D_{x} Δ t}{{(Δ x)}^{2}} + \frac{2 D_{y} Δ t}{{(Δ y)}^{2}} + \frac{2 D_{z} Δ t}{{(Δ z)}^{2}} + R Δ t) \\ S_{4} = (\frac{u Δ t}{2 Δ x} - \frac{D_{x} Δ t}{{(Δ x)}^{2}}) / (1 + \frac{2 D_{x} Δ t}{{(Δ x)}^{2}} + \frac{2 D_{y} Δ t}{{(Δ y)}^{2}} + \frac{2 D_{z} Δ t}{{(Δ z)}^{2}} + R Δ t) \end{matrix}

\begin{matrix} S_{5} = (\frac{v Δ t}{2 Δ y} - \frac{D_{y} Δ t}{{(Δ y)}^{2}}) / (1 + \frac{2 D_{x} Δ t}{{(Δ x)}^{2}} + \frac{2 D_{y} Δ t}{{(Δ y)}^{2}} + \frac{2 D_{z} Δ t}{{(Δ z)}^{2}} + R Δ t) \\ S_{6} = (\frac{w Δ t}{2 Δ z} - \frac{D_{z} Δ t}{{(Δ z)}^{2}}) / (1 + \frac{2 D_{x} Δ t}{{(Δ x)}^{2}} + \frac{2 D_{y} Δ t}{{(Δ y)}^{2}} + \frac{2 D_{z} Δ t}{{(Δ z)}^{2}} + R Δ t) \\ S_{7} = 1 / (1 + \frac{2 D_{x} Δ t}{{(Δ x)}^{2}} + \frac{2 D_{y} Δ t}{{(Δ y)}^{2}} + \frac{2 D_{z} Δ t}{{(Δ z)}^{2}} + R Δ t) \\ S_{8} = Δ t / (1 + \frac{2 D_{x} Δ t}{{(Δ x)}^{2}} + \frac{2 D_{y} Δ t}{{(Δ y)}^{2}} + \frac{2 D_{z} Δ t}{{(Δ z)}^{2}} + R Δ t) . \end{matrix}

Case 1 : We consider the three-dimensional advection-diffusion equation (47). By taking

w = 0, R = 0, Q = 0, L = 1, W = 1

and

H = 1,

this reduces to a two-dimensional advection three-dimensional diffusion equation for transport of pollutants in the street tunnel problem discussed in [12,25]

\begin{matrix} \frac{\partial C}{\partial t} + u \frac{\partial C}{\partial x} + v \frac{\partial C}{\partial y} = D_{x} \frac{\partial^{2} C}{\partial x^{2}} + D_{y} \frac{\partial^{2} C}{\partial y^{2}} + D_{z} \frac{\partial^{2} C}{\partial z^{2}}, 0 < t < T . \end{matrix}

(51)

This system of partial differential equation and initial condition (49) and boundary condition Table 1 comes from a model of pollution distribution in a street tunnel, where the wind flows steadily in the x and y directions and there is no flux of pollutant through the solid side-walls or the solid base and roof of the tunnel. We use the nondimensionalised parameters

Δ x = Δ y = Δ z = 0.1, Δ t = 0.005, D_{x} = D_{y} = D_{z} = 0.02, 0.06, 0.2, 0.5, u = 0.6, v = 0.4

and time

T = 20

. The numerical solutions and the contour plots of the problem obtained using the implicit FTCS scheme in Eq. (48) are plotted in Figure 1, Figure 2, Figure 3 and Figure 4. We have tested the accuracy of our numerical solutions by comparing them with that of Kusuma et al [12] and found that there was good agreement.

Case 2: In this problem, We consider the three-dimensional advection-diffusion equation (47). By taking

w = 0, L = 1, W = 1

and

H = 1,

with a decay rate (R) and source term (Q). The advection-diffusion-reaction equation that we consider is as follows.

\begin{matrix} \frac{\partial C}{\partial t} + u \frac{\partial C}{\partial x} + v \frac{\partial C}{\partial y} = D_{x} \frac{\partial^{2} C}{\partial x^{2}} + D_{y} \frac{\partial^{2} C}{\partial y^{2}} + D_{z} \frac{\partial^{2} C}{\partial z^{2}} - R C + Q (x, y, t) . \end{matrix}

(52)

This case of partial differential equation and initial condition (49) and boundary condition Table 1 comes from a model of pollution distribution in a street tunnel, where the wind flows steadily in the x and y directions and there is no flux of pollutant through the solid side-walls or the solid base and roof of the tunnel. We use the nondimensionalised parameters

Δ x = Δ y = Δ z = 0.1, Δ t = 0.005, D_{x} = D_{y} = D_{z} = 0.2, u = 0.6, v = 0.4

and time

T = 20

We have computed the numerical solutions of the above system using the implicit FTCS method for the following cases.

$R = 0.05$ and $Q = 0, 0.007, - e^{- t}, 0.001 sin (x t)$ . The results for this case are shown in Figure 5 (a) - (d) at a height $z = 0.2$ meters for $Q = 0, 0.007, - e^{- t}$ and $Q = 0.001 sin (x t)$ , respectively.
$R = 0.1$ and $Q = 0, 0.007, - e^{- t}, 0.001 sin (x t)$ . The results for this case are shown in Figure 6 (a) - (d) at a height $z = 0.2$ meters for $Q = 0, 0.007, - e^{- t}$ and $Q = 0.001 sin (x t)$ , respectively.
$R = 0.5$ and $Q = 0, 0.007, - e^{- t}, 0.001 sin (x t)$ . The results for this case are shown in Figure 7 (a) - (d) at a height $z = 0.2$ meters for $Q = 0, 0.007, - e^{- t}$ and $Q = 0.001 sin (x t)$ , respectively.
$z = 0.2$ m for $R = 0.05$ and $Q = 0, 0.007, - e^{- t}, 0.001 sin (x t)$ . The results for this case are shown in Figure 8 (a) at $y = 0.5$ m and (b) at $x = 0.5$ m.
$z = 0.2$ m for $Q = 0.07$ and $R = 0.05, 0.1, 0.5$ . The results for this case are shown in Figure 9 (a) at $y = 0.5$ m and (b) at $x = 0.5$ m.

It can be seen that the numerical solutions for the different Q and R values show similar qualitative behavior but with differences in detail.

5. Conclusions

In this work, we have studied a finite difference method (FDM) for obtaining numerical solutions of advection-diffusion-reaction equations (ADRE) and applied the methods to solve air pollutant problems. For the finite difference method we have used the implicit forward time central space (FTCS) scheme. For the applications, we have solved 3-D air pollutant concentration problems for tunnels. We have investigated the stability, consistency and convergence of the finite difference method. We applied the implicit FTCS finite difference method to solve 3-D ADRE models for three problems of air pollution in tunnels for a range of different source terms Q and reaction terms R. For case 1, we checked the accuracy of our solutions by comparing them with previously published results. We assumed that

Q = R = 0

, that there was advection in the x and y directions and diffusion in the x, y and z directions. We then compared our numerical solutions of

D_{x} = D_{y} = D_{z} = 0.2

with Kusuma et al [12] and found that there was good agreement. For case 2, we assumed that Q and R were non-zero, that there was advection in the

x -

and

y -

direction and diffusion in the x, y and z directions. We found that the effects of the Q and R were to rapidly reduce the pollutant concentrations in the tunnel.

Author Contributions

Conceptualization, S.S., K.P., E.J.M., S.S. and S.P.; methodology, S.S., K.P., E.J.M., S.S. and S.P.; software, S.S., K.P. and E.J.M.; validation, S.S., K.P. and E.J.M.; formal analysis, S.S., K.P., E.J.M. and S.S.; writing-review and editing, S.S., K.P. and E.J.M. All authors have read and agreed to the published version of manuscript.

Funding

The research was funded by Faculty of Applied Science, King Mongkut’s University of Technology North Bangkok, Thailand contract no. 661132.

Acknowledgments

The authors would like to express their sincere thanks to the Department of Mathematics, Faculty of Applied Science at King Mongkut’s University of Technology North Bangkok for supporting us to do this research.

Conflicts of Interest

The authors declare no conflict of interest.

References

Di Pietro, D.A.; Ern, A. Mathematical Aspects of Discontinuous Galerkin Methods 2012.
Savović, S.; Djordjevich, A. Finite difference solution of the one-dimensional advection–diffusion equation with variable coefficients in semi-infinite media. International Journal of Heat and Mass Transfer 2012, 55, 4291–4294. [Google Scholar] [CrossRef]
Appadu, A.R.; Gidey, H. Time-splitting procedures for the numerical solution of the 2D advection-diffusion equation. Mathematical Problems in Engineering 2013, 2013. [Google Scholar] [CrossRef]
Bause, M.; Schwegler, K. Higher order finite element approximation of systems of convection–diffusion–reaction equations with small diffusion. Journal of computational and applied mathematics 2013, 246, 52–64. [Google Scholar] [CrossRef]
Mojtabi, A.; Deville, M.O. One-dimensional linear advection–diffusion equation: Analytical and finite element solutions. Computers & Fluids 2015, 107, 189–195. [Google Scholar]
Gharehbaghi, A. Explicit and implicit forms of differential quadrature method for advection–diffusion equation with variable coefficients in semi-infinite domain. Journal of Hydrology 2016, 541, 935–940. [Google Scholar] [CrossRef]
Gyrya, V.; Lipnikov, K. The arbitrary order mimetic finite difference method for a diffusion equation with a non-symmetric diffusion tensor. Journal of Computational Physics 2017, 348, 549–566. [Google Scholar] [CrossRef]
Bahar, E.; Gürarslan, G. Numerical solution of advection-diffusion equation using operator splitting method. International Journal of Engineering and Applied Sciences 2017, 9, 76–88. [Google Scholar] [CrossRef]
Oyjinda, P.; Pochai, N. Numerical Simulation to Air Pollution Emission Control near an Industrial Zone. Advances in Mathematical Physics 2017, 2017. [Google Scholar] [CrossRef]
Suebyat, K.; Pochai, N. Numerical Simulation for a Three-Dimensional Air Pollution Measurement Model in a Heavy Traffic Area under the Bangkok Sky Train Platform. Abstract and Applied Analysis. Hindawi, 2018, Vol. 2018.
Al-Jawary, M.A.; Azeez, M.M.; Radhi, G.H. Analytical and numerical solutions for the nonlinear Burgers and advection–diffusion equations by using a semi-analytical iterative method. Computers & Mathematics with Applications 2018, 76, 155–171. [Google Scholar]
Kusuma, J.; Ribal, A.; Mahie, A.G. On FTCS Approach for Box Model of Three-Dimension Advection-Diffusion Equation. International Journal of Differential Equations 2018, 2018. [Google Scholar] [CrossRef]
Lou, J.; Li, L.; Luo, H.; Nishikawa, H. Reconstructed discontinuous Galerkin methods for linear advection–diffusion equations based on first-order hyperbolic system. Journal of Computational Physics 2018, 369, 103–124. [Google Scholar] [CrossRef]
Bhatt, H.; Khaliq, A.; Wade, B. Efficient Krylov-based exponential time differencing method in application to 3D advection-diffusion-reaction systems. Applied Mathematics and Computation 2018, 338, 260–273. [Google Scholar] [CrossRef]
Pananu, K.; Sungnul, S.; Sirisubtawee, S.; Phongthanapanich, S. Convergence and Applications of the Implicit Finite Difference Method for Advection-Diffusion-Reaction Equations. IAENG International Journal of Computer Science 2020, 47, 645–663. [Google Scholar]
Cheng, H.; Zheng, G. Analyzing 3D advection-diffusion problems by using the improved element-free Galerkin method. Mathematical Problems in Engineering 2020, 2020. [Google Scholar] [CrossRef]
Cruz-Quintero, E.; Jurado, F. Boundary Control for a Certain Class of Reaction-Advection-Diffusion System. Mathematics 2020, 8, 1854. [Google Scholar] [CrossRef]
Para, K.; Jitsom, B.; Eymard, R.; Sungnul, S.; Sirisubtawee, S.; Phongthanapanich, S. An accuracy comparison of piecewise linear reconstruction techniques for the characteristic finite volume method for two-dimensional convection-diffusion equation. ZAMM-Journal of Applied Mathematics and Mechanics/Zeitschrift für Angewandte Mathematik und Mechanik 2021, 101, e201900245. [Google Scholar] [CrossRef]
Sun, H.; Xu, Y.; Lin, J.; Zhang, Y. A space-time backward substitution method for three-dimensional advection-diffusion equations. Computers & Mathematics with Applications 2021, 97, 77–85. [Google Scholar]
Hidayat, M.I.P. Meshless finite difference method with B-splines for numerical solution of coupled advection-diffusion-reaction problems. International Journal of Thermal Sciences 2021, 165, 106933. [Google Scholar] [CrossRef]
Shahid, N.; Rehman, M.A.u.; Khalid, A.; Fatima, U.; Shaikh, T.S.; Ahmed, N.; Alotaibi, H.; Rafiq, M.; Khan, I.; Nisar, K.S. Mathematical analysis and numerical investigation of advection-reaction-diffusion computer virus model. Results in Physics 2021, 26, 104294. [Google Scholar] [CrossRef]
Murillo-García, O.F.; Jurado, F. Adaptive Boundary Control for a Certain Class of Reaction–Advection–Diffusion System. Mathematics 2021, 9, 2224. [Google Scholar] [CrossRef]
Para, K.; Sungnul, S.; Sirisubtawee, S.; Phongthanapanich, S. A Comparison of Numerical Solutions for Advection-Diffusion-Reaction Equations between Finite Volume and Finite Difference Methods. Engineering Letters 2022, 30. [Google Scholar]
Tannehill, J.C.; Anderson, D.; Pletcher, R.H. Computational Fluid Mechanics and Heat Transfer, 2nd ed.; Taylor & Francis, 1997. [Google Scholar]
Thongmoon, M.; McKibbin, R.; Tangmanee, S. Numerical solution of a 3-D advection-dispersion model for pollutant transport. Thai Journal of Mathematics 2012, 5, 91–108. [Google Scholar]

Figure 1. Pollutant distribution (a) and contour plot (b) for case 1 in a street tunnel with advection only in x and y directions at

D_{x} = D_{y} = D_{z} = 0.02

Figure 1. Pollutant distribution (a) and contour plot (b) for case 1 in a street tunnel with advection only in x and y directions at

D_{x} = D_{y} = D_{z} = 0.02

Figure 2. Pollutant distribution (a) and contour plot (b) for case 1 in a street tunnel with advection only in x and y directions at

D_{x} = D_{y} = D_{z} = 0.06

Figure 2. Pollutant distribution (a) and contour plot (b) for case 1 in a street tunnel with advection only in x and y directions at

D_{x} = D_{y} = D_{z} = 0.06

Figure 3. Pollutant distribution (a) and contour plot (b) for case 1 in a street tunnel with advection only in x and y directions at

D_{x} = D_{y} = D_{z} = 0.2

Figure 3. Pollutant distribution (a) and contour plot (b) for case 1 in a street tunnel with advection only in x and y directions at

D_{x} = D_{y} = D_{z} = 0.2

Figure 4. Pollutant distribution (a) and contour plot (b) for case 1 in a street tunnel with advection only in x and y directions at

D_{x} = D_{y} = D_{z} = 0.5

Figure 4. Pollutant distribution (a) and contour plot (b) for case 1 in a street tunnel with advection only in x and y directions at

D_{x} = D_{y} = D_{z} = 0.5

Figure 5. Numerical solution of air pollutant concentration and contour plot for case 2 for

R = 0.05

and (a)

Q = 0

, (b)

Q = 0.007

, (c)

Q = - e^{- t}

and (d)

Q = 0.001 sin (x t)

Figure 5. Numerical solution of air pollutant concentration and contour plot for case 2 for

R = 0.05

and (a)

Q = 0

, (b)

Q = 0.007

, (c)

Q = - e^{- t}

and (d)

Q = 0.001 sin (x t)

Figure 6. Numerical solution of air pollutant concentration and contour plot for case 2 for

R = 0.1

and (a)

Q = 0

, (b)

Q = 0.007

, (c)

Q = - e^{- t}

and (d)

Q = 0.001 sin (x t)

Figure 6. Numerical solution of air pollutant concentration and contour plot for case 2 for

R = 0.1

and (a)

Q = 0

, (b)

Q = 0.007

, (c)

Q = - e^{- t}

and (d)

Q = 0.001 sin (x t)

Figure 7. Numerical solution of air pollutant concentration and contour plot for case 2 for

R = 0.5

and (a)

Q = 0

, (b)

Q = 0.007

, (c)

Q = - e^{- t}

and (d)

Q = 0.001 sin (x t)

Figure 7. Numerical solution of air pollutant concentration and contour plot for case 2 for

R = 0.5

and (a)

Q = 0

, (b)

Q = 0.007

, (c)

Q = - e^{- t}

and (d)

Q = 0.001 sin (x t)

Figure 8. Comparison of numerical solutions for Case 2 at

z = 0.2

m for

R = 0.05

and

Q = 0, 0.007, - e^{- t}, 0.001 sin (x t)

at (a)

y = 0.5

m and (b)

x = 0.5

Figure 8. Comparison of numerical solutions for Case 2 at

z = 0.2

m for

R = 0.05

and

Q = 0, 0.007, - e^{- t}, 0.001 sin (x t)

at (a)

y = 0.5

m and (b)

x = 0.5

Figure 9. Comparison of numerical solutions for Case 2 at

z = 0.2

m for

Q = 0.07

R = 0.05, 0.1, 0.5

at (a)

y = 0.5

m and (b)

x = 0.5

Figure 9. Comparison of numerical solutions for Case 2 at

z = 0.2

m for

Q = 0.07

R = 0.05, 0.1, 0.5

at (a)

y = 0.5

m and (b)

x = 0.5

Table 1. Boundary conditions.

Boundary	Boundary condition	Value
Entrance gate : $x = 0, 0 \leq y < 0.5, 0 \leq z \leq 1$	$C (0, y, z, t)$	0
Entrance gate : $x = 0, 0.5 \leq y \leq 1, 0 \leq z \leq 1$	$C (0, y, z, t)$	1
Exit gate : $x = 1, 0 \leq y \leq 1, 0 \leq z \leq 1$	$\frac{\partial C}{\partial x} (1, y, z, t)$	0
Right side wall : $0 \leq x < 0.3, y = 0, 0 \leq z \leq 1$	$C (x, 0, z, t)$	0
Right side wall : $0.3 \leq x \leq 0.6, y = 0, 0 \leq z \leq 1$	$C (x, 0, z, t)$	1
Right side wall : $0.6 < x \leq 1, y = 0, 0 \leq z \leq 1$	$C (x, 0, z, t)$	0
Left side wall : $0 < x < 1, y = 1, 0 \leq z \leq 1$	$\frac{\partial C}{\partial y} (x, 1, z, t)$	0
Ground : $0 < x < 1, 0 < y < 1, z = 0$	$\frac{\partial C}{\partial z} (x, y, 0, t)$	0
Ceiling : $0 < x < 1, 0 < y < 1, z = 1$	$\frac{\partial C}{\partial z} (x, y, 1, t)$	0

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

© 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

Copyright: This open access article is published under a Creative Commons CC BY 4.0 license, which permit the free download, distribution, and reuse, provided that the author and preprint are cited in any reuse.

Submitted:

29 August 2023

Posted:

18 September 2023

You are already at the latest version

Alerts

Abstract

Keywords:

Subject: Computer Science and Mathematics - Computational Mathematics

1. Introduction

\begin{matrix} \frac{\partial ϕ}{\partial t} + \nabla \cdot (v ϕ - ε \nabla ϕ) + κ ϕ = q, \end{matrix}

(1)

where

ϕ

is a scalar quantity,

v = v (x)

is a given advection velocity vector,

ε \geq 0

is a diffusion coefficient,

κ

is a reaction coefficient,

q = q (x, t)

is a prescribed source term,

x

is a position vector and time

t \in [0, T]

, where

T > 0

2. Convergence Theory

We consider the following initial time, space boundary value problem (IBVP) for the advection-diffusion-reaction equation

L ϕ (x, t) = \frac{\partial ϕ (x, t)}{\partial t} + \nabla \cdot (v (x) ϕ (x, t) - ε \nabla ϕ (x, t)) + κ ϕ (x, t) = q (x, t),

(2)

where

(x, t) \in Ω \times (0, T]

, and with the initial time condition

ϕ (x, 0) = ψ (x) for x \in Ω,

(3)

and either Dirichlet or Neumann boundary space conditions

\begin{matrix} Dirichlet & ϕ (x, t) = ξ_{D} (x, t) \\ Neumann & ε \frac{\partial ϕ (x, t)}{\partial n} = ξ_{N} (x, t) for (x, t) \in \partial Ω \times [0, T] . \end{matrix}

(4)

In Eqs. (2)–(4), T is a finite time,

\partial Ω

is the boundary of the spatial domain

Ω

ξ_{D}

and

ξ_{N}

are given functions,

v (x)

is an advection vector,

ε > 0

is a diffusion coefficient,

κ

is a reaction constant and

q (x, t)

is a source term. In addition, L is a linear partial differential operator from a space of continuous functions F to a space of continuous functions H, where the dependent variable

ϕ (x, t) \in F

, and where

L ϕ (x, t) \in H

and

q (x, t) \in H

The first step in developing the finite difference method is to discretize

Ω

\partial Ω

and

[0, T]

. We assume that there are

N_{1}

step sizes of length

Δ x = \frac{L_{1}}{N_{1}}

in the x direction,

N_{2}

step sizes of length

Δ y = \frac{L_{2}}{N_{2}}

in the y direction and

N_{3}

step sizes of length

Δ z = \frac{L_{3}}{N_{3}}

in the z direction. For simplicity, we will assume that

Δ x = Δ y = Δ z = h

and therefore

L_{1} = N_{1} h

L_{2} = N_{2} h

and

L_{3} = N_{3} h

. Then the discretized version of

Ω

is the spatial grid

Ω_{h} = {(x_{i}, y_{j}, z_{k}) | x_{i} = i h, y_{j} = j h, z_{k} = k h, i \in J_{1}, j \in J_{2}, k \in J_{3}},

(5)

where

J_{1} = {0, 1, 2, . . ., N_{1}}

J_{2} = {0, 1, 2, . . ., N_{2}}

J_{3} = {0, 1, 2, . . ., N_{3}}

. The boundary of the spatial grid is

\begin{matrix} \partial Ω_{h} & = & {x_{0} = 0; y_{j} = j h, j \in J_{2}, z_{k} = k h, k \in J_{3}} \\ \cup {x_{N_{1}} = L_{1}; y_{j} = j h, j \in J_{2}, z_{k} = k h, k \in J_{3}} \\ \cup {x_{i} = i h, i \in J_{1}; y_{0} = 0; z_{k} = k h, k \in J_{3}} \\ \cup {x_{i} = i h, i \in J_{1}; y_{N_{2}} = L_{2}; z_{k} = k h, k \in J_{3}} \\ \cup {x_{i} = i h, i \in J_{1}, y_{j} = j h, j \in J_{2}; z_{0} = 0} \\ \cup {x_{i} = i h, i \in J_{1}, y_{j} = j h, j \in J_{2}; z_{N_{3}} = L_{3}} . \end{matrix}

(6)

Similarly, the interval

[0, T]

is discretized with N step sizes of size

Δ t = \frac{T}{N}

. The set of discretized times is then

T_{K} = {t^{n} | t^{n} = n Δ t, n \in K}, K = {0, 1, 2, \dots, N} .

(7)

Finally, we define

F_{h}

and

H_{h}

as discretized function spaces with domain

Ω_{h} \times T_{K}

. We also define a projection

P_{h} (ϕ) \in F_{h}

of the exact solution

ϕ (x, t)

of (2) onto

Ω_{h} \times T_{K}

and a projection

P_{h} (q) \in H_{h}

of the source term

q (x, t)

onto

Ω_{h} \times T_{K}

We can then write a finite difference scheme (FDS) corresponding to the continuous partial differential equation (2) in the form

L_{h} ϕ^{h} (x, t) = q^{h} (x, t), (x, t) \in Ω_{h} \times T_{K},

(8)

where

L_{h} : F_{h} \to H_{h}

is a finite difference operator, corresponding to the partial differential operator L, which acts from the discrete function space

F_{h}

to the discrete function space

H_{h}

. The initial condition in (3) can then be rewritten as the discretized initial condition

ϕ^{h} (x, 0) = ψ^{h} (x), x \in Ω_{h},

(9)

and, for the Dirichlet case, the boundary conditions in (4) can be rewritten as the discretized boundary conditions

ϕ^{h} (x, t) = ξ^{h} (x, t) for (x, t) \in \partial Ω_{h} \times T_{K} .

(10)

Lax’s equivalence theorem [24]

Using the definition given above of the projection

P_{h}

of a function in the continuous space H onto the discrete space

H_{h}

, we define the truncation error (T.E.) between the finite difference system Eq. (8) and the IBVP Eq. (2) as follows

T . E . = P_{h} (L ϕ (x, t)) - L_{h} ϕ^{h} (x, t), (x, t) \in Ω_{h} \times T_{K} .

(11)

Introducing the norm in the set of discrete functions as the infinity norm

\begin{matrix} \begin{matrix} ∥ ϕ^{h} ∥_{F_{h}} & = \max_{\begin{matrix} i, j . k \in J \\ n \in K \end{matrix}} | ϕ^{h} (x_{i}, y_{j}, z_{k}, t^{n}) |, \end{matrix} \end{matrix}

(12)

we have the definition of convergence for numerical methods as follows

Definition 1.

Convergence

A solution

ϕ^{h}

of the FDS in Eq. (8)-(10) converges to the solution ϕ of the IBVP in Eq. (2)-(4) if

∥ P_{h} (ϕ) - ϕ^{h} ∥_{F_{h}} \to 0 a s h \to 0 .

(13)

Definition 2.

Convergence with order m

The FDS in Eq. (8) converges with order m if

∥ P_{h} (ϕ) - ϕ^{h} ∥_{F_{h}} \leq C h^{m},

(14)

where C is a positive constant that does not depend on h.

1) Consistency : A finite difference representation of a PDE is said to be consistent [24] if the difference between the PDE and its difference representation vanishes as the mesh is refined, i.e.,

Definition 3.

The FDS (8) is a consistent approximation of the IBVP (2)-(4) with order m if

{∥ T . E . ∥}_{H_{h}} \leq C h^{m},

(15)

where C is a positive constant that does not depend on h.

2) Stability : The finite difference scheme defined by (8) with linear operator

L_{h}

will be called stable, if there exists

h_{0} > 0

such that for arbitrary

h < h_{0}

and for

∥ f^{h} ∥_{H_{h}} = \max_{\begin{matrix} i, j, k \in J \end{matrix}} | ψ (x_{i}, y_{j}, z_{k}, t^{n}) | + \max_{\begin{matrix} i, j, k \in J \\ n \in K - {0} \end{matrix}} | q (x_{i}, y_{j}, z_{k}, t^{n}) |,

(16)

the solution

ϕ^{h}

of the FDS,

L_{h} ϕ^{h} = q^{h}

in Eq. (8), exists and is unique and satisfies the inequality

∥ ϕ^{h} ∥_{F_{h}} \leq C {∥ f^{h} ∥}_{H_{h}},

(17)

where the positive constant C does not depend on h.

The inequality (17) has to be true for any initial conditions

ψ (x_{i}, y_{j}, z_{,} 0)

and source terms

q (x_{i}, y_{j}, z_{k}, t^{n})

of Eq. (2). In particular, if

q (x_{i}, y_{j}, z_{k}, t^{n}) = 0

, then the condition (17) reduces to the necessary condition for stability of the homogeneous version of (8).

For the special case of Fourier or Von Neumann Analysis [24], a solution of the FDS (8) for

q (x_{i}, y_{j}, z_{k}, t_{n}) = 0

can be written in the form

ϕ^{h} (x_{i}, y_{j}, z_{k}, t^{n}) = λ^{n} e^{I (α i + β j + η k)}, (i, j, k) \in J, n \in K, I = \sqrt{- 1},

(18)

where

α

β

and

η

are wave numbers and

e^{I (α i + β j + η k)}

are eigenfunctions corresponding to eigenvalues

λ

L_{h}

. The necessary condition for stability of FDS (8) for

q (x_{i}, y_{j}, z_{k}, t^{n}) = 0

will then hold for all

α, β, η \in ℜ

if the following inequality holds:

| λ | \leq 1 .

(19)

3. Consistency and Stability Analysis

Consider the following initial-boundary value problem

\begin{matrix} \frac{\partial ϕ}{\partial t} + u \frac{\partial ϕ}{\partial x} + v \frac{\partial ϕ}{\partial y} + w \frac{\partial ϕ}{\partial z} - ε (\frac{\partial^{2} ϕ}{\partial x^{2}} + \frac{\partial^{2} ϕ}{\partial y^{2}} + \frac{\partial^{2} ϕ}{\partial z^{2}}) + κ ϕ = q (x, y, z, t), \end{matrix}

\begin{matrix} (x, y, z) \in (0, L_{1}) \times (0, L_{2}) \times (0, L_{3}), t \in (0, T], \\ subject to & ϕ (x, y, z, 0) = ψ (x, y, z), \\ ϕ (0, y, z, t) = α_{1}, ϕ (L_{1}, y, z, t) = α_{2}, \\ ϕ (x, 0, z, t) = β_{1}, ϕ (x, L_{2}, z, t) = β_{2}, \end{matrix}

(20)

\begin{matrix} ϕ (x, y, 0, t) = η_{1}, ϕ (x, y, L_{3}, t) = η_{2}, \end{matrix}

(21)

in which the above partial differential equation is the component form of (1). Here

ϕ = ϕ (x, y, z, t)

is the dependent variable and

v = (u, v, w)

is the advection velocity where

u, v

and w are the velocities in the x-, y- and z-directions, respectively.

In this section, we will analyze the convergence of the finite difference method with the implicit FTCS scheme applied to (20). We discretize the domain of the problem:

(x_{i}, y_{j}, z_{k}, t^{n}) \in Ω_{h} \times T_{K}

, where the spatial grid

Ω_{h}

is defined in (5) and the time grid

T_{K} = [0, T]

has points

t^{n} = n Δ t

with

n \in K = {0, 1, 2, . . ., N} .

\begin{matrix} L_{h} ϕ^{h} = & \frac{ϕ_{i, j, k}^{n + 1} - ϕ_{i, j, k}^{n}}{Δ t} + u \frac{ϕ_{i + 1, j, k}^{n + 1} - ϕ_{i - 1, j, k}^{n + 1}}{2 Δ x} + v \frac{ϕ_{i, j + 1, k}^{n + 1} - ϕ_{i, j - 1, k}^{n + 1}}{2 Δ y} + w \frac{ϕ_{i, j, k + 1}^{n + 1} - ϕ_{i, j, k - 1}^{n + 1}}{2 Δ z} \\ - ε [\frac{ϕ_{i + 1, j, k}^{n + 1} - 2 ϕ_{i, j, k}^{n + 1} + ϕ_{i - 1, j, k}^{n + 1}}{{(Δ x)}^{2}} + \frac{ϕ_{i, j + 1, k}^{n + 1} - 2 ϕ_{i, j, k}^{n + 1} + ϕ_{i, j - 1, k}^{n + 1}}{{(Δ y)}^{2}} \end{matrix}

(22)

\begin{matrix} + \frac{ϕ_{i, j, k + 1}^{n + 1} - 2 ϕ_{i, j, k}^{n + 1} + ϕ_{i, j, k - 1}^{n + 1}}{{(Δ z)}^{2}}] + κ ϕ_{i, j, k}^{n + 1} = q_{i, j, k}^{n}, \\ subject to & ϕ (x_{i}, y_{j}, z_{k}, 0) = ψ (x_{i}, y_{j}, z_{k}), \\ ϕ (0, y_{j}, z_{k}, t^{n}) = α_{1}, ϕ (L_{1}, y_{j}, z_{k}, t^{n}) = α_{2}, \\ ϕ (x_{i}, 0, z_{k}, t^{n}) = β_{1}, ϕ (x_{i}, L_{2}, z_{k}, t^{n}) = β_{2}, \end{matrix}

(23)

\begin{matrix} ϕ (x_{i}, y_{j}, 0, t^{n}) = η_{1}, ϕ (x_{i}, y_{j}, L_{3}, t^{n}) = η_{2}, \end{matrix}

where

ϕ_{i, j, k}^{n} = ϕ^{h} (x_{i}, y_{j}, z_{k}, t^{n})

and

q_{i, j, k}^{n} = q (x_{i}, y_{j}, z_{k}, t^{n})

The truncation error for (22) obtained using the method is

\begin{matrix} T . E . & = & - \frac{Δ t}{2!} \frac{\partial^{2} ϕ}{\partial t^{2}} - u \frac{{(Δ x)}^{2}}{3!} \frac{\partial^{3} ϕ}{\partial x^{3}} - v \frac{{(Δ y)}^{2}}{3!} \frac{\partial^{3} ϕ}{\partial y^{3}} - w \frac{{(Δ z)}^{2}}{3!} \frac{\partial^{3} ϕ}{\partial z^{3}} \\ + ε (\frac{2 {(Δ x)}^{2}}{4!} \frac{\partial^{4} ϕ}{\partial x^{4}} + \frac{2 {(Δ y)}^{2}}{4!} \frac{\partial^{4} ϕ}{\partial y^{4}} + \frac{2 {(Δ z)}^{2}}{4!} \frac{\partial^{4} ϕ}{\partial z^{4}}) \\ + O ({(Δ t)}^{2}, {(Δ x)}^{4}, {(Δ y)}^{4}, {(Δ z)}^{4}) . \end{matrix}

(24)

Therefore, from (24), we have

\begin{matrix} T . E . = O ((Δ t), {(Δ x)}^{2}, {(Δ y)}^{2}, {(Δ z)}^{2}) . \end{matrix}

Taking

(Δ t, Δ x, Δ y, Δ z) \to 0

, then we have

\begin{matrix} lim_{(Δ t, Δ x, Δ y, Δ z) \to 0} ∥ T . E . ∥ = lim_{(Δ t, Δ x, Δ y, Δ z) \to 0} ∥ O ((Δ t), {(Δ x)}^{2}, {(Δ y)}^{2}, {(Δ z)}^{2}) ∥ = 0 . \end{matrix}

τ = Δ t, h = Δ x = Δ y = Δ z

Case I: Homogeneous equation

Setting

q_{i, j, k}^{n} = 0

in equation (22), we obtain the discretized homogeneous version of (22). Assuming a discretized solution of the resulting equation as

ϕ_{i, j, k}^{n} = λ^{n} e^{I (α i + β j + η k)}

I = \sqrt{- 1}

[24], we obtain the following equation for the eigenvalue

λ

\begin{matrix} \frac{λ^{n + 1} e^{I (α i + β j + η k)} - λ^{n} e^{I (α i + β j + η k)}}{τ} + u \frac{λ^{n + 1} e^{I (α (i + 1) + β j + η k)} - λ^{n + 1} e^{I (α (i - 1) + β j + η k)}}{2 h} \\ + v \frac{λ^{n + 1} e^{I (α i + β (j + 1) + η k)} - λ^{n + 1} e^{I (α i + β (j - 1) + η k)}}{2 h} \\ + w \frac{λ^{n + 1} e^{I (α i + β j + η (k + 1))} - λ^{n + 1} e^{I (α i + β j + η (k - 1))}}{2 h} \end{matrix}

\begin{matrix} - ε [\frac{λ^{n + 1} e^{I (α (i + 1) + β j + η k)} - 2 λ^{n + 1} e^{I (α i + β j + η k)} + λ^{n + 1} e^{I (α (i - 1) + β j + η k)}}{h^{2}} \\ + \frac{λ^{n + 1} e^{I (α i + β (j + 1) + η k)} - 2 λ^{n + 1} e^{I (α i + β j + η k)} + λ^{n + 1} e^{I (α i + β (j - 1) + η k)}}{h^{2}} \\ + \frac{λ^{n + 1} e^{I (α i + β j + η (k + 1))} - 2 λ^{n + 1} e^{I (α i + β j + η k)} + λ^{n + 1} e^{I (α i + β j + η (k - 1))}}{h^{2}}] \\ + κ λ^{n + 1} e^{I (α i + β j + η k)} = 0 . \end{matrix}

(25)

Then, we can solve (25) for the eigenvalue

λ

as follows.

\begin{matrix} λ^{n} e^{I (α i + β j + η k)} [\frac{λ - 1}{τ}] - λ^{n + 1} e^{I (α i + β j + η k)} [- u \frac{e^{I α} - e^{- I α}}{2 h} - v \frac{e^{I β} - e^{- I β}}{2 h} \\ - w \frac{e^{I η} - e^{- I η}}{2 h} + ε \frac{e^{I α} - 2 + e^{- I α}}{h^{2}} + ε \frac{e^{I β} - 2 + e^{- I β}}{h^{2}} + ε \frac{e^{I η} - 2 + e^{- I η}}{h^{2}} - κ] = 0, \\ λ = \frac{1}{1 + \frac{τ}{h} [\frac{4 ε}{h} ({sin}^{2} \frac{α}{2} + {sin}^{2} \frac{β}{2} + {sin}^{2} \frac{η}{2}) + I (u sin α + v sin β + w sin η) + κ h]} . \end{matrix}

Letting

γ = \frac{4 ε τ}{h^{2}}

and

μ = \frac{τ}{h}

, we then obtain,

\begin{matrix} λ = \frac{1}{1 + γ ({sin}^{2} \frac{α}{2} + {sin}^{2} \frac{β}{2} + {sin}^{2} \frac{η}{2}) + μ I (u sin α + v sin β + w sin η) + κ τ} . \end{matrix}

Consequently, the magnitude of

λ

\begin{matrix} | λ | & = & \frac{1}{\sqrt{{[1 + γ ({sin}^{2} \frac{α}{2} + {sin}^{2} \frac{β}{2} + {sin}^{2} \frac{η}{2}) + κ τ]}^{2} + {[μ (u sin α + v sin β + w sin η)]}^{2}}}, \\ \leq & 1 . \end{matrix}

Therefore, condition (19) for Von Neumann stability, which is independent of

τ

and h, has been proved.

Case II: Nonhomogeneous equation

We now prove that condition (17) is satisfied for the nonhomogeneous equation (22) with

q_{i, j, k}^{n} \neq 0

by finding a solution

ϕ_{i, j, k}^{n + 1}

for the linear system of difference equations (22) and (23) if

ϕ_{i, j, k}^{n}

is assumed known. Substituting

Δ x = Δ y = Δ z = h

and

Δ t = τ

in (22) and (23), we obtain

\begin{matrix} L_{h} ϕ^{h} & = & \frac{ϕ_{i, j, k}^{n + 1} - ϕ_{i, j, k}^{n}}{τ} + u \frac{ϕ_{i + 1, j, k}^{n + 1} - ϕ_{i - 1, j, k}^{n + 1}}{2 h} + v \frac{ϕ_{i, j + 1, k}^{n + 1} - ϕ_{i, j - 1, k}^{n + 1}}{2 h} + w \frac{ϕ_{i, j, k + 1}^{n + 1} - ϕ_{i, j, k - 1}^{n + 1}}{2 h} \\ - ε [\frac{ϕ_{i + 1, j, k}^{n + 1} - 2 ϕ_{i, j, k}^{n + 1} + ϕ_{i - 1, j, k}^{n + 1}}{h^{2}} + \frac{ϕ_{i, j + 1, k}^{n + 1} - 2 ϕ_{i, j, k}^{n + 1} + ϕ_{i, j - 1, k}^{n + 1}}{h^{2}} \end{matrix}

\begin{matrix} + \frac{ϕ_{i, j, k + 1}^{n + 1} - 2 ϕ_{i, j, k}^{n + 1} + ϕ_{i, j, k - 1}^{n + 1}}{h^{2}} + \frac{ϕ_{i, j, k + 1}^{n + 1} - 2 ϕ_{i, j, k}^{n + 1} + ϕ_{i, j, k - 1}^{n + 1}}{h^{2}}] + κ ϕ_{i, j, k}^{n + 1} = q_{i, j, k}^{n} \\ subject to ϕ_{0, j, k}^{n} = α_{1}, ϕ_{N_{1}, j, k}^{n} = α_{2}; j = 0, 1, . . ., N_{2}; k = 0, 1, . . ., N_{3} \\ ϕ_{i, 0, k}^{n} = β_{1}, ϕ_{i, N_{2}, k}^{n} = β_{2}; i = 0, 1, . . ., N_{1}; k = 0, 1, . . ., N_{3} \end{matrix}

(26)

\begin{matrix} ϕ_{i, j, 0}^{n} = η_{1}, ϕ_{i, j, N_{3}}^{n} = η_{2}; i = 0, 1, . . ., N_{1}; j = 0, 1, . . ., N_{2} \end{matrix}

(27)

After rearranging, we can rewrite (26) in the form

a ϕ_{i - 1, j, k}^{n + 1} + b ϕ_{i, j - 1, k}^{n + 1} + c ϕ_{i, j, k - 1}^{n + 1} + d ϕ_{i, j, k}^{n + 1} + e ϕ_{i + 1, j, k}^{n + 1} + f ϕ_{i, j + 1, k}^{n + 1} + g ϕ_{i, j, k + 1}^{n + 1} = ϕ_{i, j, k}^{n},

(28)

\begin{matrix} where & a = (\frac{u τ}{2 h} + \frac{ε τ}{h^{2}}), b = (\frac{v τ}{2 h} + \frac{ε τ}{h^{2}}), c = (\frac{w τ}{2 h} + \frac{ε τ}{h^{2}}), \\ d = (- 1 - \frac{6 ε τ}{h^{2}} - κ τ), e = (- \frac{u τ}{2 h} + \frac{ε τ}{h^{2}}), f = (- \frac{v τ}{2 h} + \frac{ε τ}{h^{2}}), \\ g = (- \frac{w τ}{2 h} + \frac{ε τ}{h^{2}}), ϕ_{i, j, k}^{n} = - ϕ_{i, j, k}^{n} - τ q_{i, j, k}^{n}, \end{matrix}

and where the step sizes

τ

and h can be chosen so that the coefficients

a, b, c, d, e, f, g

satisfy the conditions

| d | > | a | + | b | + | c | + | e | + | f | + | g | + δ, δ > 0 .

(29)

Lemma 1.

If the coefficients of Eq.(28) satisfy the conditions Eq.(29) then the solution of Eq.(28) exists and is unique and satisfies the inequality

| ϕ_{i, j, k}^{n + 1} | \leq \max \{| α_{1} |, | α_{2} |, | β_{1} |, | β_{2} |, | η_{1} |, | η_{2} |, \frac{1}{δ} \max_{p, q, r} | ϕ_{p, q, r}^{n} |\} .

(30)

Proof

First, we prove inequality (30). Assume that

| ϕ_{p, q, r}^{n + 1} | = \max {| ϕ_{i, j, k}^{n + 1} |, i = 0, 1, . . ., N_{1}, j = 0, 1, . . ., N_{2}, k = 0, 1, . . ., N_{3}}

Let

0 < p < N_{1}, 0 < q < N_{2}, 0 < r < N_{3}

. Then

\begin{matrix} | d | | ϕ_{p, q, r}^{n + 1} | & = & | - a ϕ_{p - 1, q, r}^{n + 1} - b ϕ_{p, q - 1, r}^{n + 1} - c ϕ_{p, q, r - 1}^{n + 1} - e ϕ_{p + 1, q, r}^{n + 1} - f ϕ_{p, q + 1, r}^{n + 1} - g ϕ_{p, q, r + 1}^{n + 1} + ϕ_{p, q, r}^{n} | \\ \leq & | a ϕ_{p - 1, q, r}^{n + 1} | + | b ϕ_{p, q - 1, r}^{n + 1} | + | c ϕ_{p, q, r - 1}^{n + 1} | + | e ϕ_{p + 1, q, r}^{n + 1} | + | f ϕ_{p, q + 1, r}^{n + 1} | + | g ϕ_{p, q, r + 1}^{n + 1} | \\ + | ϕ_{p, q, r}^{n} | \\ = & | a | | ϕ_{p - 1, q, r}^{n + 1} | + | b | | ϕ_{p, q - 1, r}^{n + 1} | + | c | | ϕ_{p, q, r - 1}^{n + 1} | + | e | | ϕ_{p + 1, q, r}^{n + 1} | + | f | | ϕ_{p, q + 1, r}^{n + 1} | \\ + | g | | ϕ_{p, q, r + 1}^{n + 1} | + | ϕ_{p, q, r}^{n} | \\ \leq & (| a | + | b | + | c | + | e | + | f | + | g |) | ϕ_{p, q, r}^{n + 1} | + | ϕ_{p, q, r}^{n} |, \\ | ϕ_{p, q, r}^{n + 1} | & \leq & \frac{| ϕ_{p, q, r}^{n} |}{| d | - | a | - | b | - | c | - | e | - | f | - | g |} \\ \leq & \frac{| ϕ_{p, q, r}^{n} |}{δ} . \end{matrix}

This completes the proof.

Now we will prove the stability of the implicit FTCS scheme by showing that inequality (17) is satisfied.

Proof

Let

i^{*}

j^{*}

and

k^{*}

be the smallest non-negative integers such that

| ϕ_{i^{*}, j^{*}, k^{*}}^{n + 1} | = \max_{\begin{matrix} (i, j, k) \in J \end{matrix}} | ϕ_{i, j, k}^{n + 1} |,

(31)

where

J = {(i, j, k) | i = 0, 1, 2, . . ., N_{1}, j = 0, 1, 2, . . ., N_{2}, k = 0, 1, 2, . . ., N_{3}} .

Now let

\bar{J}

be the discretized boundary set of J, i.e.,

\begin{matrix} \bar{J} & = & {(i, j, k) \in J | i = 0, N_{1}} \cup {(i, j, k) \in J | i \neq 0, i \neq N_{1}, j = 0, N_{2}} \\ \cup {(i, j, k) \in J | i \neq 0, i \neq N_{1}, j \neq 0, j \neq N_{2}, k = 0, N_{3}} . \end{matrix}

(32)

It is obvious that if

(i^{*}, j^{*}, k^{*}) \in \bar{J},

then we have

\max_{(i, j, k) \in J} | ϕ_{i, j, k}^{n + 1} | \leq \max {| α_{1} |, | α_{2} |, | β_{1} |, | β_{2} |, | η_{1} |, | η_{2} |} .

(33)

(i^{*}, j^{*}, k^{*}) \in J = J - \bar{J},

then replacing

(i, j, k)

in (26) with

(i^{*}, j^{*}, k^{*})

we obtain

\begin{matrix} - ϕ_{i^{*}, j^{*}, k^{*}}^{n} - τ q_{i^{*}, j^{*}, k^{*}}^{n} = (\frac{u τ}{2 h} + \frac{ε τ}{h^{2}}) ϕ_{i^{*} - 1, j^{*}, k^{*}}^{n + 1} + (\frac{v τ}{2 h} + \frac{ε τ}{h^{2}}) ϕ_{i^{*}, j^{*} - 1, k^{*}}^{n + 1} + (\frac{w τ}{2 h} + \frac{ε τ}{h^{2}}) ϕ_{i^{*}, j^{*}, k^{*} - 1}^{n + 1} \\ - (1 + \frac{6 ε τ}{h^{2}} + κ τ) ϕ_{i^{*}, j^{*}, k^{*}}^{n + 1} - (\frac{u τ}{2 h} - \frac{ε τ}{h^{2}}) ϕ_{i^{*} + 1, j^{*}, k^{*}}^{n + 1} \\ - (\frac{v τ}{2 h} - \frac{ε τ}{h^{2}}) ϕ_{i^{*}, j^{*} + 1, k^{*}}^{n + 1} - (\frac{w τ}{2 h} - \frac{ε τ}{h^{2}}) ϕ_{i^{*}, j^{*}, k^{*} + 1}^{n + 1} . \end{matrix}

(34)

Next we make the following hypothesis:

\begin{matrix} (H 1) : \begin{matrix} u [ϕ_{i^{*} - 1, j^{*}, k^{*}}^{n + 1} - ϕ_{i^{*} + 1, j^{*}, k^{*}}^{n + 1}] \leq 0, \\ v [ϕ_{i^{*}, j^{*} - 1, k^{*}}^{n + 1} - ϕ_{i^{*}, j^{*} + 1, k^{*}}^{n + 1}] \leq 0, \\ w [ϕ_{i^{*}, j^{*}, k^{*} - 1}^{n + 1} - ϕ_{i^{*}, j^{*}, k^{*} + 1}^{n + 1}] \leq 0 . \end{matrix} \end{matrix}

(35)

Without loss of generality we assume that

ϕ_{i^{*}, j^{*}, k^{*}}^{n + 1} > 0

and further assume that the hypothesis (H1) holds. We can estimate the right side of (34) as

\begin{matrix} \frac{ε τ}{h^{2}} [ϕ_{i^{*} - 1, j^{*}, k^{*}}^{n + 1} - ϕ_{i^{*}, j^{*}, k^{*}}^{n + 1}] + \frac{ε τ}{h^{2}} [ϕ_{i^{*}, j^{*} - 1, k^{*}}^{n + 1} - ϕ_{i^{*}, j^{*}, k^{*}}^{n + 1}] + \frac{ε τ}{h^{2}} [ϕ_{i^{*}, j^{*}, k^{*} - 1}^{n + 1} - ϕ_{i^{*}, j^{*}, k^{*}}^{n + 1}] \\ + \frac{ε τ}{h^{2}} [ϕ_{i^{*} + 1, j^{*}, k^{*}}^{n + 1} - ϕ_{i^{*}, j^{*}, k^{*}}^{n + 1}] + \frac{ε τ}{h^{2}} [ϕ_{i^{*}, j^{*} + 1}^{n + 1} - ϕ_{i^{*}, j^{*}, k^{*}}^{n + 1}] + \frac{ε τ}{h^{2}} [ϕ_{i^{*}, j^{*}, k^{*} + 1}^{n + 1} - ϕ_{i^{*}, j^{*}, k^{*}}^{n + 1}] \\ + \frac{u τ}{2 h} [ϕ_{i^{*} - 1, j^{*}, k^{*}}^{n + 1} - ϕ_{i^{*} + 1, j^{*}, k^{*}}^{n + 1}] + \frac{v τ}{2 h} [ϕ_{i^{*}, j^{*} - 1, k^{*}}^{n + 1} - ϕ_{i^{*}, j^{*} + 1, k^{*}}^{n + 1}] \\ + \frac{w τ}{2 h} [ϕ_{i^{*}, j^{*}, k^{*} - 1}^{n + 1} - ϕ_{i^{*}, j^{*}, k^{*} + 1}^{n + 1}] - ϕ_{i^{*}, j^{*}, k^{*}}^{n + 1} - κ τ ϕ_{i^{*}, j^{*}, k^{*}}^{n + 1} \\ \leq - ϕ_{i^{*}, j^{*}, k^{*}}^{n + 1} . \end{matrix}

(36)

Then, substituting (36) into the right hand side of (34), we obtain

ϕ_{i^{*}, j^{*}, k^{*}}^{n + 1} \leq ϕ_{i^{*}, j^{*}, k^{*}}^{n} + τ q_{i^{*}, j^{*}, k^{*}}^{n} .

(37)

Therefore,

\begin{matrix} \max_{\begin{matrix} (i, j, k) \in J \end{matrix}} | ϕ_{i, j, k}^{n + 1} | \leq \max_{\begin{matrix} (i, j, k) \in J \end{matrix}} | ϕ_{i, j, k}^{n} + τ q_{i, j, k}^{n} | \leq \max_{\begin{matrix} (i, j, k) \in J \end{matrix}} | ϕ_{i, j, k}^{n} | + τ \max_{\begin{matrix} (i, j, k) \in J \\ n \in K \end{matrix}} | q_{i, j, k}^{n} |, \end{matrix}

(38)

where

K = {0, 1, 2, . . ., N} .

We now establish the following inequality representing the maximum principle as

\begin{matrix} \max_{\begin{matrix} (i, j, k) \in J \end{matrix}} | ϕ_{i, j, k}^{n + 1} | \leq \max \{| α_{1} |, | α_{2} |, | β_{1} |, | β_{2} |, | η_{1} |, | η_{2} |, \max_{\begin{matrix} (i, j, k) \in J \end{matrix}} | ϕ_{i, j, k}^{n} | + τ \max_{\begin{matrix} (i, j, k) \in J \\ n \in K \end{matrix}} | q_{i, j, k}^{n} |\} . \end{matrix}

(39)

We will now prove that for the given initial and boundary conditions in (27), the maximum solution

\max_{(i, j, k) \in J} | ϕ_{i, j, k}^{n} |

is bounded for all

n \leq N

for finite

t^{N} = T

At time

t = t^{0} = 0

, the maximum

ϕ (x_{i}, y_{j}, z_{k}, 0)

\max_{(i, j, k) \in J} ψ (x_{i}, y_{j}, z_{k}, 0)

Then, from (39), the maximum at time

t = t^{1} = τ

\max_{\begin{matrix} (i, j, k) \in J \end{matrix}} | ϕ_{i, j, k}^{1} | \leq \max \{| α_{1} |, | α_{2} |, | β_{1} |, | β_{2} |, | η_{1} |, | η_{2} |, \max_{\begin{matrix} (i, j, k) \in J \end{matrix}} | ψ_{i, j, k} | + τ \max_{\begin{matrix} (i, j, k) \in J \\ n \in K \end{matrix}} | q_{i, j, k}^{0} |\} .

(40)

Similarly, the maximum at time

t = t^{2} = 2 τ

\begin{matrix} \max_{\begin{matrix} (i, j, k) \in J \end{matrix}} | ϕ_{i, j, k}^{2} | & \leq & \max \{| α_{1} |, | α_{2} |, | β_{1} |, | β_{2} |, | η_{1} |, | η_{2} |, \max_{\begin{matrix} (i, j, k) \in J \end{matrix}} | ϕ_{i, j, k}^{1} | + τ \max_{\begin{matrix} (i, j, k) \in J \\ n \in K \end{matrix}} | q_{i, j, k}^{1} |\} \\ = & \max {| α_{1} |, | α_{2} |, | β_{1} |, | β_{2} |, | η_{1} |, | η_{2} |, \\ \max_{\begin{matrix} (i, j, k) \in J \end{matrix}} | ψ_{i, j, k} | + τ (\max_{\begin{matrix} (i, j, k) \in J \\ n \in K \end{matrix}} | q_{i, j, k}^{0} | + \max_{\begin{matrix} (i, j, k) \in J \\ n \in K \end{matrix}} | q_{i, j, k}^{1} |)}, \\ \leq & \max {| α_{1} |, | α_{2} |, | β_{1} |, | β_{2} |, | η_{1} |, | η_{2} |, \\ \max_{\begin{matrix} (i, j, k) \in J \end{matrix}} | ψ_{i, j, k} | + 2 τ (\max_{\begin{matrix} (i, j, k) \in J \\ n \in K \end{matrix}} | q_{i, j, k}^{0} | + \max_{\begin{matrix} (i, j, k) \in J \\ n \in K \end{matrix}} | q_{i, j, k}^{1} |)} . \end{matrix}

(41)

Then, by iteration we find that the maximum at time

t = t^{n}

\begin{matrix} \max_{\begin{matrix} (i, j, k) \in J \end{matrix}} | ϕ_{i, j, k}^{n} | & \leq & \max {| α_{1} |, | α_{2} |, | β_{1} |, | β_{2} |, | η_{1} |, | η_{2} |, \\ \max_{\begin{matrix} (i, j, k) \in J \end{matrix}} | ψ_{i, j, k} | + n τ \max (\max_{\begin{matrix} (i, j, k) \in J \\ n \in K \end{matrix}} | q_{i, j, k}^{0} |, \dots, \max_{\begin{matrix} (i, j, k) \in J \\ n \in K \end{matrix}} | q_{i, j, k}^{n - 1} |)}, \end{matrix}

(42)

and the maximum at time

t = t^{N} = T

\begin{matrix} \max_{\begin{matrix} (i, j, k) \in J \end{matrix}} | ϕ_{i, j, k}^{N} | & \leq & \max {| α_{1} |, | α_{2} |, | β_{1} |, | β_{2} |, | η_{1} |, | η_{2} |, \\ \max_{\begin{matrix} (i, j, k) \in J \end{matrix}} | ψ_{i, j, k} | + N τ \max (\max_{\begin{matrix} (i, j, k) \in J \\ n \in K \end{matrix}} | q_{i, j, k}^{0} |, \dots, \max_{\begin{matrix} (i, j, k) \in J \\ n \in K \end{matrix}} | q_{i, j, k}^{N - 1} |)}, \end{matrix}

(43)

Then, as in Eq. (16), for

n \in K

, we can define

\begin{matrix} {∥ f ∥}_{H_{h}} & = & \max_{n \in K} {| α_{1} |, | α_{2} |, | β_{1} |, | β_{2} |, | η_{1} |, | η_{2} |, \\ \max_{\begin{matrix} (i, j, k) \in J \end{matrix}} | ψ_{i, j} | + n τ \max {\max_{\begin{matrix} (i, j, k) \in J \end{matrix}} | q_{i, j, k}^{0} |, \dots, \max_{\begin{matrix} (i, j, k) \in J \end{matrix}} | q_{i, j, k}^{n - 1} |}} . \end{matrix}

(44)

The above inequality is true for all

n \in K

, and hence we have that

∥ ϕ^{h} ∥_{Ω_{h}} \leq C {∥ f^{h} ∥}_{H_{h}},

(45)

where C = 1 and ∥ ϕ^{h} ∥_{Ω_{h}} = \max_{\begin{matrix} i, j, k \in J \\ n \in K \end{matrix}} | ϕ_{i, j, k}^{n} | .

(46)

Hence, inequality (17) is satisfied and the solution is stable. This completes the proof.

4. Numerical Results

\begin{matrix} \frac{\partial C}{\partial t} + u \frac{\partial C}{\partial x} + v \frac{\partial C}{\partial y} + w \frac{\partial C}{\partial z} = D_{x} \frac{\partial^{2} C}{\partial x^{2}} + D_{y} \frac{\partial^{2} C}{\partial y^{2}} + D_{z} \frac{\partial^{2} C}{\partial z^{2}} - R C (x, y, z, t) + Q (x, y, z, t), \end{matrix}

(47)

where

C = C (x, y, z, t)

[kg/m

^{3}

] is the air pollutant concentration at

(x, y, z)

[m] and time t s, u, v, w are wind velocity components [m/s] in

x -

y -

and

z -

directions, respectively,

D_{x}

and

D_{y}

are constant diffusion coefficients in the horizontal direction [m

^{2}

/s],

D_{z}

is a constant diffusion coefficient in the

z -

direction (vertical) [m

^{2}

/s], R is a reaction coefficient, and

Q (x, y, z, t)

is the rate of change of concentrations of sources or sinks of air pollutants [kg/m

^{3} \cdot

s].

In this section, we use an implicit forward time central space (FTCS) scheme to solve the problems. Using the results from section 2, we can write the finite diference equation as

\begin{matrix} \frac{C_{i, j, k}^{n + 1} - C_{i, j, k}^{n}}{Δ t} + u \frac{C_{i + 1, j, k}^{n + 1} - C_{i - 1, j, k}^{n + 1}}{2 Δ x} + v \frac{C_{i, j + 1, k}^{n + 1} - C_{i, j - 1, k}^{n + 1}}{2 Δ y} + w \frac{C_{i, j, k + 1}^{n + 1} - C_{i, j, k - 1}^{n + 1}}{2 Δ z} \\ - D_{x} [\frac{C_{i + 1, j, k}^{n + 1} - 2 C_{i, j, k}^{n + 1} + C_{i - 1, j, k}^{n + 1}}{{(Δ x)}^{2}}] - D_{y} [\frac{C_{i, j + 1, k}^{n + 1} - 2 C_{i, j, k}^{n + 1} + C_{i, j - 1, k}^{n + 1}}{{(Δ y)}^{2}}] \\ - D_{z} [\frac{C_{i, j, k + 1}^{n + 1} - 2 C_{i, j, k}^{n + 1} + C_{i, j, k - 1}^{n + 1}}{{(Δ z)}^{2}}] + R C_{i, j, k}^{n + 1} = Q_{i, j, k}^{n} . \end{matrix}

(48)

We assume that the initial condition is

C (x, y, z, 0) = 0, 0 \leq x \leq 1; 0 \leq y \leq 1; 0 \leq z \leq 1,

(49)

and that the boundary conditions are as shown in Table 1.

After rearrangement, we can rewrite (48) in the form

\begin{matrix} C_{i, j, k}^{n + 1} + S_{1} C_{i - 1, j, k}^{n + 1} + S_{2} C_{i, j - 1, k}^{n + 1} + S_{3} C_{i, j, k - 1}^{n + 1} + S_{4} C_{i + 1, j, k}^{n + 1} + S_{5} C_{i, j + 1, k}^{n + 1} + S_{6} C_{i, j, k + 1}^{n + 1} \\ = S_{7} C_{i, j, k}^{n} + S_{8} Q_{i, j, k}^{n}, \end{matrix}

(50)

where

\begin{matrix} S_{1} = (- \frac{u Δ t}{2 Δ x} - \frac{D_{x} Δ t}{{(Δ x)}^{2}}) / (1 + \frac{2 D_{x} Δ t}{{(Δ x)}^{2}} + \frac{2 D_{y} Δ t}{{(Δ y)}^{2}} + \frac{2 D_{z} Δ t}{{(Δ z)}^{2}} + R Δ t) \\ S_{2} = (- \frac{v Δ t}{2 Δ y} - \frac{D_{y} Δ t}{{(Δ y)}^{2}}) / (1 + \frac{2 D_{x} Δ t}{{(Δ x)}^{2}} + \frac{2 D_{y} Δ t}{{(Δ y)}^{2}} + \frac{2 D_{z} Δ t}{{(Δ z)}^{2}} + R Δ t) \\ S_{3} = (- \frac{w Δ t}{2 Δ z} - \frac{D_{z} Δ t}{{(Δ z)}^{2}}) / (1 + \frac{2 D_{x} Δ t}{{(Δ x)}^{2}} + \frac{2 D_{y} Δ t}{{(Δ y)}^{2}} + \frac{2 D_{z} Δ t}{{(Δ z)}^{2}} + R Δ t) \\ S_{4} = (\frac{u Δ t}{2 Δ x} - \frac{D_{x} Δ t}{{(Δ x)}^{2}}) / (1 + \frac{2 D_{x} Δ t}{{(Δ x)}^{2}} + \frac{2 D_{y} Δ t}{{(Δ y)}^{2}} + \frac{2 D_{z} Δ t}{{(Δ z)}^{2}} + R Δ t) \end{matrix}

\begin{matrix} S_{5} = (\frac{v Δ t}{2 Δ y} - \frac{D_{y} Δ t}{{(Δ y)}^{2}}) / (1 + \frac{2 D_{x} Δ t}{{(Δ x)}^{2}} + \frac{2 D_{y} Δ t}{{(Δ y)}^{2}} + \frac{2 D_{z} Δ t}{{(Δ z)}^{2}} + R Δ t) \\ S_{6} = (\frac{w Δ t}{2 Δ z} - \frac{D_{z} Δ t}{{(Δ z)}^{2}}) / (1 + \frac{2 D_{x} Δ t}{{(Δ x)}^{2}} + \frac{2 D_{y} Δ t}{{(Δ y)}^{2}} + \frac{2 D_{z} Δ t}{{(Δ z)}^{2}} + R Δ t) \\ S_{7} = 1 / (1 + \frac{2 D_{x} Δ t}{{(Δ x)}^{2}} + \frac{2 D_{y} Δ t}{{(Δ y)}^{2}} + \frac{2 D_{z} Δ t}{{(Δ z)}^{2}} + R Δ t) \\ S_{8} = Δ t / (1 + \frac{2 D_{x} Δ t}{{(Δ x)}^{2}} + \frac{2 D_{y} Δ t}{{(Δ y)}^{2}} + \frac{2 D_{z} Δ t}{{(Δ z)}^{2}} + R Δ t) . \end{matrix}

Case 1 : We consider the three-dimensional advection-diffusion equation (47). By taking

w = 0, R = 0, Q = 0, L = 1, W = 1

and

H = 1,

this reduces to a two-dimensional advection three-dimensional diffusion equation for transport of pollutants in the street tunnel problem discussed in [12,25]

\begin{matrix} \frac{\partial C}{\partial t} + u \frac{\partial C}{\partial x} + v \frac{\partial C}{\partial y} = D_{x} \frac{\partial^{2} C}{\partial x^{2}} + D_{y} \frac{\partial^{2} C}{\partial y^{2}} + D_{z} \frac{\partial^{2} C}{\partial z^{2}}, 0 < t < T . \end{matrix}

(51)

Δ x = Δ y = Δ z = 0.1, Δ t = 0.005, D_{x} = D_{y} = D_{z} = 0.02, 0.06, 0.2, 0.5, u = 0.6, v = 0.4

and time

T = 20

Case 2: In this problem, We consider the three-dimensional advection-diffusion equation (47). By taking

w = 0, L = 1, W = 1

and

H = 1,

with a decay rate (R) and source term (Q). The advection-diffusion-reaction equation that we consider is as follows.

\begin{matrix} \frac{\partial C}{\partial t} + u \frac{\partial C}{\partial x} + v \frac{\partial C}{\partial y} = D_{x} \frac{\partial^{2} C}{\partial x^{2}} + D_{y} \frac{\partial^{2} C}{\partial y^{2}} + D_{z} \frac{\partial^{2} C}{\partial z^{2}} - R C + Q (x, y, t) . \end{matrix}

(52)

Δ x = Δ y = Δ z = 0.1, Δ t = 0.005, D_{x} = D_{y} = D_{z} = 0.2, u = 0.6, v = 0.4

and time

T = 20

We have computed the numerical solutions of the above system using the implicit FTCS method for the following cases.

$R = 0.05$ and $Q = 0, 0.007, - e^{- t}, 0.001 sin (x t)$ . The results for this case are shown in Figure 5 (a) - (d) at a height $z = 0.2$ meters for $Q = 0, 0.007, - e^{- t}$ and $Q = 0.001 sin (x t)$ , respectively.
$R = 0.1$ and $Q = 0, 0.007, - e^{- t}, 0.001 sin (x t)$ . The results for this case are shown in Figure 6 (a) - (d) at a height $z = 0.2$ meters for $Q = 0, 0.007, - e^{- t}$ and $Q = 0.001 sin (x t)$ , respectively.
$R = 0.5$ and $Q = 0, 0.007, - e^{- t}, 0.001 sin (x t)$ . The results for this case are shown in Figure 7 (a) - (d) at a height $z = 0.2$ meters for $Q = 0, 0.007, - e^{- t}$ and $Q = 0.001 sin (x t)$ , respectively.
$z = 0.2$ m for $R = 0.05$ and $Q = 0, 0.007, - e^{- t}, 0.001 sin (x t)$ . The results for this case are shown in Figure 8 (a) at $y = 0.5$ m and (b) at $x = 0.5$ m.
$z = 0.2$ m for $Q = 0.07$ and $R = 0.05, 0.1, 0.5$ . The results for this case are shown in Figure 9 (a) at $y = 0.5$ m and (b) at $x = 0.5$ m.

It can be seen that the numerical solutions for the different Q and R values show similar qualitative behavior but with differences in detail.

5. Conclusions

Q = R = 0

, that there was advection in the x and y directions and diffusion in the x, y and z directions. We then compared our numerical solutions of

D_{x} = D_{y} = D_{z} = 0.2

with Kusuma et al [12] and found that there was good agreement. For case 2, we assumed that Q and R were non-zero, that there was advection in the

x -

and

y -

direction and diffusion in the x, y and z directions. We found that the effects of the Q and R were to rapidly reduce the pollutant concentrations in the tunnel.

Author Contributions

Funding

The research was funded by Faculty of Applied Science, King Mongkut’s University of Technology North Bangkok, Thailand contract no. 661132.

Acknowledgments

Conflicts of Interest

The authors declare no conflict of interest.

References

Di Pietro, D.A.; Ern, A. Mathematical Aspects of Discontinuous Galerkin Methods 2012.
Savović, S.; Djordjevich, A. Finite difference solution of the one-dimensional advection–diffusion equation with variable coefficients in semi-infinite media. International Journal of Heat and Mass Transfer 2012, 55, 4291–4294. [Google Scholar] [CrossRef]
Appadu, A.R.; Gidey, H. Time-splitting procedures for the numerical solution of the 2D advection-diffusion equation. Mathematical Problems in Engineering 2013, 2013. [Google Scholar] [CrossRef]
Bause, M.; Schwegler, K. Higher order finite element approximation of systems of convection–diffusion–reaction equations with small diffusion. Journal of computational and applied mathematics 2013, 246, 52–64. [Google Scholar] [CrossRef]
Mojtabi, A.; Deville, M.O. One-dimensional linear advection–diffusion equation: Analytical and finite element solutions. Computers & Fluids 2015, 107, 189–195. [Google Scholar]
Gharehbaghi, A. Explicit and implicit forms of differential quadrature method for advection–diffusion equation with variable coefficients in semi-infinite domain. Journal of Hydrology 2016, 541, 935–940. [Google Scholar] [CrossRef]
Gyrya, V.; Lipnikov, K. The arbitrary order mimetic finite difference method for a diffusion equation with a non-symmetric diffusion tensor. Journal of Computational Physics 2017, 348, 549–566. [Google Scholar] [CrossRef]
Bahar, E.; Gürarslan, G. Numerical solution of advection-diffusion equation using operator splitting method. International Journal of Engineering and Applied Sciences 2017, 9, 76–88. [Google Scholar] [CrossRef]
Oyjinda, P.; Pochai, N. Numerical Simulation to Air Pollution Emission Control near an Industrial Zone. Advances in Mathematical Physics 2017, 2017. [Google Scholar] [CrossRef]
Suebyat, K.; Pochai, N. Numerical Simulation for a Three-Dimensional Air Pollution Measurement Model in a Heavy Traffic Area under the Bangkok Sky Train Platform. Abstract and Applied Analysis. Hindawi, 2018, Vol. 2018.
Al-Jawary, M.A.; Azeez, M.M.; Radhi, G.H. Analytical and numerical solutions for the nonlinear Burgers and advection–diffusion equations by using a semi-analytical iterative method. Computers & Mathematics with Applications 2018, 76, 155–171. [Google Scholar]
Kusuma, J.; Ribal, A.; Mahie, A.G. On FTCS Approach for Box Model of Three-Dimension Advection-Diffusion Equation. International Journal of Differential Equations 2018, 2018. [Google Scholar] [CrossRef]
Lou, J.; Li, L.; Luo, H.; Nishikawa, H. Reconstructed discontinuous Galerkin methods for linear advection–diffusion equations based on first-order hyperbolic system. Journal of Computational Physics 2018, 369, 103–124. [Google Scholar] [CrossRef]
Bhatt, H.; Khaliq, A.; Wade, B. Efficient Krylov-based exponential time differencing method in application to 3D advection-diffusion-reaction systems. Applied Mathematics and Computation 2018, 338, 260–273. [Google Scholar] [CrossRef]
Pananu, K.; Sungnul, S.; Sirisubtawee, S.; Phongthanapanich, S. Convergence and Applications of the Implicit Finite Difference Method for Advection-Diffusion-Reaction Equations. IAENG International Journal of Computer Science 2020, 47, 645–663. [Google Scholar]
Cheng, H.; Zheng, G. Analyzing 3D advection-diffusion problems by using the improved element-free Galerkin method. Mathematical Problems in Engineering 2020, 2020. [Google Scholar] [CrossRef]
Cruz-Quintero, E.; Jurado, F. Boundary Control for a Certain Class of Reaction-Advection-Diffusion System. Mathematics 2020, 8, 1854. [Google Scholar] [CrossRef]
Para, K.; Jitsom, B.; Eymard, R.; Sungnul, S.; Sirisubtawee, S.; Phongthanapanich, S. An accuracy comparison of piecewise linear reconstruction techniques for the characteristic finite volume method for two-dimensional convection-diffusion equation. ZAMM-Journal of Applied Mathematics and Mechanics/Zeitschrift für Angewandte Mathematik und Mechanik 2021, 101, e201900245. [Google Scholar] [CrossRef]
Sun, H.; Xu, Y.; Lin, J.; Zhang, Y. A space-time backward substitution method for three-dimensional advection-diffusion equations. Computers & Mathematics with Applications 2021, 97, 77–85. [Google Scholar]
Hidayat, M.I.P. Meshless finite difference method with B-splines for numerical solution of coupled advection-diffusion-reaction problems. International Journal of Thermal Sciences 2021, 165, 106933. [Google Scholar] [CrossRef]
Shahid, N.; Rehman, M.A.u.; Khalid, A.; Fatima, U.; Shaikh, T.S.; Ahmed, N.; Alotaibi, H.; Rafiq, M.; Khan, I.; Nisar, K.S. Mathematical analysis and numerical investigation of advection-reaction-diffusion computer virus model. Results in Physics 2021, 26, 104294. [Google Scholar] [CrossRef]
Murillo-García, O.F.; Jurado, F. Adaptive Boundary Control for a Certain Class of Reaction–Advection–Diffusion System. Mathematics 2021, 9, 2224. [Google Scholar] [CrossRef]
Para, K.; Sungnul, S.; Sirisubtawee, S.; Phongthanapanich, S. A Comparison of Numerical Solutions for Advection-Diffusion-Reaction Equations between Finite Volume and Finite Difference Methods. Engineering Letters 2022, 30. [Google Scholar]
Tannehill, J.C.; Anderson, D.; Pletcher, R.H. Computational Fluid Mechanics and Heat Transfer, 2nd ed.; Taylor & Francis, 1997. [Google Scholar]
Thongmoon, M.; McKibbin, R.; Tangmanee, S. Numerical solution of a 3-D advection-dispersion model for pollutant transport. Thai Journal of Mathematics 2012, 5, 91–108. [Google Scholar]

Figure 1. Pollutant distribution (a) and contour plot (b) for case 1 in a street tunnel with advection only in x and y directions at

D_{x} = D_{y} = D_{z} = 0.02

Figure 1. Pollutant distribution (a) and contour plot (b) for case 1 in a street tunnel with advection only in x and y directions at

D_{x} = D_{y} = D_{z} = 0.02

Figure 2. Pollutant distribution (a) and contour plot (b) for case 1 in a street tunnel with advection only in x and y directions at

D_{x} = D_{y} = D_{z} = 0.06

Figure 2. Pollutant distribution (a) and contour plot (b) for case 1 in a street tunnel with advection only in x and y directions at

D_{x} = D_{y} = D_{z} = 0.06

Figure 3. Pollutant distribution (a) and contour plot (b) for case 1 in a street tunnel with advection only in x and y directions at

D_{x} = D_{y} = D_{z} = 0.2

Figure 3. Pollutant distribution (a) and contour plot (b) for case 1 in a street tunnel with advection only in x and y directions at

D_{x} = D_{y} = D_{z} = 0.2

Figure 4. Pollutant distribution (a) and contour plot (b) for case 1 in a street tunnel with advection only in x and y directions at

D_{x} = D_{y} = D_{z} = 0.5

Figure 4. Pollutant distribution (a) and contour plot (b) for case 1 in a street tunnel with advection only in x and y directions at

D_{x} = D_{y} = D_{z} = 0.5

Figure 5. Numerical solution of air pollutant concentration and contour plot for case 2 for

R = 0.05

and (a)

Q = 0

, (b)

Q = 0.007

, (c)

Q = - e^{- t}

and (d)

Q = 0.001 sin (x t)

Figure 5. Numerical solution of air pollutant concentration and contour plot for case 2 for

R = 0.05

and (a)

Q = 0

, (b)

Q = 0.007

, (c)

Q = - e^{- t}

and (d)

Q = 0.001 sin (x t)

Figure 6. Numerical solution of air pollutant concentration and contour plot for case 2 for

R = 0.1

and (a)

Q = 0

, (b)

Q = 0.007

, (c)

Q = - e^{- t}

and (d)

Q = 0.001 sin (x t)

Figure 6. Numerical solution of air pollutant concentration and contour plot for case 2 for

R = 0.1

and (a)

Q = 0

, (b)

Q = 0.007

, (c)

Q = - e^{- t}

and (d)

Q = 0.001 sin (x t)

Figure 7. Numerical solution of air pollutant concentration and contour plot for case 2 for

R = 0.5

and (a)

Q = 0

, (b)

Q = 0.007

, (c)

Q = - e^{- t}

and (d)

Q = 0.001 sin (x t)

Figure 7. Numerical solution of air pollutant concentration and contour plot for case 2 for

R = 0.5

and (a)

Q = 0

, (b)

Q = 0.007

, (c)

Q = - e^{- t}

and (d)

Q = 0.001 sin (x t)

Figure 8. Comparison of numerical solutions for Case 2 at

z = 0.2

m for

R = 0.05

and

Q = 0, 0.007, - e^{- t}, 0.001 sin (x t)

at (a)

y = 0.5

m and (b)

x = 0.5

Figure 8. Comparison of numerical solutions for Case 2 at

z = 0.2

m for

R = 0.05

and

Q = 0, 0.007, - e^{- t}, 0.001 sin (x t)

at (a)

y = 0.5

m and (b)

x = 0.5

Figure 9. Comparison of numerical solutions for Case 2 at

z = 0.2

m for

Q = 0.07

R = 0.05, 0.1, 0.5

at (a)

y = 0.5

m and (b)

x = 0.5

Figure 9. Comparison of numerical solutions for Case 2 at

z = 0.2

m for

Q = 0.07

R = 0.05, 0.1, 0.5

at (a)

y = 0.5

m and (b)

x = 0.5

Table 1. Boundary conditions.

Boundary	Boundary condition	Value
Entrance gate : $x = 0, 0 \leq y < 0.5, 0 \leq z \leq 1$	$C (0, y, z, t)$	0
Entrance gate : $x = 0, 0.5 \leq y \leq 1, 0 \leq z \leq 1$	$C (0, y, z, t)$	1
Exit gate : $x = 1, 0 \leq y \leq 1, 0 \leq z \leq 1$	$\frac{\partial C}{\partial x} (1, y, z, t)$	0
Right side wall : $0 \leq x < 0.3, y = 0, 0 \leq z \leq 1$	$C (x, 0, z, t)$	0
Right side wall : $0.3 \leq x \leq 0.6, y = 0, 0 \leq z \leq 1$	$C (x, 0, z, t)$	1
Right side wall : $0.6 < x \leq 1, y = 0, 0 \leq z \leq 1$	$C (x, 0, z, t)$	0
Left side wall : $0 < x < 1, y = 1, 0 \leq z \leq 1$	$\frac{\partial C}{\partial y} (x, 1, z, t)$	0
Ground : $0 < x < 1, 0 < y < 1, z = 0$	$\frac{\partial C}{\partial z} (x, y, 0, t)$	0
Ceiling : $0 < x < 1, 0 < y < 1, z = 1$	$\frac{\partial C}{\partial z} (x, y, 1, t)$	0

MDPI Initiatives

Important Links

Choose an area of interest and we will send you notifications of new preprints at your preferred frequency.

Disclaimer

Convergence Analysis of Numerical Solutions of Advection-Diffusion-Reaction Equations Using a Finite Difference Method an Application to Air Pollution Problems

Abstract

1. Introduction

2. Convergence Theory

Lax’s equivalence theorem [24]

3. Consistency and Stability Analysis

4. Numerical Results

5. Conclusions

Author Contributions

Funding

Acknowledgments

Conflicts of Interest

References

Abstract

1. Introduction

2. Convergence Theory

Lax’s equivalence theorem [24]

3. Consistency and Stability Analysis

4. Numerical Results

5. Conclusions

Author Contributions

Funding

Acknowledgments

Conflicts of Interest

References

MDPI Initiatives

Important Links

Subscribe