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

Convection, diusion, 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-diusion-reaction (CDR) equation in 2D ow. e eciency of solving computationally partial dierential equations can be illustrated by using a precise numerical method that yields remarkable precision at a low cost.e accuracy and computational eciency of two secondorder nite dierence methods were investigated. e results were compared to a nite volume technique, which has a memory advantage and conserves mass, momentum, and energy even on coarse grids. For various diusion coecient values, numerical simulation of unsteady CDR equation are also performed. e techniques were examined for consistency and convergence. e eectiveness and accuracy of these approaches for solving CDR equations are demonstrated by simulation results. Eciency is measured using L2 and L∞, and the estimated results are compared to the corresponding analytical solution.


Introduction
Di erent types of partial di erential equations govern mathematical models of physical, chemical, biological, and environmental events. ese models illustrate the conservation laws, whose solutions can be invariably discontinuous and have a large gradient, making precise numerical solutions extremely di cult to nd. CDR equations appear frequently to model the phenomena in the eld of biology, chemistry, ecology, physics, and engineering. e basic convection-di usion-reaction model exhibits the time evolution of chemical or biological species. e mathematical equations representing this evolution are partial di erential equations that can be derived from mass balances [1]. Convection-di usion-reaction equations appear in modeling of the phenomenon of cancer and tumor [2][3][4][5][6]. Water above the ground, in pores and crevices of the ground, and the unsaturated zone are all a ected by advection and di usion, which have an impact on the transfer of contaminants. It can be signi cantly in uenced by the reaction. e numerical features of modeling advection di usion transport have been the subject of extensive research, which is still ongoing. As a result, a huge amount of literature are available for advection di usion transport [7]. e mathematical models which are used to examine the pollution in water of river, stream, and lakes are also based on the convection-di usion-reaction equation [8][9][10].
For 2D transient scalar CDR problem, the governing partial di erential equation is expressed as follows: Reaction, di usion, and convection are terms in the governing (1) containing the unknown variable's zero, second, and rst 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 uid movement in a noninertial environment.
is event may be expressed in a di usive-convective-reactive form, with inertial forces within the convection, viscous e ects in di usive section, and Coriolis forces in the reaction term [11]. e 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. e flow characteristics are frequently characterized by high gradients, necessitating particular numerical treatment. e majority of classical schemes have spurious oscillations, which result in excessive numerical dispersion [13][14][15][16][17][18].
e study of CDR equation has garnered considerable attention recently . 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 NO 3 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. e 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. is 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]. e 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. e 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. e 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]. e explicit and implicit approaches are two types of computational strategies for solving the PDEs. e 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. e amount of memory available is a critical factor to consider when simulating large-scale industrial operations.

Governing Equation
For the unsteady scalar CDR problem in two dimensions, the governing partial differential equation is expressed as follows: where (x, y, t) ∈ Γ × (0 T] alongside initial conditions u(x, y, 0) � u0(x, y), (x, y) ∈Γ. Given the Dirichlet boundary condition can be denoted by the following expression: It is considered that the functions f 1 , f 2 , u 0 , g 1 , and g 2 are smooth enough to simulate the conduction, diffusion, and so on.
u 0 , f 1 , f 2 , g 1 , and g 2 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.

Analytical Solution.
Analytical solution [26,27] of the two dimensional CDR (2) is as follows: 2 International Journal of Differential Equations Within the computations the constrain term h � dx � dy, the initial condition u 0 and the boundary conditions are set from the exact solution. Here, 3) T 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 .

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 (x p , y q ), where x p � a + p.dx for each p � 0, 1, 2, . . . P and y q � c + q.dy for each q � 0, 1, 2, . . . Q, also the lines x � x p and y � y q are grid lines, and their junctions are grid's mesh. For every work point interior of the lattice, (x p , y q ), 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.

2nd
Order Crank-Nicolson Implicit Scheme. By applying the Crank-Nicolson method to (2), by integrating in the compact way, we obtain the following: by substituting the above expressions into (2), the 2nd order implicit method is given by the following: International Journal of Differential Equations e implicit scheme appears that the accuracy is of . Moreover, we use the notation U p,q and u p,q representing the exact and approximate solution of mesh points (x p , y q , t n ) 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. e iterative method is utilized for a linear system. e CN scheme is computationally expensive.

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: In simplified form, the (6) is as follows: Initial Node in Horizontal Direction Final Node in Horizontal Direction where In simplified form the (10) is as follows: Initial Node in Vertical Direction Final Node in Vertical Direction where is scheme has accuracy of 2nd order both in space and time along with the unconditional stability.
e linear system in the horizontal direction is as follows: where m � 1, . . ., M − 1 e linear system in the vertical direction is as follows: where l � 1, . . ., L − 1 where a 1 , c 1 create tridiagonal matrix and the array b 1 , d 1 are dependent upon p and q e reactive term in horizontal-direction So, also for the reactive term in the vertical direction, it is as follows: G(ku * y,l ) � [kG(u * 1,M− 1 , . . . , kG(u * 1,1 )] e notations u * p,q and u n+1 p,q are used in same manner. e scheme gives us 2nd order accuracy in time and space . For stability analysis, this scheme exhibits unconditional stability. e stability of alternating direction implicit strategy (6) and (10) can be obtained by using the von Neumann criteria for stability. e regular stability investigation criteria indicate toward the amplification factor G of (6) and (10) may be observed as by substituting: u n � λ n e i k x rh x +k y sh y .
It is obvious that the amplification factor is always less of equal to |G| ≤ (1). Hence, it has unconditional stability. Eventually, the methodology results in a tridiagonal linear system. Iterative strategies were applied to unravel the abovementioned framework. e trick utilized in developing the ADI scheme is to divide the time step into two, and apply International Journal of Differential Equations 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][29][30][31][32]. To get the approximate solution, we have to find the solution of a nonlinear tridiagonal framework with each passing step of time.

Finite Volume Method.
e finite volume method has a lot of flexibility in solving the transport phenomena equations [33][34][35][36][37][38][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 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.

2nd
Order Finite Volume Scheme. By solving (2) using the finite volume method, we have the following: where Bottom left corner boundary cell Bottom right corner boundary cell Top left corner boundary cell Top right corner boundary cell 6 International Journal of Differential Equations Bottom non-corner boundary cell Left non-corner boundary cell Top non-corner boundary cell Right non-corner boundary cell e scheme gives us an accuracy of second order. For stability criteria, we know that P e � advection strength/ diffusion strength Figure 1. To define the stability of 2nd order finite volume scheme, P e must be less than or equal to 2. e 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.

Algorithm.
Clearly, the framework is tridiagonal and can be fathomed with the omas algorithm. e measurement of H is p × q. In common a tridiagonal framework can be composed as follows: e above system can be expressed as follows in a matrixvector form: where H represents a coefficient matrix. e column vector on the right side is known. Our primary objective is to discover the resultant vector u. Presently, we have the following: · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · · a n− 1 b n− 1 c n− 1 · · · · · · · · · · · · a n b n By comparing both sides of the Hu � S, it gives us the components of the matrices L and U. e computational steps for the execution of the omas calculation are appeared in comes about, taken from particular illustrations.

Error Norms
e precision and consistency of the methods is studied in terms of error norms specifically L 2 and L ∞ which can be expressed as follows: where u (x, y, t) and U (x, y, t) represent the approximate and exact solutions at the grid point (x p , y q , t n ). In this method, ρ(U p,q − u p,q ) � max(λ) and λ is an eigen value of (U p,q − u p,q ).

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. e 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. e convective and diffusive fluxes at the control volume faces are calculated using a finite volume approach that employs consistent formulations. e article's main goal is to identify the accuracy and efficiency of second order numerical methods.
e algorithms developed are compared to (3) analytical findings. e parameters σ � 1, A → � (2, 3) T , 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. e 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. e truncation error is computed using the L 2 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  International Journal of Differential Equations of ϵ. On the other hand, all these approaches have computational problems for smaller values of ϵ ≤ 10 − 3 . Furthermore, it is observed that by reducing step sizes h values in Table 5 with grid size 20 × 20 and ϵ � 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. e 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. e 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,  Error estimation at ϵ � 0.01, k � 0.0001, gridsize � 10 × 10.     Error estimation at grid size � 50 × 50 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.  On 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. For k � 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. e 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. e findings demonstrate that these schemes are reliable for solving PDEs.
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.

Conclusion
e 2D transient scalar CDR problem is solved numerically in this paper. eoretical formulations of the proposed schemes were fully addressed. e three scheme's efficiency,

Exact Solution
Numerical Solution using ADI Scheme International Journal of Differential Equations 13 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. e FVM, on the other hand, ensures that conservation is maintained across the domain. e stated techniques might be used to address similar mathematical challenges in biology, physics, and engineering.
Data Availability e article contains the supporting data for the study's findings.

Conflicts of Interest
e authors declare that they have no conflicts of interest.