Finite element-based hybrid techniques for advection-diffusion-reaction processes

Article history: Received: 31 January 2017 Accepted: 24 January 2018 Available Online: 4 February 2018 In this paper, numerical solutions of the advectio n-diffusion-reaction (ADR) equation are investigated using the Galerkin, collo ati n and TaylorGalerkin cubic B-spline finite element methods in strong form of spatial elements using an α-family optimization approach for time variation. T he main objective of this article is to capture effective results of the fini te element techniques with Bspline basis functions under the consideration of t he ADR processes. All produced results are compared with the exact soluti on and the literature for various versions of problems including pure advecti on, pure diffusion, advectiondiffusion, and advection-diffusionreaction equations. It is proved that the present methods have good agreement with the exact solution and the literature.


Introduction
Consider the following advection-diffusion-reaction equation with given initial and boundary conditions: (4) Many quantities are encountered in various field of science such as mass, heat, energy, velocity, and concentration represented in the advection-diffusionreaction (ADR) equation as the dependent variable ‫.ܥ‬The ADR equation has great importance in different areas, especially those involving fluid flow [1,2].The ADR equation models various physical and chemical processes, as stated in the literature [3], such as heat transfer in draining film, dispersion of tracers in porous media, the spread of pollutants in rivers and stream, the dispersion of dissolved material in estuaries and coastal sea, etc.When the advection is dominant to the diffusion in the equation, the exact solutions mostly fail and thus diverge.In these cases, the effective numerical methods need to be constructed to obtain accurate and stable results of the model equation.Nowadays, B-spline basis functions are main interest of many researchers to find out effective numerical solutions of partial differential equations [4,5].
Various versions of finite element methods have profoundly been analyzed in the literature.For instance; Least-squares B-spline finite element method was used by Dag et al. [6], a cubic B-spline collocation method was introduced by Goh [7], an upwind finite element method was organized by Ramakrishman [8], the quartic and quintic B-spline methods were used by Korkmaz and Dag [9] for their own problems.In the study of Irk et al. [3], an extended cubic B-spline collocation method was also considered.In addition to finite element-based methods, some other numerical methods were also taken into consideration in dealing with the ADR processes [1,10].This study discovers some finite element based hybrid techniques to analyze the model problems encountered in broad range of science.To integrate the resulted system of ordinary differential equations α-family of time approximation is performed and fully discrete algebraic equations are obtained in terms of the parameters.Note that the strong form of the ADR equation ( 1) is accepted, as opposed to the weak form commonly used in the literature, since the strong form leads to computationally more economic and more accurate results.All produced results are compared with the literature and exact solutions.Various test problems involving pure advection, pure diffusion, advection-diffusion and advection-diffusion-reaction are demonstrated with quantitatively and qualitatively produced results.

Galerkin method
To solve equation (1) with given boundary conditions (3)-( 4) and initial condition (2), the Galerkin cubic Bspline finite element method is used for spatial approximation.The selection of these types of basis functions is very suitable and advantageous.The wellknown advantages of using cubic B-splines are the continuity of the approximate solution and the first and second order-derivatives at all region.The interval [ܽ, ܾ] is partitioned into ܰ finite elements.Each element has equal length ℎ and element nodes are discretized as ܽ = ‫ݔ‬ < ‫ݔ‬ ଵ < ⋯ < ‫ݔ‬ ே = ܾ where ‫ݔ‬ ାଵ = ‫ݔ‬ + ℎ (݈ = 0,1.… , ܰ − 1).Let ߮ be the cubic B-spline basis functions [11] as follows The corresponding cubic B-spline basis functions include the set of splines {߮ ିଵ , ߮ , … , ߮ ேାଵ } and the global approximation function ‫ܥ‬ ሚ ே ‫,ݔ(‬ ‫)ݐ‬ can be expressed as ) where ߚ ሺ‫ݐ‬ሻ is the time part of approximation function ܿ̃ேሺ‫,ݔ‬ ‫ݐ‬ሻ and is to be determined from the time approximation.To compute element matrices, it is required to use local coordinate system considering (5) and ߪ = ‫ݔ‬ − ‫ݔ‬ where 0 ≤ ߪ ≤ ℎ , the basis functions are expressed in the following form Each finite element [x ୪ , x ୪ାଵ ] is covered by the set of four cubic B-splines {߮ ିଵ , ߮ , ߮ ାଵ , ߮ ାଶ }.
Now it is time to apply the Galerkin approach.By considering element ‫ݔ[‬ , ‫ݔ‬ ାଵ ], let us consider the strong form of equation ( 1) over the interval ‫ݔ[‬ , ‫ݔ‬ ାଵ ], one can then write The test function ‫ݒ‬ is selected to be equal to the cubic B-spline basis functions.This approach is known as the Galerkin approach in the finite element society.Use of (8) and local coordinate system (7) transforms equation (10) to the following relation or in a matrix notation where In equation (11); ‫ܯ‬ , ‫ܭ‬ and ‫ܮ‬ are (4 × 4) time independent matrices.After the assembling process of each element and imposing the boundary conditions the matrix form will finally be where ‫ܯ‬ * , ‫ܮ‬ * and ‫ܭ‬ * are (ܰ + 1) × (ܰ + 1) matrices, ߚ = (ߚ , … , ߚ ே ) ் is the unknown time approximation vector and R is an ((ܰ + 1) × 1) time dependent vector.By considering Dirichlet boundary conditions ሺ3ሻ − ሺ4ሻ, ܴ is defined as follows Note that ‫ܯ‬ * * , ‫ܮ‬ * * and ‫ܭ‬ * * are the assembled matrices before imposing the boundary conditions.Thus, equation ( 14) is a system of ordinary differential equations, which is integrated using α-family of time approximation.

Taylor-Galerkin method
The third approximation method in solving equation ( 1) is the Taylor-Galerkin method being effective for many problems represented by differential equations.
The main idea of the method is that the time approximation based on Taylor series expansion is performed before the spatial discretization.After performing the time discretization, the Galerkin method is used for the spatial approximation by utilizing B-splines basis functions (5).The order of the TGFEM schemes can be determined by the truncation error of the Taylor expansion.In this study, we prefer to use the second order TGFEM schemes for the numerical solution of equation (1).Use of the Taylor expansion of the function ‫ܥ‬ with respect to ‫ݐ‬ gives rise to, Taking derivative of equation ( 1) with respect to ‫ݐ‬ leads to, .This selection has the same as the Galerkin method.For the sake of brevity, it is preferred to use a different way to approximate the second time derivative of function ‫.ܥ‬It is noticeable that, the first order time derivative of the function, i.e. ‫ܥ‬ ௧ , can be rewritten by using ADR equation (1) itself.Euler time stepping is used for the diffusion term of equation (25) while the rest of the terms are dealt with the considered technique.Thus, equation ( 25) can be re-expressed as follows, Substitution of equations ( 1) and (26) into equation ( 24) and doing some mathematical manipulations lead to the following semi-discrete system, Using the Galerkin approach for equation ( 27) under the consideration of equations ( 5)-( 8) and doing some algebraic operations, one obtains the following iteration system, and After assembling the procedure and imposing the Dirichlet boundary conditions, the matrix equation will then be Note that ‫ܯ‬ * * and ‫ܭ‬ * * are the assembled matrices before prescribing the boundary conditions.Equation (30) is a recursive relation between ߚ and ߚ ାଵ .By obtaining ߚ we can calculate the solution vector for each time step.

α-family of time approximation
To solve the ODE systems ( 14) and ( 21), α-family of time approximation is preferred since the method is easy to implement, satisfies the unconditional stability by the dependence of the selection of the parameter α and has the required accuracy.As stated in [12], the αfamily of approximation can be defined as where 0 ≤ ߙ ≤ 1, ‫ݐ‬ ௦ାଵ − ‫ݐ‬ ௦ = ‫ݐ݀‬ and ߚ ሶ stands for the time differentiation.Using the same steps of the procedure given in [12], both equation ( 14) and equation ( 21) give where ‫ݏ‬ represents the time index.{ߚ} can be obtained under the consideration of initial condition (2).Then, by using recursive relation (34), the other solutions are computed.

Numerical illustrations
This section is devoted to numerical illustration of the various test problems for the advection-diffusionreaction processes by considering quantitative and qualitative results.Accuracy and stability of the obtained results are figured out by demonstrating error norms and pointwise solutions.Produced results are compared with the literature and exact solutions.To evaluate error norms of the present results we prefer to use the following norm definitions, .
where ߩ is the real problem parameter.This solution construct a transportation of an initial concentration of 10 height units whose peak value is at ‫ݔ‬ initially along an infinitely long channel as well as it maintains its own shape during the propagation.The parameters of the problem are taken to be ܸ ൌ 0.5, ‫ݔ‬ ൌ 2000 ݉ and ߩ ൌ 264.To compare with the literature [6], all parameters are taken to be equal.The final propagation time is ‫ݐ‬ ൌ 9600 ‫ݏ‬ while the initial and boundary conditions are as follows ‫ܥ‬ሺ0, ‫ݐ‬ሻ ൌ 0, In Figure 1, we demonstrate the propagation of the initial pulses up to 9600 s by considering the Galerkin method for the parameters ݄ ൌ 50 and ‫ݐ݀‬ ൌ 10.The comparison of the absolute errors produced by the Galerkin, collocation and Taylor-Galerkin methods are given in Figure 2 for the values of ܸ ൌ 0.5, ݄ ൌ 100 ݉, ‫ݐ݀‬ ൌ 50 ‫ݏ‬ and ‫ݐ‬ ൌ 9600 ‫.ݏ‬As seen in Figure 2, the Galerkin method has been seen to be more accurate than the rest of the considered methods.
The error norms and peak location of the concentration are compared with the literature [6] and exact peak location in Table 2 for various values of the Courant numbers, i.e. ‫ݎܥ‬ ൌ ሺܸ݀‫ݐ‬ሻ/݄.Because of the stability condition in Table 2, the Taylor-Galerkin method is not preferred to use.The method is not stable for higher Courant numbers.with the following exact solution [13], ‫,ݔ‪ሺ‬ܥ‬ ‫ݐ‬ሻ ൌ expሺെ‫ݐ‬ሻ sin ሺߨ‫ݔ‬ሻ.
(38) The problem has homogenous Dirichlet boundary conditions and initial condition can be taken from the exact solution (38).In Table 3, we compare maximum error norms of the present schemes with the literature [13] and among each other.As seen in Section 2, the same discretized equations have been obtained for the Galerkin and the Taylor-Galerkin methods in case of pure diffusion.As realized in Table 3, the present Galerkin scheme produces better accuracy comparison to the literature [13] and the current collocation scheme.Thus the computed results have been seen to represent the related physical problem.Yet, comparison of absolute errors has been seen both qualitatively and quantitatively in Figure 3 for ݄ ൌ 0.01 ݉, ‫ݐ݀‬ ൌ 0.0001 ‫ݏ‬ and ‫ݐ‬ ൌ 1 ‫.ݏ‬ Figure 4 illustrates the diffusion process of the initial concentration with the diffusion constant ‫ܦ‬ ൌ 1/ߨ ଶ by using the Galerkin approach for ݄ ൌ 0.01 ݉ and ‫ݐ݀‬ ൌ 0.0001 ‫.ݏ‬ The boundary conditions of the problem can be written from exact solution (40) as follows: 2 ( 0.5) ( ) exp 0.00125 0.025 (0.5 ) ( ) exp 0.00125 0.04 0.000625 0.02 0.025 (1.5 ) ( ) exp 0.00125 0.04 0.000625 0.02 The produced results and exact solution are compared in Table 4 for the various values of the Courant number and various values of α at the peak location of the concentration, ‫ݔ‬ ൌ 0.5 and ‫ݐ‬ ൌ 1.The present solutions are compared with the work of Kadalbajoo and Arora [14] in Table 5 under the consideration of the values ൌ 0.1, ݄ ൌ 0.01, ܲ݁ ൌ 1 and ‫ݐ‬ ൌ 1.As seen in Tables 4-5, the Galerkin method seems to be more accurate than the other suggested methods.It is also seen that the Galerkin and the Taylor-Galerkin methods are preferable comparison to the literature [14].For various choices of the parameters, qualitative behavior of the model problem is sketched in Figures 5-6.There has been seen to be good agreement among the suggested methods (see Figure 6).Thus the produced results showed that the current problem has represented the physical behaviour very well.It is important to note that the collocation method takes less computational time than both the Galerkin and the Taylor-Galerkin methods while the collocation method is of a bit cruder results than the others.For the sake of simplicity, the upper bound of spatial domain is taken to be ‫ܮ‬ ൌ 1.For various values of the parameters, absolute errors of the considered numerical methods are compared in Figures 7-8.

Finite
element-based hybrid techniques for advection-diffusion-reaction processes 131In the ADR equation, there are two important problem parameters need to be considered, the Peclet and the Courant numbers.Non-dimensional parameter, Courant ‫)ݎܥ(‬ number, gives the fractional distance relative to the grid spacing travelled due to advection in a single time step ‫ݎܥ‬ ൌ ሺܸ݀‫ݐ‬ሻ/݄.This parameter especially plays an important role when we need to determine stability conditions of the considered numerical approaches.The Peclet number is another crucial non-dimensional parameter which compares the characteristic time for dispersion and diffusion given a length scale with the characteristic time for advection, i.e. ܲ݁ ൌ ሺܸ݄ሻ/‫ܦ‬ where the parameters are as in equation (1).Problem 1[6] Pure advection in an infinitely long channel: Initially, we consider pure advection problem, i.e. ‫ܦ‬ ൌ 0 and ߠ ൌ 0. The analytic solution of the problem of interest is as follows[6]

Figure 1 .
Figure 1.Propagation of initial pulse with constant wave speed ܸ ൌ 0.5 for various time values up to 9600 s.

Problem 4 .
Advection-diffusion-reaction process on a finite lineLet us now consider the ADR equation with arbitrary values of ܸ, ߠ and with homogeneous Dirichlet boundary conditions.Taking following initial condition ‫,ݔ‪ሺ‬ܥ‬ 0ሻ ൌ sin ሺߨ‫ݔ‬ሻ (42) leads us to derive the following Fourier series exact solution,

Table 1 .
Values of approximate function and its derivatives at the end points of the element

Table 3 .
Comparison of the error norms produced with various values of the Peclet number and ‫ݎܥ‬ ൌ 0 in Problem 2.

Table 2 .
Comparison of the error norms produced with various values of the Courant number and ܲ݁ ൌ ∞ in Problem 1.

Table 5 .
[14]mparison of present solutions with the literature[14]and exact solution at various nodal points and at ‫ݐ‬ ൌ 1 in Problem 3.