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 using and , and the estimated results are compared to the corresponding analytical solution.
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:
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 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:where alongside initial conditions , (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) = 1 (x, y, t), u (x, d, t) = 2 (x, y, t),
where (x, y, t) , in , time interval is . It is considered that the functions , , , , and are smooth enough to simulate the conduction, diffusion, and so on.
, , , , and are the given sufficiently smooth functions and u (x, y, t) may represent heat, diffusion, etc. Here, ; diffusion coefficient is denoted by while reaction coefficient is represented by and denotes the convection coefficient, respectively.
2.1. Analytical Solution
Analytical solution [26, 27] of the two dimensional CDR (2) is as follows:Within the computations the constrain term , the initial condition and the boundary conditions are set from the exact solution. Here, , and . Moreover, to see the impact of diffusion coefficient on the solution, values of diffusion coefficient is selected such that .
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 and to express step sizes, and in respective horizontal and vertical directions. Dividing the interval into equidistant parts of width and the interval into equidistant parts of width . Illustrate a grid by mapping horizontal and vertical axes connecting the coordinate points , where for each and for each , also the lines and are grid lines, and their junctions are grid’s mesh. For every work point interior of the lattice, , for and . We apply different calculations to approximate the numerical arrangement to the issue in this equation (2). We expect , where 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:by substituting the above expressions into (2), the 2nd order implicit method is given by the following:
The implicit scheme appears that the accuracy is of . Moreover, we use the notation and representing the exact and approximate solution of mesh points 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 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 to , and later from to the . 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 In simplified form, the (6) is as follows: Initial Node in Horizontal Direction Final Node in Horizontal Direction where Sweep in Vertical Direction In simplified form the (10) is as follows: Initial Node in Vertical Direction Final Node in Vertical Direction where
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:where m = 1, …, M − 1
The linear system in the vertical direction is as follows:where l = 1, …, L − 1
where create tridiagonal matrix and the array are dependent upon p and q
The reactive term in horizontal-direction
So, also for the reactive term in the vertical direction, it is as follows:
The notations are used in same manner. The scheme gives us 2nd order accuracy in time and space . 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 of (6) and (10) may be observed as by substituting:
It is obvious that the amplification factor is always less of equal to . 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 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.
3.3.1. 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
Bottom non-corner boundary cell
Left non-corner boundary cell
Top non-corner boundary cell
Right non-corner boundary cell
The scheme gives us an accuracy of second order. For stability criteria, we know that Figure 1. To define the stability of 2nd order finite volume scheme, 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.
3.4. Algorithm
Clearly, the framework is tridiagonal and can be fathomed with the Thomas algorithm. The measurement of H is . In common a tridiagonal framework can be composed as follows:The above system can be expressed as follows in a matrix-vector form:where represents a coefficient matrix. The column vector on the right side is known. Our primary objective is to discover the resultant vector . Presently, we have the following:
The method is clarified within the taking after steps,and
By comparing both sides of the , it gives us the components of the matrices and . 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 specifically and 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 . In this method, and is an eigen value of .
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 , , and changing 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 linked to two variables at with a fixed grid size of is found in Table 1. The truncation error is computed using the and norms with a fixed grid size of in Table 2. Table 2 demonstrates that modifying the time with the same grid size at , improves the accuracy of the error norm. Table 3 shows that error decreases as grid size increases at . 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 .
Furthermore, it is observed that by reducing step sizes values in Table 5 with grid size 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 and time level with a fixed grid size of is presented using error analysis norm. Results can be seen using the norm. Variation of the parameter toward smaller values results in an increase in accuracy. When we changed the value of 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.
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 , , and grid size 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 and , 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 and grid size . 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 with smaller time space. Figure 8 shows efficient and reliable numerical results using the ADI technique with grid size , and . 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 along the 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.
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.
Data Availability
The article contains the supporting data for the study’s findings.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Authors’ Contributions
All authors contributed equally.