Numerical Solution of Two-Dimensional Nonlinear Unsteady Advection-Diffusion-Reaction Equations with Variable Coefficients

Tsega, Endalew Getnet

doi:https://doi.org/10.1155/2024/5541066

International Journal of Mathematics and Mathematical Sciences

On this page

Abstract Introduction Discussion Conclusions Data Availability Conflicts of Interest References Copyright Related Articles

Research Article | Open Access

Volume 2024 | Article ID 5541066 | https://doi.org/10.1155/2024/5541066

Numerical Solution of Two-Dimensional Nonlinear Unsteady Advection-Diffusion-Reaction Equations with Variable Coefficients

Endalew Getnet Tsega¹

Academic Editor: Saranya Shekar

Received29 Jan 2024

Revised12 Mar 2024

Accepted21 Mar 2024

Published31 Mar 2024

Abstract

The advection-diffusion-reaction (ADR) equation is a fundamental mathematical model used to describe various processes in many different areas of science and engineering. Due to wide applicability of the ADR equation, finding accurate solution is very important to better understand a physical phenomenon represented by the equation. In this study, a numerical scheme for solving two-dimensional unsteady ADR equations with spatially varying velocity and diffusion coefficients is presented. The equations include nonlinear reaction terms. To discretize the ADR equations, the Crank–Nicolson finite difference method is employed with a uniform grid. The resulting nonlinear system of equations is solved using Newton’s method. At each iteration of Newton’s method, the Gauss–Seidel iterative method with sparse matrix computation is utilized to solve the block tridiagonal system and obtain the error correction vector. The consistency and stability of the numerical scheme are investigated. MATLAB codes are developed to implement this combined numerical approach. The validation of the scheme is verified by solving a two-dimensional advection-diffusion equation without reaction term. Numerical tests are provided to show the good performances of the proposed numerical scheme in simulation of ADR problems. The numerical scheme gives accurate results. The obtained numerical solutions are presented graphically. The result of this study may provide insights to apply numerical methods in solving comprehensive models of physical phenomena that capture the underlying situations.

1. Introduction

The ADR equation is a partial differential equation used to represent various mathematical models in science and engineering. The equation describes transport processes driven by advection, diffusion, and reaction. Advection involves the collective displacement of species propelled by a flow field. Diffusion, characterized by the movement of species driven by concentration gradients, serves to homogenize concentration distributions over time. Reaction process occurring within the system may either generate or deplete the transported species. In practical scenarios, the coupling of advection, diffusion, and reaction processes results in intricate interactions, leading to temporal and spatial alterations within the ADR system. The ADR equation is used in modeling contaminant transport [1, 2], fluid flow [3], mass and heat transfer [4], chemical reaction process [5], diffusion across a cell membrane [6], tumor cell expansion [7], population dynamics [8], and so on.

Because of its importance in many fields, obtaining accurate solution for ADR equation has been the interest of numerous researchers. There have been a large number of solutions for ADR equations presented in literature using analytic [1, 9, 10] and numerical [11–15] techniques. Due to the intricate nature of many advection-diffusion-reaction problems, numerical techniques are commonly employed to obtain solutions for the ADR equation. Various numerical techniques have been employed in the literature to compute numerical solutions for two-dimensional nonlinear ADR equations. For instance, Mesgarani et al. [16] used radial basis functions to solve time-dependent nonlinear ADR equations with variable coefficients. They transformed the ADR equation into system ordinary differential equations and used the fourth-order Runge–Kutta method to compute solution of the system. Similarly, Ali et al. [17] applied a finite difference method incorporating with Lucas and Fibonacci polynomials to obtain the solution of one- and two-dimensional nonlinear ADR equations. Ngondiep [18] employed a six-level time-split Leap-Frog/Crank–Nicolson approach to approximate solutions of two-dimensional nonlinear time-dependent ADR equations. He deduced that the numerical scheme is fast, spatial fourth-order convergent, and temporal second-order accurate. Zhang et al. [19] used fourth-order fully implicit finite difference scheme to solve unsteady nonlinear ADR equations. Tambue [20] presented a finite volume method combined with exponential time differencing to solve the ADR equation involving the nonlinear reaction term.

Numerous researchers have solved ADR equations that model real-realistic problems. Yu et al. [21] utilized the homotopy analysis method (HAM) to simulate a contaminant transport model governed by ADR equations, presenting contaminant concentration profiles at various time points based on simulation outcomes. Para et al. [12] investigated a water pollution scenario represented by ADR equations, employing a finite volume method. Their numerical solution was compared with results obtained using a finite difference method, revealing satisfactory agreement between the solutions. Gurarslan et al. [22] obtained a numerical solution for contaminant transport described by a one-dimensional advection-diffusion equation, employing a finite difference scheme in space and a Runge–Kutta method in time. The integrated scheme yielded accurate solutions. Putri et al. [23] conducted numerical simulations of advection-diffusion equations representing oxygen demand concentration in waste stabilization ponds using a finite difference method. Pananu et al. [24] utilized a finite difference method scheme to investigate a water pollutant dispersion problem.

A Crank–Nicolson method is an implicit finite difference method for solving partial differential equations involving time and space derivatives [25]. In solving partial differential equations using the Crank–Nicolson method, both the time and space derivatives are approximated by central differences. The method is unconditionally stable and second-order accurate both in space and time variables for solving the diffusion equation. However, the method requires more computing time to solve the resulting system of algebraic equations. Specially, when the method is applied to solve nonlinear and multidimensional problems, linearization techniques like Newton’s method [26] are required. Solving the system of nonlinear equations introduced through discretization demands much extra time. Any method solving system of linear equations (matrix-inversion method, Gauss–Seidel method, Lu decomposition method, and so on) can be used to solve the block tridiagonal system to obtain the error correction vector required in Newton’s method.

The objective of this paper is to devise a numerical scheme based on finite difference approximation to solve ADR equations with all of the following properties: (i) unsteady, (ii) two-dimensional, (iii) nonlinear reaction term, and (iv) spatially varying velocity and diffusion coefficients. The intention is to employ the developed numerical scheme to simulate the distribution of concentration of a contaminant within a flowing fluid.

2. Problem Description

In this study, the ADR equation of the following form is considered [27]:with initial conditionand boundary conditionwhere is contaminant concentration or temperature at position (x, y) and time t, is the advection velocity vector, is a tensor expressing the diffusivity of , is the reaction term assumed to be nonlinear function of making the ADR equation nonlinear, and is the boundary of a bounded domain in . We assume contaminant is released in an incompressible fluid so that the continuity equation is satisfied. If we consider and , where u, , , and are velocity components and diffusion coefficients in x- and y-directions, the ADR equation in (1) can be rewritten as

For this work, the spatial domain Ω = [0, 1] × [0, 1] and time domain [0, T] are used.

3. Discretization of ADR Equation

To apply finite difference discretization of (4), we divide the spatial domain into uniform meshes of step sizes and in x- and y-directions, as shown in Figure 1, and the time domain with step size . The grid points are defined as

In equation (4), approximating the time derivative by central difference at point , the derivatives of the diffusion coefficients by central difference, the derivatives of C by the averages of central difference at and time levels, and reaction term using the average values of R at and time levels, we get

Note that in discretization (6), the Crank–Nicolson Method is used. Rearranging equation (6), we obtainwhere , and

For , , equation (7) yields a system of nonlinear equations with unknowns .

4. Newton’s Method

For , define equationswhere the functions include variable terms from time level values and reaction terms and constant terms from the boundary values and time level values. Let be a vector function whose components are the functions in (9). Let be the vector of unknowns in (7). Let be the initial approximation of . Using Newton’s method, an error correction vector such that can be obtained iteratively by solving the transformed linear system:where . The iterative process for solving h in (10) at each time level continues until we obtain at this time level with the required accuracy. The computation is performed until , where Tol is the tolerance of approximation for Newton’s method. The matrix A in (10) is a block tridiagonal matrix. For instance, for sample mesh indicated in Figure 1 with 6 × 6 grid points, A has the following form:

To avoid computation with zero entries in matrix A, the Gauss–Seidel iterative method with sparse matrix computation is used to solve the system in (10). A termination criterion is used in the Gauss–Seidel iterative method, where tol is a prescribed number to bound the error. The values of at k + 2, k + 3, …, time levels are solved similarly.

5. Convergence and Accuracy of the Numerical Scheme

Let L and be the differential operator and the finite difference operator representing equations (4) and (6), respectively. If the diffusion coefficient is constant, the truncation error T.E. of the finite difference scheme in (6) is

The truncation error vanishes as . Hence, the finite difference scheme is consistent with the partial differential equation.

Let us analyse the stability of the finite difference scheme (6) using the Von Neumann analysis method for some cases [25]. Assume velocity and diffusion coefficients are constants and the reaction term is zero. Suppose the discretization of C in (6) at grid point (i, j, k) has the form [13]where , A, , and are a constants, and is the amplification factor. Substituting (12) in (6) and dividing both sides of equations by , we get

Using the facts that and and solving for from (13), we havewhere

For the stability, we must have . From (14), .

Hence, the finite difference scheme (6) is unconditionally stable. According to Lax’s equivalence theorem, the finite difference scheme (6) gives a convergent solution for the considered cases.

In the proposed numerical scheme, the linearized system of equations at each iteration is solved using the Gauss–Seidel Iterative method to obtain error correction vector. Thus, if A is diagonally dominant, the Gauss–Seidel iterative method is convergent. A is diagonally dominant if the following inequality is satisfied for each row in the matrix:where is the initial guess to be taken for iteration and m indicates the component of . In (16), the following indexes have to be used: for all terms in the inequality; for the term in the left side and for the first and second terms in the right side; for the third term; and for the fourth term in the right side.

The accuracy of the numerical scheme is calculated by the absolute maximum error formula given by [28]where and are the approximate and exact solutions, respectively. The rate of convergence for the scheme is calculated using the formula as follows:

Here, and are maximum absolute errors with sets of number of grid points and , respectively, where the set is obtained using half of the spatial and temporal step sizes used in . If the exact solution is not known in (17), the absolute maximum error can be estimated using two approximate solutions asand the rate of convergence can be computed using (19) and (18) accordingly.

6. Validation

To validate the convergent, accuracy, and computational efficiency of the proposed numerical scheme, the advection-diffusion equation in [29] with constant diffusion coefficient and velocity components and is taken. The equation iswith initial conditionand boundary conditions

An approximate solution of (20) is obtained by the numerical scheme with spatial step sizes and time step size and . A tolerance of 0.000001 is employed for both Newton’s method and the Gauss–Seidel iterative method. Figure 2 illustrates the comparison between the approximate and exact solutions. In this computation, the absolute maximum error of the results obtained using the numerical scheme for and is . Furthermore, Figure 3 presents a comparison of the exact and approximate solutions for at x = 0.5 and y = 0.5 and for at y = 0.5 and t = 0.5. Notably, there is a strong agreement between the numerical and exact solutions, as evidenced by the figures. In addition, Table 1 shows the maximum absolute errors, convergence rates, and CPU time. The validation process demonstrates that the proposed numerical scheme yields accurate results and can be utilized for addressing similar partial differential equations. The scheme requires increased computational time when employing very small spatial step sizes.

(a)

(b)

7. Numerical Examples and Discussion

In this section, numerical demonstrations are provided to illustrate the practical utilization of the proposed numerical scheme in solving ADR problems. The computed results are illustrated in graphs, and the estimated maximum absolute errors are presented in tabular form. For all examples, a tolerance of 0.000001 is used for Newton’s method and Gauss–Seidel iterative method in the numerical scheme. All units of quantities in examples are in SI units. Exact solutions are not available for the considered examples.

Example 1. Consider the ADR problem in (4) with , (see velocity streamlines in Figure 4), , [1] (see graphs of Dx and Dy in Figure 5), [30], initial conditionand boundary condition for and Figure 6 displays the solution of Example 1 using the numerical scheme with spatial step sizes and time step size and at and . In solving the error correction vector h at each iteration, the matrix A for this problem is diagonally dominant and hence the numerical scheme yields convergent solution. The maximum absolute error obtained in the computation at t = 2 using (16) is 0.1459. It is observed that as the maximum absolute error decreases, we take smaller step sizes. There is a rapid change of concentration distribution profile.

Example 2. Let us take the ADR problem in (4) with diffusion coefficient and (see Graphs of Dx and Dy in Figure 7) having the same velocity components, reaction term, and initial and boundary conditions as Example 1.
The numerical solution of Example 2 at and is displayed in Figure 8. The same grid size is used as Example 1 for the computation. This example targets to see the effect of diffusion coefficient on the contaminant concentration distribution. As we can observe in Figure 8, higher values of the diffusion coefficient in the contamination transport give wider coverage area of concentration distribution for the same advection velocity which agrees with situation observed in [12]. The maximum value of the concentration at t = 2 is attained near the entrance gate.

Example 3. Consider the ADR problem in (4) with and (see the velocity streamline in Figure 4) and having the same diffusion coefficients, reaction term, and initial and boundary conditions as Example 1.
The numerical solution of Example 3 at and is depicted in Figure 9. The same grid sizes are used as Example 1 to obtain the solution. For this problem, the matrix A that appears at each iteration in solving error correction vector h is also diagonally dominant. Hence, the numerical scheme is convergent. The maximum absolute error obtained in the computation is 0.2336. There is relatively slow change of concentration distribution profile between t = 1 and t = 3 as it can be observed from the graphs. The contaminant concentration distribution profile follows the advection velocity streamlines (see Figures 4 and 9).

Example 4. Let , (see the velocity streamlines in Figure 4), and [29] in (4) having the same diffusion coefficient and initial and boundary conditions as Example 1.
In Example 4, the reaction term is a third degree polynomial. The numerical solution of this problem is presented in Figure 10 at t = 1 and t = 4. The same grid sizes are used as Example 1 in the computation. The numerical scheme yields a convergent solution. The maximum absolute error obtained in the computation is 0.1615.
For all examples, the maximum absolute errors are computed for the four sets of number of grid points , , , and for and using (16) as indicated in Table 1. For all problems, the maximum absolute error decreases as grid sizes decrease which shows the accuracy and the convergent of the proposed numerical scheme in solving the problems. Figure 4 shows the velocity stream lines of the three variable velocities in the examples. These streamlines are plotted using the streamline function of MATLAB. From the numerical results of problems in Examples 1, 3, and 4, it is observed that the contaminant concentration distribution profiles follow the fluid flow velocity streamlines which agreed with the condition realized in [1] for small diffusion coefficient. This reveals the applicability of the proposed numerical scheme in simulations of ADR problems.

8. Conclusions

In this study, a numerical scheme is employed to obtain the solution of two-dimensional unsteady nonlinear advection-diffusion-reaction equations with velocity and diffusion coefficients that vary spatially. By developing and implementing the comprehensive numerical scheme, which integrates the Crank–Nicolson Method, Newton’s method, and the Gauss–Seidel iterative method, nonlinear ADR problems are successfully solved. The absolute maximum errors and convergence rate of the numerical scheme are estimated to investigate the accuracy and efficiency of the numerical scheme. The findings demonstrate the scheme’s capability to accurately approximate solutions for ADR problems and its suitability for simulating related two-dimensional nonlinear partial differential equations. The scheme requires substantial CPU time during computation when a large number of spatial grid points are utilized.

Data Availability

No data were used to support the findings of this study.

Conflicts of Interest

The author declares that there are no conflicts of interest.

References

C. Zoppou and J. H. Knight, “Analytical solution of a spatially variable coefficient advection–diffusion equation in up to three dimensions,” Applied Mathematical Modelling, vol. 23, no. 9, pp. 667–685, 1999.
View at: Publisher Site | Google Scholar
L. Li, Z. Jiang, and Z. Yin, “Fourth-order compact finite difference method for solving two-dimensional convection–diffusion equation,” Advances in Difference Equations, vol. 2018, no. 1, p. 234, 2018.
View at: Publisher Site | Google Scholar
D. A. Garzón-Alvarado, C. H. Galeano, and J. M. Mantilla, “Computational examples of reaction–convection–diffusion equations solution under the influence of fluid flow: first example,” Applied Mathematical Modelling, vol. 36, no. 10, pp. 5029–5045, 2012.
View at: Publisher Site | Google Scholar
M. Arshad, A. Hussain, A. Hassan et al., “Heat and mass transfer analysis above an unsteady infinite porous surface with chemical reaction,” Case Studies in Thermal Engineering, vol. 36, Article ID 102140, 2022.
View at: Publisher Site | Google Scholar
M. R. Eid, “Chemical reaction effect on MHD boundary-layer flow of two-phase nanofluid model over an exponentially stretching sheet with a heat generation,” Journal of Molecular Liquids, vol. 220, pp. 718–725, 2016.
View at: Publisher Site | Google Scholar
W. Stein, Transport and Diffusion across Cell Membranes, Elsevier, Amsterdam, Netherlands, 2012.
N. Meghdadi, M. Soltani, H. Niroomand-Oscuii, and N. Yamani, “Personalized image-based tumor growth prediction in a convection–diffusion–reaction model,” Acta Neurologica Belgica, vol. 120, no. 1, pp. 49–57, 2020.
View at: Publisher Site | Google Scholar
S. Komala Sheshachala and R. Codina, “Finite element modeling of nonlinear reaction–diffusion–advection systems of equations,” International Journal of Numerical Methods for Heat and Fluid Flow, vol. 28, no. 11, pp. 2688–2715, 2018.
View at: Publisher Site | Google Scholar
K. Zamani and F. A. Bombardelli, “Analytical solutions of nonlinear and variable-parameter transport equations for verification of numerical solvers,” Environmental Fluid Mechanics, vol. 14, no. 4, pp. 711–742, 2014.
View at: Publisher Site | Google Scholar
J. Pérez Guerrero, E. M. Pontedeiro, M. T. van Genuchten, and T. H. Skaggs, “Analytical solutions of the one-dimensional advection–dispersion solute transport equation subject to time-dependent boundary conditions,” Chemical Engineering Journal, vol. 221, pp. 487–491, 2013.
View at: Publisher Site | Google Scholar
E. Ngondiep, “A novel three-level time-split approach for solving two-dimensional nonlinear unsteady convection-diffusion-reaction equation,” The Journal of Mathematics and Computer Science, vol. 26, no. 3, pp. 222–248, 2021.
View at: Publisher Site | Google Scholar
K. Para, S. Sungnul, S. Sirisubtawee, and S. Phongthanapanich, “A comparison of numerical solutions for advection-diffusion-reaction equations between finite volume and finite difference methods,” Engineering Letters, vol. 30, no. 2, 2022.
View at: Google Scholar
J. Wei, Y. Ge, and Y. Wang, “High-order compact difference method for solving two and three-dimensional unsteady convection diffusion reaction equations,” Axioms, vol. 11, no. 3, p. 111, 2022.
View at: Publisher Site | Google Scholar
A. Singh, S. Das, S. H. Ong, and H. Jafari, “Numerical solution of nonlinear reaction–advection–diffusion equation,” Journal of Computational and Nonlinear Dynamics, vol. 14, no. 4, Article ID 041003, 2019.
View at: Publisher Site | Google Scholar
U. Ahmed, D. S. Mashat, and D. A. Maturi, “Finite volume method for a time-dependent convection-diffusion-reaction equation with small parameters,” International Journal of Differential Equations, vol. 2022, Article ID 3476309, 15 pages, 2022.
View at: Publisher Site | Google Scholar
H. Mesgarani, M. Kermani, and M. Abbaszadeh, “Application of SPD-RBF method of lines for solving nonlinear advection–diffusion–reaction equation with variable coefficients,” International Journal of Numerical Methods for Heat and Fluid Flow, vol. 32, no. 2, pp. 850–886, 2022.
View at: Publisher Site | Google Scholar
I. Ali, S. Haq, K. S. Nisar, and S. U. Arifeen, “Numerical study of 1D and 2D advection-diffusion-reaction equations using Lucas and Fibonacci polynomials,” Arabian Journal of Mathematics, vol. 10, no. 3, pp. 513–526, 2021.
View at: Publisher Site | Google Scholar
E. Ngondiep, “A six-level time-split leap-frog/crank–nicolson approach for two-dimensional nonlinear time-dependent convection diffusion reaction equation,” International Journal of Computational Methods, vol. 20, no. 8, Article ID 2250064, 2023.
View at: Publisher Site | Google Scholar
L. Zhang, Y. Ge, X. Yang, and Y. Zhang, “High accuracy compact difference and multigrid methods for two-dimensional time-dependent nonlinear advection-diffusion-reaction problems,” International Journal of Computer Mathematics, vol. 100, no. 7, pp. 1552–1579, 2023.
View at: Publisher Site | Google Scholar
A. Tambue, “An exponential integrator for finite volume discretization of a reaction–advection–diffusion equation,” Computers and Mathematics with Applications, vol. 71, no. 9, pp. 1875–1897, 2016.
View at: Publisher Site | Google Scholar
C. Yu, A. Deng, J. Ma, X. Cai, and C. Wen, “Semi-analytical solutions for two-dimensional convection–diffusion–reactive equations based on homotopy analysis method,” Environmental Science and Pollution Research, vol. 25, no. 34, pp. 34720–34729, 2018.
View at: Publisher Site | Google Scholar
G. Gurarslan, H. Karahan, D. Alkaya, M. Sari, and M. Yasar, “Numerical solution of advection-diffusion equation using a sixth order compact finite difference method,” Mathematical Problems in Engineering, vol. 2013, Article ID 672936, 7 pages, 2013.
View at: Publisher Site | Google Scholar
G. A. Putri, S. Sunarsih, and S. Hariyanto, “Numerical simulation of advection-diffusion on flow in waste stabilization ponds (1-dimension) with finite difference method forward time central space scheme,” Environmental Engineering Research, vol. 23, no. 4, pp. 442–448, 2018.
View at: Publisher Site | Google Scholar
K. Pananu, S. Sungnul, S. Sirisubtawee, and S. Phongthanapanich, “Convergence and applications of the implicit finite difference method for advection-diffusion-reaction equations,” IAENG International Journal of Computer Science, vol. 47, no. 4, pp. 645–663, 2020.
View at: Google Scholar
R. S. Gupta, Elements of Numerical Analysis, Cambridge University Press, Cambridge, UK, 2nd edition, 2015.
S. S. Ray, Numerical Analysis with Algorithms and Programming, Chapman and Hall/CRC, Nework, Boca Raton, FL, USA, 1st edition, 2016.
B. McLaughlin, J. Peterson, and M. Ye, “Stabilized reduced order models for the advection–diffusion–reaction equation using operator splitting,” Computers and Mathematics with Applications, vol. 71, no. 11, pp. 2407–2420, 2016.
View at: Publisher Site | Google Scholar
P. Roul and V. Rohil, “A high-order numerical scheme based on graded mesh and its analysis for the two-dimensional time-fractional convection-diffusion equation,” Computers and Mathematics with Applications, vol. 126, pp. 1–13, 2022.
View at: Publisher Site | Google Scholar
S. A. Al-Bayati, “Boundary element analysis for convection-diffusion-reaction problems combining dual reciprocity and radial integration methods,” Brunel University, London, UK, 2018, Ph.D. thesis.
View at: Google Scholar
V. S. Aswin and A. Awasthi, “A robust numerical scheme for the simulation of nonlinear convection–diffusion–reaction equation,” International Journal for Computational Methods in Engineering Science and Mechanics, vol. 20, no. 5, pp. 347–357, 2019.
View at: Publisher Site | Google Scholar

Copyright

Copyright © 2024 Endalew Getnet Tsega. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

PDF Download Citation

Download other formats

Order printed copies

Views

216

Downloads

109

Citations