# Finite Volume Method for a Time-Dependent Convection-Diffusion-Reaction Equation with Small Parameters

Authors: Uzair Ahmed; Daoud Suleiman Mashat; Dalal Adnan Maturi
Journal: International Journal of Differential Equations
 (2022)

---

## Abstract

Convection, diffusion, and reaction mechanisms are characteristics of transient mass-transfer phenomena that occur in natural and industrial systems. In this article, we contemplate a passive scalar transport governed by the convection-diffusion-reaction (CDR) equation in 2D flow. The efficiency of solving computationally partial differential equations can be illustrated by using a precise numerical method that yields remarkable precision at a low cost. The accuracy and computational efficiency of two second-order finite difference methods were investigated. The results were compared to a finite volume technique, which has a memory advantage and conserves mass, momentum, and energy even on coarse grids. For various diffusion coefficient values, numerical simulation of unsteady CDR equation are also performed. The techniques were examined for consistency and convergence. The effectiveness and accuracy of these approaches for solving CDR equations are demonstrated by simulation results. Efficiency is measured usingL2 and L∞, and the estimated results are compared to the corresponding analytical solution.


---

## Body

## 1. Introduction

Different types of partial differential equations govern mathematical models of physical, chemical, biological, and environmental events. These models illustrate the conservation laws, whose solutions can be invariably discontinuous and have a large gradient, making precise numerical solutions extremely difficult to find. CDR equations appear frequently to model the phenomena in the field of biology, chemistry, ecology, physics, and engineering. The basic convection-diffusion-reaction model exhibits the time evolution of chemical or biological species. The mathematical equations representing this evolution are partial differential equations that can be derived from mass balances [1]. Convection-diffusion-reaction equations appear in modeling of the phenomenon of cancer and tumor [2–6]. Water above the ground, in pores and crevices of the ground, and the unsaturated zone are all affected by advection and diffusion, which have an impact on the transfer of contaminants. It can be significantly influenced by the reaction. The numerical features of modeling advection diffusion transport have been the subject of extensive research, which is still ongoing. As a result, a huge amount of literature are available for advection diffusion transport [7]. The mathematical models which are used to examine the pollution in water of river, stream, and lakes are also based on the convection-diffusion-reaction equation [8–10].For 2D transient scalar CDR problem, the governing partial differential equation is expressed as follows:(1)ϕt+u⟶ϕx+ϕy−ϵϕxx+ϕyy+σϕ=h.Reaction, diffusion, and convection are terms in the governing (1) containing the unknown variable’s zero, second, and first order derivatives. Electrophoresis segregation process simulation or the operation of many chemical reactors may be modeled using this type of equation as a elementary model. Both the intensity and the heating rate are unknown variables in such mechanisms. Another interesting use is the simulating fluid movement in a non-inertial environment. This event may be expressed in a diffusive-convective-reactive form, with inertial forces within the convection, viscous effects in diffusive section, and Coriolis forces in the reaction term [11]. The precise solution of the CDR equation may include abrupt strata where the gradient of the solution varies quickly whether the system is substantially governed by advection or reaction. In such instances, some traditional numerical methods produce erroneous oscillations that pollute the entire domain [12]. Because of the convection term, numerical simulation is challenging for forecasting transport processes based on transient CDR equation. The flow characteristics are frequently characterized by high gradients, necessitating particular numerical treatment. The majority of classical schemes have spurious oscillations, which result in excessive numerical dispersion [13–18].The study of CDR equation has garnered considerable attention recently [1–24]. Idelsohn et al. [11] used a Petrov–Galerkin methodology to solve CDR scenarios in 1995. Later in 1999, a semi-discretization approach and an ADI scheme, proposed by Polezhaev, were applied to inspect the CDR equation in 1D and 2D, respectively, by Sheu et al. [1]. It is possible for finite difference approximations to provide cyclic results when used to model convective-diffusive-reactive mass transfer. Akram Hossain [7] presented an explicit scheme based on finite difference approach to investigate advective-dispersive transport with reaction in 1999.Kachiashvili et al. [9] used finite difference approximations to analyze the river Khobistskali contaminated by NO3 in 1D, 2D, and 3D CDR equation. Phongthanapanich and Dechaumphai presented a hybrid finite element and finite volume model to discuss a 2D CDR equation. The CDR equation was discretized using the FVM, and the gradient variables at cell sides were defined using the finite element scheme [16]. Moreover, in 2010, Phongthanapanich and Dechaumphai [19] suggested an explicit scheme of FVM to study the CDR system with triangle based meshes. Furthermore, Sarra [6] worked on complexly shaped domains based on advection-diffusion-reaction model using a local radial basis function method in 2012.Kaya [12] also proposed a FDM to solve CDR equations. Using the self-cleansing system of the river, Simon and Koya developed a 1D Streeter–Phelps model. This model also counts the amount of pollutant downstream since the pollutants move downstream with the stream velocity. Based on a few practical presumptions, the Streeter–Phelps demonstrate is adjusted and split into a system of coupled advection-diffusion-reaction equations expressing biochemical oxygen demand and the concentrations of dissolved oxygen [8]. Subsequently, in 2017, a mathematical model describing contaminants and suspended particles in a river or waterway was examined by Samalerk and Pochai [10]. The proposed water quality model was governed by advection-diffusion-reaction equation using the Saulyev finite difference scheme. In 2018, Maghdadi et al. [4] implemented a customized cancer progression forecasting model based on CDR equation. In same year, Al-Bayati and Wrobel [20] introduced a new scheme for an unsteady convection-diffusion-reaction equation relying on the element method. Soon after, they introduced another new scheme built on the formulation of radial integration BEM in 2019 [21]. Yousefi et al. [22] devised a novel and more efficient numerical approach for solving advection-diffusion equations. The approach was based on the integration of the advection-diffusion equation. Following that, the integral equation would be converted to a system of linear algebraic equations using a HCBFs. Later, Yousefi et al. [23] used Chebyshev polynomials to solve two-dimensional transient advection-diffusion equations using the spectral integral equation approach. To approximate the solution of one- and two-dimensional conservation laws, Yousefi et al. [24] presented a higher-order scheme based on a spectral volume technique. The study revealed that the Chebyshev spectral volume approach leads to the conventional spectral volume method when used in conjunction with the novel limiter. Recently, in 2020, an unconditionally convergent implicit forward time central space scheme of FDM was applied to study the 2D advection-diffusion-reaction equation by Pananu et al. [25].The explicit and implicit approaches are two types of computational strategies for solving the PDEs. The earlier method is prevalent due to its simplicity. However, CFL stability criteria restrict the technique, so a minimal time step must be taken to maintain the stability of computation process [19]. When compared to the explicit method, the implicit method’s solution accuracy degrades with time, thus larger time steps cannot be used. It takes long time to invert the coefficient matrix, which is another issue with this problem. The amount of memory available is a critical factor to consider when simulating large-scale industrial operations.


## 2. Governing Equation

For the unsteady scalar CDR problem in two dimensions, the governing partial differential equation is expressed as follows:(2)ut+A⟶ux+uy−ϵuxx+uyy+σu=h,where x,y,t∈Γ×0T alongside initial conditions ux,y,0=u0x,y, (x, y) ∈Γ. Given the Dirichlet boundary condition can be denoted by the following expression:u (a, y, t) = f1 (x, y, t), u (b, y, t) = f2 (x, y, t)u (x, c, t) = g1 (x, y, t), u (x, d, t) = g2 (x, y, t),where (x, y, t) ∈Γ×0T, Γ=x,y:a≤x≤b,c≤y≤d in R2, time interval is 0,T . It is considered that the functions f1, f2, u0, g1, and g2 are smooth enough to simulate the conduction, diffusion, and so on.u0, f1, f2, g1, and g2 are the given sufficiently smooth functions and u (x, y, t) may represent heat, diffusion, etc. Here, h=dx=dy; diffusion coefficient is denoted by ϵ while reaction coefficient is represented by σ and A⟶ denotes the convection coefficient, respectively.


### 2.1. Analytical Solution

Analytical solution [26, 27] of the two dimensional CDR (2) is as follows:(3)Ut,x,y=16sinπtx1−xy1−y×12+arctan2ϵ−1/20.252−x−0.52−y−0.52π.Within the computations the constrain term h=dx=dy, the initial condition u0 and the boundary conditions are set from the exact solution. Here, Γ×T=x,y,t:0≤x,y≤1,0<t≤T, A⟶=2,3T and σ=1. Moreover, to see the impact of diffusion coefficient on the solution, values of diffusion coefficient is selected such that ϵ=10,10−1,10−2,10−3.


## 2.1. Analytical Solution

Analytical solution [26, 27] of the two dimensional CDR (2) is as follows:(3)Ut,x,y=16sinπtx1−xy1−y×12+arctan2ϵ−1/20.252−x−0.52−y−0.52π.Within the computations the constrain term h=dx=dy, the initial condition u0 and the boundary conditions are set from the exact solution. Here, Γ×T=x,y,t:0≤x,y≤1,0<t≤T, A⟶=2,3T and σ=1. Moreover, to see the impact of diffusion coefficient on the solution, values of diffusion coefficient is selected such that ϵ=10,10−1,10−2,10−3.


## 3. Numerical Methods

In this section, numerical schemes, to investigate the approximate solution of 2D CDR equation in a finite domainΓ, will be explored theoretically. Initially, we determine integers P and Q to express step sizes, dx=b−a/P and dy=d−c/Q in respective horizontal and vertical directions. Dividing the interval a,b into P equidistant parts of width dx and the interval c,d into Q equidistant parts of width dy. Illustrate a grid by mapping horizontal and vertical axes connecting the coordinate points xp,yq, where xp=a+p.dx for each p=0,1,2,…P and yq=c+q.dy for each q=0,1,2,…Q, also the lines x=xp and y=yq are grid lines, and their junctions are grid’s mesh. For every work point interior of the lattice, xp,yq, for p=1,2,…,P−1 and q=1,2,…,Q−1. We apply different calculations to approximate the numerical arrangement to the issue in this equation (2). We expect tn=nk, n=0,1,…,NT where t represents the time.


### 3.1. 2nd Order Crank–Nicolson Implicit Scheme

By applying the Crank–Nicolson method to (2), by integrating in the compact way, we obtain the following:(4)ut=up,qn+1−up,qnk,u=up,qn+1+up,qn2δxu=up+1,qn+1−up−1,qn+1+up+1,qn−up−1,qn4h,uy=up,q+1n+1−up,q−1n+1+up,q+1n−up,q−1n4hδx2u=up+1,qn+1−2up,qn+1+up−1,qn+1+up+1,qn−2up,qn+up−1,qn2h2δy2u=up,q+1n+1−2up,qn+1+up,q−1n+1+up,q+1n−2up,qn+up,q−1n2h2.by substituting the above expressions into (2), the 2nd order implicit method is given by the following:(5)up,qn+1−up,qn−P1up+1,qn+1−2up,qn+1+up−1,qn+1+up+1,qn−2up,qn+up−1,qn−P1up,q+1n+1−2up,qn+1+up,q−1n+1+up,q+1n−2up,qn+up,q−1n+P2up+1,qn+1−up−1,qn+1+up+1,qn−up−1,qn+P3up,q+1n+1−up,q−1n+1+up,q+1n−up,q−1n+σ2up,qn+1+up,qn=0whereP1=ε2kh2,P2=Ax4h,P3=Ay4h,The implicit scheme appears that the accuracy is ofOk2+hx2+hy2. Moreover, we use the notation Up,q and up,q representing the exact and approximate solution of mesh points xp,yq,tn respectively. Furthermore, we utilize the von Neumann stability investigation of method (4) to determine the stability of the second order implicit method, which reveals the unconditional stability of scheme. A network is presently pentadiagonal, even though results are good, but tragically because of huge iteration size, it amplifies from the inclining at slightest n sections absent in each path, but an alternative strategy is that schemes can be utilized to deal similar issues, because of the huge transfer speed, expanding network focuses the calculation ended up more troublesome. To overcome this trouble, another numerical scheme is needed due to large number of calculations. The iterative method is utilized for a linear system. The CN scheme is computationally expensive.


### 3.2. Computationally Efficient Implicit Scheme

To find a time effective substitute, we investigated that the CN implicit scheme for 2D (2) and discovered that efficiency in term of time is not a prominent feature for this scheme. To get excellent results along with time effectiveness, the efficient ADI scheme can be applied. Here, the finite difference conditions are composed at two levels. In any case, two diverse limited contrast approximations were utilized then again, initially to perform operations from n to n∗, and later from n∗ to the n+1. Similar constraints are utilized for this scheme as portrayed before. For the determination of the alternating direction implicit method, following steps have been taken:Horizontal Sweep(6)2kup,q∗−up,qn=Q1up+1,q∗−2up,q∗+up−1,q∗+up+1,qn−2up,qn+up−1,qn+Q2up,q+1n−2up,qn+up,q−1n−Q3up+1,qn−up−1,qn−Q4up,q+1n−up,q−1n−σup,qn.In simplified form, the (6) is as follows:(7)−Q1up−1,q∗+2k+Q2up,q∗−Q1up+1,q∗=Q1up+1,qn−2up,qn+up−1,qn+Q2up,q+1n−2up,qn+up,q−1n−Q3up+1,qn−up−1,qn−Q4up,q+1n−up,q−1n+Gun.Initial Node in Horizontal Direction(8)2k+Q2up,q∗−Q1up+1,q∗=Q1up+1,qn−2up,qn+up−1,qn+Q2up,q+1n−2up,qn+up,q−1n−Q3up+1,qn−up−1,qn−Q4up,q+1n−up,q−1n+Gun+Q1up−1,q∗.Final Node in Horizontal Direction(9)−Q1up−1,q∗+2k+Q2up,q∗=Q1up+1,qn−2up,qn+up−1,qn+Q2up,q+1n−2up,qn+up,q−1n−Q3up+1,qn−up−1,qn−Q4up,q+1n−up,q−1n+Gun+Q1up+1,q∗, where Q1=Γ/2h2,Q2=Γ/h2,Q3=Ax/2h,Q4=Ay/2handGun=2/k−σup,qnSweep in Vertical Direction(10)2kup,qn+1−up,q∗=Q1up,q+1n+1−2up,qn+1+up,q−1n+1+up,q−1∗−2up,q∗+up,q−1∗+Q2up+1,q∗−2up,q∗+up−1,q∗−Q3up+1,q∗−up−1,q∗−Q4up,q+1∗−up,q−1∗−σup,q∗.In simplified form the (10) is as follows:(11)−Q1up,q−1n+1+2k+Q2up,qn+1−Q1up,q+1n+1=Q1up,q+1∗−2up,q∗+up,q−1∗+Q2up+1,q∗−2up,q∗+up−1,q∗−Q3up+1,q∗−up−1,q∗−Q4up,q+1∗−up,q−1∗+2k−σup,q∗.Initial Node in Vertical Direction(12)2k+Q2up,qn+1−Q1up,q+1n+1=Q1up,q+1∗−2up,q∗+up,q−1∗+Q2up+1,q∗−2up,q∗+up−1,q∗−Q3up+1,q∗−up−1,q∗−Q4up,q+1∗−up,q−1∗+2k−σup,q∗+Q1up,q−1n+1.Final Node in Vertical Direction(13)−Q1up,q−1n+1+2k+Q2up,qn+1=Q1up,q+1∗−2up,q∗+up,q−1∗+Q2up+1,q∗−2up,q∗+up−1,q∗−Q3up+1,q∗−up−1,q∗−Q4up,q+1∗−up,q−1∗+2k−σup,q∗+Q1up,q+1n+1, where Q1=Γ/2h2,Q2=Γ/h2,Q3=Ax/2h,Q4=Ay/2handGu∗=2/k−σup,q∗This scheme has accuracy of 2nd order both in space and time along with the unconditional stability.The linear system in the horizontal direction is as follows:(14)a1ux,mu∗=b1m+Fkux,mn,where m = 1, …, M − 1The linear system in the vertical direction is as follows:(15)c1uy,mn+1=d1m+Gkuy,m∗,where l = 1, …, L − 1wherea1,c1 create tridiagonal matrix and the array b1,d1 are dependent upon p and qThe reactive term in horizontal-directionGkux,mn=kGu1,mn,…,kGu1,L−1nSo, also for the reactive term in the vertical direction, it is as follows:Gkuy,l∗=kGu1,M−1∗,…,kGu1,1∗The notationsup,q∗andup,qn+1 are used in same manner. The scheme gives us 2nd order accuracy in time and space Ok2+hx2+hy2. For stability analysis, this scheme exhibits unconditional stability. The stability of alternating direction implicit strategy (6) and (10) can be obtained by using the von Neumann criteria for stability. The regular stability investigation criteria indicate toward the amplification factor G of (6) and (10) may be observed as by substituting:(16)un=λneikxrhx+kyshy.It is obvious that the amplification factor is always less of equal toG≤1. Hence, it has unconditional stability. Eventually, the methodology results in a tridiagonal linear system. Iterative strategies were applied to unravel the abovementioned framework. The trick utilized in developing the ADI scheme is to divide the time step into two, and apply two diverse stencils in each half time step; hence, to increase time by one time step in grid point. Initially, we determine the values of these stencils, which are selected so that the arising linear framework is tridiagonal [28–32]. To get the approximate solution, we have to find the solution of a non-linear tridiagonal framework with each passing step of time.


### 3.3. Finite Volume Method

The finite volume method has a lot of flexibility in solving the transport phenomena equations [33–39]. Whenever fluid flow appears, finite volume becomes more and more useful. When we have a fluid flow, we must have the corresponding velocity field satisfying the continuity equation. By solving (2) using the finite volume method, we define the quantities as algebraic values of the variables u at the main grid points s.t. (W, E, P, N, S). To do that, we need profile assumptions. For applying the finite volume method, following steps should be needed.where(17)W=west,E=east,P=center point,N=north,S=south.First, divide the domain logically into a number of finite size sub-domains. Each sub-domains is represented by a finite number of grid points. Subsequently, it converts and integrates the partial differential equation into balancing equations for a single element. Using an integrated quadrature of a given order of precision, the process converts volumetric and area integrals to algebraic variables across components and their surfaces. As a consequence, a collection of semi-discrete equations arises. Also, at this step, interpolation profiles are defined to estimate the fluctuation of the variables within the element and link the surface variables to their data points, thereby transforming algebraic connections to equations. In respect to the two steps within the FVM process, approximation selection influences the by and large exactness of the consequent numeric.


#### 3.3.1. 2nd Order Finite Volume Scheme

By solving (2) using the finite volume method, we have the following:(18)uPn+1=C1uPn+C2DeuEn−uPn−DwuPn−uwn+C3DnuNn−uPn−DsuPn−uSn−C4uEn−uWn−C5uNn−uSn,where(19)C1=1−σΔxΔy,C2=ΔtΔx,C3=ΔtΔy,C4=ΔtAx2Δx,C5=ΔtAy2Δy,De=ϵΔxe,Dw=ϵΔxw,Dn=ϵΔyn,Ds=ϵΔys,Bottom left corner boundary cell(20)uPn+1=C1uPn+C2DeuEn−uPn−2DwuPn−uleft+C3DnuNn−uPn−2DsuPn−ubottom−C4uEn−uleft−C5uNn−ubottom,Bottom right corner boundary cell(21)uPn+1=C1uPn+C22Deuright−uPn−DwuPn−uwn+C3DnuNn−uPn−2DsuPn−ubottom−C4uright−uWn−C5uNn−ubottom,Top left corner boundary cell(22)uPn+1=C1uPn+C2DeuEn−uPn−2DwuPn−uleft+C32Dnutop−uPn−DsuPn−uSn−C4uEn−uleft−C5utop−uSn,Top right corner boundary cell(23)uPn+1=C1uPn+C22Deuright−uPn−DwuPn−uwn+C32Dnutop−uPn−DsuPn−uSn−C4uright−uWn−C5utop−uSn,Bottom non-corner boundary cell(24)uPn+1=C1uPn+C2DeuEn−uPn−DwuPn−uwn+C3DnuNn−uPn−2DsuPn−ubottom−C4uEn−uWn−C5uNn−ubottom,Left non-corner boundary cell(25)uPn+1=C1uPn+C2DeuEn−uPn−2DwuPn−uleft+C3DnuNn−uPn−DsuPn−uSn−C4uEn−uleft−C5uNn−uSn,Top non-corner boundary cell(26)uPn+1=C1uPn+C2DeuEn−uPn−DwuPn−uwn+C32Dnutop−uPn−DsuPn−uSn−C4uEn−uWn−C5utop−uSn,Right non-corner boundary cell(27)uPn+1=C1uPn+C22Deuright−uPn−DwuPn−uwn+C3DnuNn−uPn−DsuPn−uSn−C4uright−uWn−C5uNn−uSn.The scheme gives us an accuracy of second order. For stability criteria, we know thatPe=advection strength/diffusion strength Figure 1. To define the stability of 2nd order finite volume scheme, Pe must be less than or equal to 2. The method in (18) demonstrates that despite its conservativeness and convergence, it is conditionally stable. Since the abovementioned inequality establishes a straight maximum limit for the limit step, improving spatial precision can be costly at times.Figure 1
Peclet number.


### 3.4. Algorithm

Clearly, the framework is tridiagonal and can be fathomed with the Thomas algorithm. The measurement of H isp×q. In common a tridiagonal framework can be composed as follows:(28)apxp−1+bpxp+cpxp+1=Sp,witha1=cl=0,The above system can be expressed as follows in a matrix-vector form:(29)Hu=S,where H represents a coefficient matrix. The column vector on the right side is known. Our primary objective is to discover the resultant vector u. Presently, we have the following:(30)H=b1c1000⋯0a2b2c200⋯0⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯an−1bn−1cn−1⋯⋯⋯⋯anbn,u¯=u1,u2,u3,…,untS¯=s1,s2,s3,…,snt,The method is clarified within the taking after steps,(31)H=LU,L=μ10000…0β2μ2000…0α3β3μ300…0………………0……αn−1βn−1μn−10………αnβnμn,and(32)U=1δ1λ100…001δ2λ20…0001δ3λ3…0………………0……01δn−2λ2……001δn−1………001.By comparing both sides of theHu=S, it gives us the components of the matrices L and U. The computational steps for the execution of the Thomas calculation are appeared in comes about, taken from particular illustrations.


## 3.1. 2nd Order Crank–Nicolson Implicit Scheme

By applying the Crank–Nicolson method to (2), by integrating in the compact way, we obtain the following:(4)ut=up,qn+1−up,qnk,u=up,qn+1+up,qn2δxu=up+1,qn+1−up−1,qn+1+up+1,qn−up−1,qn4h,uy=up,q+1n+1−up,q−1n+1+up,q+1n−up,q−1n4hδx2u=up+1,qn+1−2up,qn+1+up−1,qn+1+up+1,qn−2up,qn+up−1,qn2h2δy2u=up,q+1n+1−2up,qn+1+up,q−1n+1+up,q+1n−2up,qn+up,q−1n2h2.by substituting the above expressions into (2), the 2nd order implicit method is given by the following:(5)up,qn+1−up,qn−P1up+1,qn+1−2up,qn+1+up−1,qn+1+up+1,qn−2up,qn+up−1,qn−P1up,q+1n+1−2up,qn+1+up,q−1n+1+up,q+1n−2up,qn+up,q−1n+P2up+1,qn+1−up−1,qn+1+up+1,qn−up−1,qn+P3up,q+1n+1−up,q−1n+1+up,q+1n−up,q−1n+σ2up,qn+1+up,qn=0whereP1=ε2kh2,P2=Ax4h,P3=Ay4h,The implicit scheme appears that the accuracy is ofOk2+hx2+hy2. Moreover, we use the notation Up,q and up,q representing the exact and approximate solution of mesh points xp,yq,tn respectively. Furthermore, we utilize the von Neumann stability investigation of method (4) to determine the stability of the second order implicit method, which reveals the unconditional stability of scheme. A network is presently pentadiagonal, even though results are good, but tragically because of huge iteration size, it amplifies from the inclining at slightest n sections absent in each path, but an alternative strategy is that schemes can be utilized to deal similar issues, because of the huge transfer speed, expanding network focuses the calculation ended up more troublesome. To overcome this trouble, another numerical scheme is needed due to large number of calculations. The iterative method is utilized for a linear system. The CN scheme is computationally expensive.


## 3.2. Computationally Efficient Implicit Scheme

To find a time effective substitute, we investigated that the CN implicit scheme for 2D (2) and discovered that efficiency in term of time is not a prominent feature for this scheme. To get excellent results along with time effectiveness, the efficient ADI scheme can be applied. Here, the finite difference conditions are composed at two levels. In any case, two diverse limited contrast approximations were utilized then again, initially to perform operations from n to n∗, and later from n∗ to the n+1. Similar constraints are utilized for this scheme as portrayed before. For the determination of the alternating direction implicit method, following steps have been taken:Horizontal Sweep(6)2kup,q∗−up,qn=Q1up+1,q∗−2up,q∗+up−1,q∗+up+1,qn−2up,qn+up−1,qn+Q2up,q+1n−2up,qn+up,q−1n−Q3up+1,qn−up−1,qn−Q4up,q+1n−up,q−1n−σup,qn.In simplified form, the (6) is as follows:(7)−Q1up−1,q∗+2k+Q2up,q∗−Q1up+1,q∗=Q1up+1,qn−2up,qn+up−1,qn+Q2up,q+1n−2up,qn+up,q−1n−Q3up+1,qn−up−1,qn−Q4up,q+1n−up,q−1n+Gun.Initial Node in Horizontal Direction(8)2k+Q2up,q∗−Q1up+1,q∗=Q1up+1,qn−2up,qn+up−1,qn+Q2up,q+1n−2up,qn+up,q−1n−Q3up+1,qn−up−1,qn−Q4up,q+1n−up,q−1n+Gun+Q1up−1,q∗.Final Node in Horizontal Direction(9)−Q1up−1,q∗+2k+Q2up,q∗=Q1up+1,qn−2up,qn+up−1,qn+Q2up,q+1n−2up,qn+up,q−1n−Q3up+1,qn−up−1,qn−Q4up,q+1n−up,q−1n+Gun+Q1up+1,q∗, where Q1=Γ/2h2,Q2=Γ/h2,Q3=Ax/2h,Q4=Ay/2handGun=2/k−σup,qnSweep in Vertical Direction(10)2kup,qn+1−up,q∗=Q1up,q+1n+1−2up,qn+1+up,q−1n+1+up,q−1∗−2up,q∗+up,q−1∗+Q2up+1,q∗−2up,q∗+up−1,q∗−Q3up+1,q∗−up−1,q∗−Q4up,q+1∗−up,q−1∗−σup,q∗.In simplified form the (10) is as follows:(11)−Q1up,q−1n+1+2k+Q2up,qn+1−Q1up,q+1n+1=Q1up,q+1∗−2up,q∗+up,q−1∗+Q2up+1,q∗−2up,q∗+up−1,q∗−Q3up+1,q∗−up−1,q∗−Q4up,q+1∗−up,q−1∗+2k−σup,q∗.Initial Node in Vertical Direction(12)2k+Q2up,qn+1−Q1up,q+1n+1=Q1up,q+1∗−2up,q∗+up,q−1∗+Q2up+1,q∗−2up,q∗+up−1,q∗−Q3up+1,q∗−up−1,q∗−Q4up,q+1∗−up,q−1∗+2k−σup,q∗+Q1up,q−1n+1.Final Node in Vertical Direction(13)−Q1up,q−1n+1+2k+Q2up,qn+1=Q1up,q+1∗−2up,q∗+up,q−1∗+Q2up+1,q∗−2up,q∗+up−1,q∗−Q3up+1,q∗−up−1,q∗−Q4up,q+1∗−up,q−1∗+2k−σup,q∗+Q1up,q+1n+1, where Q1=Γ/2h2,Q2=Γ/h2,Q3=Ax/2h,Q4=Ay/2handGu∗=2/k−σup,q∗This scheme has accuracy of 2nd order both in space and time along with the unconditional stability.The linear system in the horizontal direction is as follows:(14)a1ux,mu∗=b1m+Fkux,mn,where m = 1, …, M − 1The linear system in the vertical direction is as follows:(15)c1uy,mn+1=d1m+Gkuy,m∗,where l = 1, …, L − 1wherea1,c1 create tridiagonal matrix and the array b1,d1 are dependent upon p and qThe reactive term in horizontal-directionGkux,mn=kGu1,mn,…,kGu1,L−1nSo, also for the reactive term in the vertical direction, it is as follows:Gkuy,l∗=kGu1,M−1∗,…,kGu1,1∗The notationsup,q∗andup,qn+1 are used in same manner. The scheme gives us 2nd order accuracy in time and space Ok2+hx2+hy2. For stability analysis, this scheme exhibits unconditional stability. The stability of alternating direction implicit strategy (6) and (10) can be obtained by using the von Neumann criteria for stability. The regular stability investigation criteria indicate toward the amplification factor G of (6) and (10) may be observed as by substituting:(16)un=λneikxrhx+kyshy.It is obvious that the amplification factor is always less of equal toG≤1. Hence, it has unconditional stability. Eventually, the methodology results in a tridiagonal linear system. Iterative strategies were applied to unravel the abovementioned framework. The trick utilized in developing the ADI scheme is to divide the time step into two, and apply two diverse stencils in each half time step; hence, to increase time by one time step in grid point. Initially, we determine the values of these stencils, which are selected so that the arising linear framework is tridiagonal [28–32]. To get the approximate solution, we have to find the solution of a non-linear tridiagonal framework with each passing step of time.


## 3.3. Finite Volume Method

The finite volume method has a lot of flexibility in solving the transport phenomena equations [33–39]. Whenever fluid flow appears, finite volume becomes more and more useful. When we have a fluid flow, we must have the corresponding velocity field satisfying the continuity equation. By solving (2) using the finite volume method, we define the quantities as algebraic values of the variables u at the main grid points s.t. (W, E, P, N, S). To do that, we need profile assumptions. For applying the finite volume method, following steps should be needed.where(17)W=west,E=east,P=center point,N=north,S=south.First, divide the domain logically into a number of finite size sub-domains. Each sub-domains is represented by a finite number of grid points. Subsequently, it converts and integrates the partial differential equation into balancing equations for a single element. Using an integrated quadrature of a given order of precision, the process converts volumetric and area integrals to algebraic variables across components and their surfaces. As a consequence, a collection of semi-discrete equations arises. Also, at this step, interpolation profiles are defined to estimate the fluctuation of the variables within the element and link the surface variables to their data points, thereby transforming algebraic connections to equations. In respect to the two steps within the FVM process, approximation selection influences the by and large exactness of the consequent numeric.


### 3.3.1. 2nd Order Finite Volume Scheme

By solving (2) using the finite volume method, we have the following:(18)uPn+1=C1uPn+C2DeuEn−uPn−DwuPn−uwn+C3DnuNn−uPn−DsuPn−uSn−C4uEn−uWn−C5uNn−uSn,where(19)C1=1−σΔxΔy,C2=ΔtΔx,C3=ΔtΔy,C4=ΔtAx2Δx,C5=ΔtAy2Δy,De=ϵΔxe,Dw=ϵΔxw,Dn=ϵΔyn,Ds=ϵΔys,Bottom left corner boundary cell(20)uPn+1=C1uPn+C2DeuEn−uPn−2DwuPn−uleft+C3DnuNn−uPn−2DsuPn−ubottom−C4uEn−uleft−C5uNn−ubottom,Bottom right corner boundary cell(21)uPn+1=C1uPn+C22Deuright−uPn−DwuPn−uwn+C3DnuNn−uPn−2DsuPn−ubottom−C4uright−uWn−C5uNn−ubottom,Top left corner boundary cell(22)uPn+1=C1uPn+C2DeuEn−uPn−2DwuPn−uleft+C32Dnutop−uPn−DsuPn−uSn−C4uEn−uleft−C5utop−uSn,Top right corner boundary cell(23)uPn+1=C1uPn+C22Deuright−uPn−DwuPn−uwn+C32Dnutop−uPn−DsuPn−uSn−C4uright−uWn−C5utop−uSn,Bottom non-corner boundary cell(24)uPn+1=C1uPn+C2DeuEn−uPn−DwuPn−uwn+C3DnuNn−uPn−2DsuPn−ubottom−C4uEn−uWn−C5uNn−ubottom,Left non-corner boundary cell(25)uPn+1=C1uPn+C2DeuEn−uPn−2DwuPn−uleft+C3DnuNn−uPn−DsuPn−uSn−C4uEn−uleft−C5uNn−uSn,Top non-corner boundary cell(26)uPn+1=C1uPn+C2DeuEn−uPn−DwuPn−uwn+C32Dnutop−uPn−DsuPn−uSn−C4uEn−uWn−C5utop−uSn,Right non-corner boundary cell(27)uPn+1=C1uPn+C22Deuright−uPn−DwuPn−uwn+C3DnuNn−uPn−DsuPn−uSn−C4uright−uWn−C5uNn−uSn.The scheme gives us an accuracy of second order. For stability criteria, we know thatPe=advection strength/diffusion strength Figure 1. To define the stability of 2nd order finite volume scheme, Pe must be less than or equal to 2. The method in (18) demonstrates that despite its conservativeness and convergence, it is conditionally stable. Since the abovementioned inequality establishes a straight maximum limit for the limit step, improving spatial precision can be costly at times.Figure 1
Peclet number.


## 3.3.1. 2nd Order Finite Volume Scheme

By solving (2) using the finite volume method, we have the following:(18)uPn+1=C1uPn+C2DeuEn−uPn−DwuPn−uwn+C3DnuNn−uPn−DsuPn−uSn−C4uEn−uWn−C5uNn−uSn,where(19)C1=1−σΔxΔy,C2=ΔtΔx,C3=ΔtΔy,C4=ΔtAx2Δx,C5=ΔtAy2Δy,De=ϵΔxe,Dw=ϵΔxw,Dn=ϵΔyn,Ds=ϵΔys,Bottom left corner boundary cell(20)uPn+1=C1uPn+C2DeuEn−uPn−2DwuPn−uleft+C3DnuNn−uPn−2DsuPn−ubottom−C4uEn−uleft−C5uNn−ubottom,Bottom right corner boundary cell(21)uPn+1=C1uPn+C22Deuright−uPn−DwuPn−uwn+C3DnuNn−uPn−2DsuPn−ubottom−C4uright−uWn−C5uNn−ubottom,Top left corner boundary cell(22)uPn+1=C1uPn+C2DeuEn−uPn−2DwuPn−uleft+C32Dnutop−uPn−DsuPn−uSn−C4uEn−uleft−C5utop−uSn,Top right corner boundary cell(23)uPn+1=C1uPn+C22Deuright−uPn−DwuPn−uwn+C32Dnutop−uPn−DsuPn−uSn−C4uright−uWn−C5utop−uSn,Bottom non-corner boundary cell(24)uPn+1=C1uPn+C2DeuEn−uPn−DwuPn−uwn+C3DnuNn−uPn−2DsuPn−ubottom−C4uEn−uWn−C5uNn−ubottom,Left non-corner boundary cell(25)uPn+1=C1uPn+C2DeuEn−uPn−2DwuPn−uleft+C3DnuNn−uPn−DsuPn−uSn−C4uEn−uleft−C5uNn−uSn,Top non-corner boundary cell(26)uPn+1=C1uPn+C2DeuEn−uPn−DwuPn−uwn+C32Dnutop−uPn−DsuPn−uSn−C4uEn−uWn−C5utop−uSn,Right non-corner boundary cell(27)uPn+1=C1uPn+C22Deuright−uPn−DwuPn−uwn+C3DnuNn−uPn−DsuPn−uSn−C4uright−uWn−C5uNn−uSn.The scheme gives us an accuracy of second order. For stability criteria, we know thatPe=advection strength/diffusion strength Figure 1. To define the stability of 2nd order finite volume scheme, Pe must be less than or equal to 2. The method in (18) demonstrates that despite its conservativeness and convergence, it is conditionally stable. Since the abovementioned inequality establishes a straight maximum limit for the limit step, improving spatial precision can be costly at times.Figure 1
Peclet number.


## 3.4. Algorithm

Clearly, the framework is tridiagonal and can be fathomed with the Thomas algorithm. The measurement of H isp×q. In common a tridiagonal framework can be composed as follows:(28)apxp−1+bpxp+cpxp+1=Sp,witha1=cl=0,The above system can be expressed as follows in a matrix-vector form:(29)Hu=S,where H represents a coefficient matrix. The column vector on the right side is known. Our primary objective is to discover the resultant vector u. Presently, we have the following:(30)H=b1c1000⋯0a2b2c200⋯0⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯⋯an−1bn−1cn−1⋯⋯⋯⋯anbn,u¯=u1,u2,u3,…,untS¯=s1,s2,s3,…,snt,The method is clarified within the taking after steps,(31)H=LU,L=μ10000…0β2μ2000…0α3β3μ300…0………………0……αn−1βn−1μn−10………αnβnμn,and(32)U=1δ1λ100…001δ2λ20…0001δ3λ3…0………………0……01δn−2λ2……001δn−1………001.By comparing both sides of theHu=S, it gives us the components of the matrices L and U. The computational steps for the execution of the Thomas calculation are appeared in comes about, taken from particular illustrations.


## 4. Error Norms

The precision and consistency of the methods is studied in terms of error norms specificallyL2 and L∞ which can be expressed as follows:(33)RMS Error=∑p,q=1LUp,q−up,q2L×L,L∞=max1≤p≤P∑q=1PUp.q−up,q,L2=ρUp,q−up,qtUp,q−up,q,where u (x, y, t) and U (x, y, t) represent the approximate and exact solutions at the grid point xp,yq,tn. In this method, ρUp,q−up,q=maxλ and λ is an eigen value of Up,q−up,q.


## 5. Results and Discussion

In several application areas, such as electromagnetics, biomathematics, precise modeling of electric signals, and water modeling takes more time. Finite volume and second order implicit schemes have been devised to solve this problem and to lower the error norm that best marches the accuracy effectively. To calculate numerical derivatives effectively, we utilized a short step size along the grid line. MATLAB algorithms for the 2D CDR (2) were also constructed.The ADI scheme is used to reduce computing costs while increasing efficiency, and the FVM is applied to ensure conservation over the entire scope of the problem. The convective and diffusive fluxes at the control volume faces are calculated using a finite volume approach that employs consistent formulations. The article’s main goal is to identify the accuracy and efficiency of second order numerical methods. The algorithms developed are compared to (3) analytical findings. The parameters σ=1, A⟶=2,3T, and changing ϵ=10,10−1,10−2 for the 2D CDR equation are used to achieve this goal. To emphasize the relevance of second order schemes on such phenomenological systems, numerical and analytical solutions are contrasted and justified in terms of error norms.The calculated results are realized by changing the grid space, time space and by varyingϵ values. Furthermore, utilizing 2nd-order finite difference schemes, the comparison of analytic and numerical results with ϵ=0.01 linked to two variables p,q at t=0.1 with a fixed grid size of 10×10 is found in Table 1. The truncation error is computed using the L2 and L∞ norms with a fixed grid size of 10×10 in Table 2. Table 2 demonstrates that modifying the time t=2,1,0.5 with the same grid size at k=0.0001,epsilon=0.001, improves the accuracy of the L∞ error norm. Table 3 shows that error decreases as grid size increases at k=0.001. We also discovered that the alternating direction implicit method is more time efficient. In comparison to the Crank–Nicolson method, the findings show that ADI has accurate and excellent convergence. It is evident that as ϵ decreases, (2) becomes convection dominated and as a result, we’ll require extra components to deal with the associated inconsistencies. Table 4 shows that if treated correctly, FVM may produce results that are comparable to the analytical solution for higher and lower values of ϵ. On the other hand, all these approaches have computational problems for smaller values of ϵ≤10−3.Table 1
Evaluation of approximate and exact results of ADI and CN att = 0.1 for unknown up,q,t. Numerical solution at ϵ=0.1,k=0.001,gridsize=10×10.
Solution comparisonExactCNADIError-CNError-ADI1.5996−×10−22.8032×10−22.9011×10−2−1.2036×10−2−1.3015×10−2−1.8177×10−2−1.7837×10s−1.9226×10−31.9961×10−2−1.6393×10−2−9.8860×10−3−8.3268×10−3−8.1314×10−3−1.5592×10−3−1.7546×10−31.1234×10−21.0062×10−21.2385×10−21.1720×10−3−1.1510×10−31.8177×10−21.7006×10−21.7231×10−21.1710×10−38.5600×10−4Table 2
Error profile data by using CN & ADI att = 2 for unknown up,q,t.
Error estimation atϵ=0.01,k=0.0001,gridsize=10×10.tL2−CNL∞−CNL2−ADIL∞−ADI24.6192×10−17.1484×10−24.1595×10−17.8537×10−211.3019×10−21.2378×10−21.4931×10−21.2681×10−20.51.0178×10−22.0701×10−31.1231×10−22.0371×10−3Table 3
Errors of CN and ADI at different gird sizes withk=0.001 for unknown up,q,t.
Error estimation atϵ=0.01.Grid-SizeLinf−CNLinf−ADITotal-time CNTotal-time ADISelf-time CNSelf-time ADI5×56.7949×10−36.8337×10−30.130.050.0020.000310×106.5020×10−36.5522×10−30.1740.1090.0550.00520×206.7281×10−46.6720×10−456.91533.190.1860.09525×256.7949×10−56.7312×10−5134.680.30.2620.159Table 4
Error norm of finite volume by varyingϵ for unknown up,q,t.
Error estimation at grid size=20×20ϵL2L∞RMS103.4721×10−35.4204×10−46.6410×10−311.2739×10−29.7252×10−35.4026×10−210−14.0424×10−13.6456×10−21.4912×10−210−21.4352×10−15.4240×10−15.4212×10−1Furthermore, it is observed that by reducing step sizesh values in Table 5 with grid size 20×20andϵ=0.001 can prevent turbulence. In this situation, applying error norms to reduce the space step size improves the accuracy of the finite volume technique findings. By decreasing the space step, the root mean square value contributes to accuracy. It is a great strategy for dealing with conservative rules. In addition, the comparison of three techniques will be taken into account in Table 6. The spatial step size h and time level t=0.1,0.001,0.5 with a fixed grid size of 50×50 is presented using error analysis norm. Results can be seen using the L∞ norm. Variation of the h parameter toward smaller values results in an increase in accuracy. When we changed the value of t from high to low, we saw that the error dropped. Moreover, the use of a finite volume method yielded conservative findings that fully support the precise solution. The CN and ADI schemes solve the problem effectively because of this excellent characteristic of the solution, but the second order finite volume scheme also shows a compatible solution with exactitude. In Table 7, CPU time shows efficiency of the proposed scheme with respect to self-time, or time invested in its child unit relative to the algorithm’s overall run time, which is the time spent in its child function to the total time spent in the algorithm. When the two machines are compared by changing the grid size, it is clear that finite volume produces good results with tiny intervals. Using two machines, HP and Dell, it was discovered that Dell is more efficient than HP.Table 5
Error norm of finite volume at different step sizes for unknownup,q,t.
Numerical solution of FVM at grid size=20×20 and ϵ=0.01HL2L∞RMS0.12.4721×10−22.2954×10−22.4721×10−10.011.5709×10−29.9821×10−41.5709×10−20.0054.5319×10−32.7512×10−44.3512×10−30.0011.8429×10−41.6710×10−51.8429×10−40.00011.8657×10−61.1801×10−71.8657×10−6Table 6
Error norms of CN, ADI, and FVM at different times and space step size for unknownup,q,t.
Error estimation at grid size=50×50Time levelhL∞−CNL∞−ADIL∞−FVM0.010.0013.5078×10−43.4638×10−43.4931×10−40.10.015.8713×10−35.7141×10−35.8136×10−30.50.15.8178×10−25.2231×10−25.5373×10−2Table 7
Performance of CPU time efficiency using FVM.
Grid sizeDellHPTotal-timeSelf-timeTotal-timeSelf-time25×2599.2180.175113.8210.18850×50224.5010.206284.8660.22975×75486.1600.346640.0870.390100×1002308.5260.8053010.2641.009On a two-dimensional partial differential system, the graphical representation of numerical schemes is seen. A comparison was performed between analytic and numerical findings using a second order finite difference and finite volume methodology. Fork=0.0001, t=2, ϵ=0.1 and grid size =10×10 can be seen from Figure 2. A study of CN’s approximate and analytical findings will be taken into account in this figure. The approximate results of the CN and ADI schemes are shown in Figures 3 and 4 for k=0.001 and k=0.0001, respectively. Figure 5 shows more accurate and refined observations than Figure 2 while using the same time step. An example is shown in Figure 6, which compares the exact and numerical solution of second order finite volume scheme at time ϵ=0.1 and grid size =20×20. A detailed examination of Figure 6 reveals good agreement with the exact answer, as evidenced by the data points at various places. Figure 7 compares the exact and numerical findings of a second order finite volume method at h=0.2 with smaller time space. Figure 8 shows efficient and reliable numerical results using the ADI technique with grid size =75×75, ϵ=0.01 and k=0.001. At a higher grid size, the numerical results of the FVM are contrasted to the exact values in Figure 9. While, two-dimensional numerical methods at different time and epsilon along the p,q axis may be seen in all figures. It demonstrates that the degree of numerical solutions in different grid places approximates the analytic results. The findings demonstrate that these schemes are reliable for solving PDEs.Figure 2
Analytical and approximate results of the CN scheme at10×10 grid size, time =2, time step =0.0001, and ϵ=0.1.Figure 3
Numerical solution of the CN scheme at=30×30 grid size, k=0.01, and ϵ=0.01.Figure 4
Approximate results of ADI atk=0.001 and grid size =25×25.Figure 5
Analytical and approximate results of ADI at15×15 grid size, time =0.5, and time step =0.0001.Figure 6
Analytical and approximate results of FVM at20×20 grid size, time step =0.01.ϵ=0.1.Figure 7
Analytical and approximate results of FVM at15×15 grid size, time step =0.0001, and h=0.2.Figure 8
Analytical and approximate results of ADI at75×75 grid size, time step =0.001, and ϵ=0.01.Figure 9
Analytical and approximate results of FVM at100×100 grid size, time step =0.01, and ϵ=0.1.In conclusion, the numerical experiments demonstrated the correctness and effectiveness of the numerical methods. However, the finite volume scheme is more suitable because it outperforms the finite difference method in the way that it does not require a structured mesh, but it can be used if desired. Despite the fact that structured grids require less memory to store than unstructured grids due to their easier connection, FVM surpass the finite difference methods for extremely complicated geometries. For the stated finite volume formulation of the diffusive flow, flux consistency assures conservation throughout the whole domain. It is aspirant that the analytical and numerical solutions are strongly aligned with generation encrypting, as shown in the figures and tables.


## 6. Conclusion

The 2D transient scalar CDR problem is solved numerically in this paper. Theoretical formulations of the proposed schemes were fully addressed. The three scheme’s efficiency, precision, and stability are shown by a convincing example. Increasing the grid size or reducing the space step size resulted in error norms convergent toward zero. Although all of the schemes are consistent, the alternating direction implicit and finite volume schemes are less computationally expensive than the Crank–Nicolson (CN) scheme. The FVM, on the other hand, ensures that conservation is maintained across the domain. The stated techniques might be used to address similar mathematical challenges in biology, physics, and engineering.



---

*Source: 3476309-2022-05-17.xml*