A Comparative Study on Stabilized Finite Element Methods for the Convection-Diffusion-Reaction Problems

,


Introduction
The modeling of real world problems to predict physical quantities is an important subject in science and engineering.Many physical, chemical, and biological phenomena contain interactions of convection, diffusion, and reaction processes, which can be modeled in terms of the convection-diffusionreaction (CDR) problems [1][2][3].Explicit solutions to those problems are only known under specific circumstances; therefore numerical simulation of the problem is required.However, the numerical approximation of those equations is a challenging problem, in general, and it is not possible to get reasonable approximations with the standard numerical algorithms.Therefore developing efficient and effective computational methods for solving those problems has been drawing attention of many researchers for several decades [4][5][6][7].
The exact solutions of the convection-diffusion-reaction (CDR) problems exhibit very fine structures, so-called layers that are small subregions where the derivatives of the solution are very large.In this case, it is well-known that the Galerkin finite element method produces nonphysical oscillations polluting all over the domain.To improve the numerical approximations, a new class of the finite element methods, the so-called stabilized FEMs, have been introduced.The main idea underlying those methods is to augment the variational formulation by mesh-dependent terms in order to gain control over the derivatives of the solution.The Streamline-Upwind Petrov-Galerkin (SUPG) method proposed by Brooks and Hughes in [8] is one of the most well-known stabilized formulations in the literature.It is extensively used in the simulation of many real world problems.There are also several variants of the SUPG method in the finite element literature which are also popular [9,10].Not only can those methods be applied to many physical problems, but also their computer implementation is quite practical.However the lack of formalism in derivation of the stability parameters used in stabilized formulations has been seen as a major drawback of those strategies.
Another important and a more recent class of the stable discretization for the CDR problem is based on augmenting the finite element space by a set of suitable enriching functions.As an important and well-known example to that class, the residual-free bubbles (RFBs) method could be mentioned that is first proposed by Brezzi and Russo in [11].The RFB method is based on enlarging the polynomial finite element space by special type of enriching functions, the so-called residual-free bubbles, and then eliminating those functions from the formulation to get an improved approximation for the polynomial part of the solution.Although the RFBs form an effective platform to derive improved discretizations, they are defined by a set of local differential equations posed inside each element which may not be easier to solve than the original one, except that the problem domains are simple element geometries [12].That observation has motivated the introduction of a further option, the so-called pseudo residual-free bubble (PRFB) method.In that approach, the RFB functions are basically approximated by piecewise linear functions on an appropriately chosen simple subgrid established inside each element and then they are replaced by their approximate counterparts, the so-called pseudo residual-free bubbles, in the numerical computations.Here, the locations of the subgrid nodes are of critical importance and therefore they should be chosen specially so that the fine scale-effect of the exact solution can accurately be represented in the coarse scale numerical approximation [13][14][15].In particular, their location is determined by minimizing the residual of a local differential problem defining the bubbles, with respect to the  1 -norm.It seems that the PRFB method is quite robust and effective algorithm in the numerical approximation of the convection-diffusion-reaction equations.
In this work, we compare the performances of the Streamline-Upwind Petrov-Galerkin (SUPG), Galerkin/ Least-Squares (GLS), Subgrid Scale (SGS), Link-Cutting Bubble (LCB) [16], and pseudo residual-free bubble (PRFB) methods for the numerical solution of the convectiondiffusion-reaction equation on several benchmark problems.A wide range of problem parameters have been examined and the corresponding numerical results are presented.As it can be seen from the numerical results, the results are comparable especially in the convection-dominated regimes.However the approximations generated by the PRFB method are slightly better than the others in the reaction-dominated cases.
The outline of this paper is as follows.In Section 2, we state the model problem, the variational form, and its standard Galerkin discretization.This is followed by a brief survey to the basic ideas of the stabilizations through augmented variational forms in Section 3. In Section 4, we describe the idea of augmenting finite element space which is the starting point for deriving the pseudo residual-free bubbles.We display the computation of the PRFBs and the derivation of the subgrid points in Section 5. Finally, we perform the numerical tests for several benchmark problems in both 1D and 2D in Section 6.

Statement of the Problem
The convection-diffusion-reaction problem in two space dimensions consists of finding the scalar-valued function  in a polygonal domain Ω with a smooth boundary Γ = Ω and Dirichlet boundary conditions, such that where  > 0 is the diffusion coefficient,  is the velocity field,  > 0 is the reaction coefficient, and  is a smooth source function in  2 (Ω).This model may describe how the concentration of the substance distributed in the medium changes under the influence of interactions of convection, diffusion, and reaction processes.The variational formulation corresponding to problem (1) reads as follows: Find  ∈  1 0 (Ω) such that where is a continuous and coercive bilinear form on the Hilbert space  1 0 (Ω).Let us consider a partition T ℎ of the computational domain Ω into triangular elements  that do not overlap and let ℎ  = diam() with ℎ = max ∈T ℎ ℎ  .If we restrict ourselves to the case of continuous, piecewise linear elements which we denote by then the standard Galerkin finite element method reads as follows: Find   ∈   such that It is well-known that if the problem is convection or reactiondominated, then, unless ℎ is of the same or smaller size of , the solution of (5) possesses spurious oscillations spreading all over the domain.Therefore a stabilization is needed.

Stabilizations through Augmented Variational Forms
A major class of the stabilized finite element methods is based on enlarging the discrete variational formulation by a mesh-dependent stabilization term.The general form of those methods for problem (5) reads as follows: Find   ∈   such that where (⋅, ⋅) indicates a consistent stabilizing term added to the standard Galerkin formulation (5).The general form of the operator (⋅, ⋅) is given by where   is the stabilization parameter that needs to be tuned up by user and (⋅, ⋅)  denotes the restriction of the inner product over .Each choice of the operator L gives rise to The SUPG method was designed by Brooks and Hughes [8] and was later generalized for multidimensional problems by Hughes and Mallet [17].This method adds up additional diffusion to the Galerkin problem (5) in the direction of streamlines.A generalization of the SUPG method is the GLS method whose essential feature is to augment the Galerkin problem (5) with a least-squares form of the residuals that are based on the corresponding differential equations.The SGS method was introduced in [18] which is similar to the GLS method yet is derived from a more theoretical setting.There are also other stabilized methods used in the literature and a thorough comparison of those methods can be found in [19][20][21][22][23].
One of the important points in the implementation of the stabilized methods is the choice of the stabilization parameter   .The parameter   is constructed by user based on the problem and/or algorithm properties (the discrete maximum principle, convergence and stability analysis, etc.).Some of the possibilities that have been proposed for   in the literature [21,24,25] are as follows: where An appropriate choice from the stability parameters ( 9) can be used for the selected method through (8) in computations.In the numerical calculations herein, we will employ the classical stabilized finite element formulations with the parameter    in ( 9) for all test problems under investigation, as the others produce comparable results.
The fact that the stabilized methods can be applied to many physical problems and their computer implementation is quite practical made those methods favorable to many researchers.However the lack of general method in derivation of the stability parameters has been seen as a major drawback of those strategies.Therefore alternative algorithms have still been of interest to researchers.

Stabilization through Augmented Spaces
Another successful class of the stabilized finite element methods is based on enlarging the polynomial finite element space   with more rich functions inspired by the original differential problem.Employing special type of enriching functions, the so-called bubble functions, it is possible to get improved approximations without increasing the size of the discrete problem.The method basically consists of defining the enriched finite element space  ℎ =   ⊕   where the enriching space   is defined by and solving problem (5) on  ℎ .In this approach, it is required that the bubble component   of  ℎ satisfies the original differential equation in each  strongly; that is, With that choice, it is possible to eliminate the bubbles from the formulation, keeping the size of the computation matrices the same, yet getting an improved approximation.This approach is known as the residual-free bubbles (RFBs) method in the finite element literature [11,12], which, for the present problem, reads as follows: The term (  , V  ) is responsible for the stabilization of the numerical method and the bubble component   should be computed before we solve (13) for its linear part   .However, since   is identified by the linear part   and the source function  through (12), its solution is not an easy task, in general.
On the other hand, regarding the configuration of the problem and the simple element geometry, it is possible to compute a reasonable approximation to   that provides a similar stabilizing effect into the numerical method (13).Indeed, in a recent approach based on the RFB method, the bubble functions are replaced by their suitable approximate counterparts, the so-called pseudo bubbles [14,15].The development and the details of that approach are given in the following section.

The Pseudo Residual-Free Bubbles and Stabilizing Subgrids
In this section, we shall display the computation of the pseudo residual-free bubbles and the construction of the associated subgrid nodes on which the bubble functions are approximated (see [15] for details).We note that the configuration of the subgrid nodes is of crucial importance and therefore should be determined in a special manner.Before we outline the algorithm, we introduce the notations used throughout the presentation.Denote the vertices of  by   , the edges of  by   opposite to   , the length of   by |  |, the midpoint of edge   by   , the outward unit normal to   by   , ]  = |  |  ,  ]  = (, ]  ), and the area of  by || (see Figure 1).The actual numbering of the vertices will be chosen according to the direction of .Now let us assume that the location of   is along the median from   ; that is, To specify the problem regimes, we further introduce Consider the bubble functions   ,  = 1, 2, 3, defined by where   are the restrictions of the piecewise linear basis functions for   to  and Now if where   ∈ R,  = 1, 2, 3. Thus, ( 12) is automatically satisfied with the present choice of bubble functions.Using the element geometry and the problem properties, it is possible to construct the approximate bubbles, say  *  , that have the same qualitative behavior with their continuous counterpart.To start with, let   be a piecewise linear function with   (  ) = 0 and   (  ) = 1 ∀,  = 1, 2, 3 and  *  () =     () be the classical Galerkin approximation of   through (16); that is, Here the choice of the internal point   is critical and the main criteria used to determine its location are to minimize  1 norm of the residual coming out from the bubble equation ( 16) in the critical case where a layer structure exists.In other words, we choose   such that is minimum.The criteria (20) are employed to determine the locations of subgrid nodes   in convection and reactiondominated cases separately.

Convection-Dominated Flows.
In order to present the explicit descriptions of   , we have to distinguish among the following two cases, in which the convection-dominated regime is considered with respect to the number of inflow edges (see Figure 2).

One Inflow Edge.
The problem is assumed to be convection-dominated if From a previous calculation [15], the explicit location of  1 along the median from  1 is given by where . The choice of other points  2 and  3 should be consistent with the physics of the problem.Thus we take

Two Inflow Edges.
The problem is assumed to be convection-dominated if From [15], the position of  2 along the median from  2 is given by where , and the position of  3 along the median from  3 is given by where The choice of  1 should be consistent with the physics of the problem.Thus we take With the present choice, it is noteworthy to observe the self-adjustment of subgrid nodes as the problem evolves from the diffusion-dominated regime into the convectiondominated one, so that the pseudo bubbles contribute to the stability of the numerical method at their most (see Figure 3).

Reaction-Dominated Flows.
In order to present the explicit descriptions of   in the reaction-dominated regime, we have to distinguish among the following two cases with respect to the number of inflow edges (see Figure 4).

One Inflow Edge. The problem is assumed to be reaction-dominated if
The location of  1 is taken from (22) and the position of  2 along the median from  2 is given by where , and the position of  3 along the median from  3 is given by where [15]).Thus the positions of   ,  = 1, 2, 3, are determined by taking

Two Inflow Edges.
In this case, the problem is assumed to be reaction-dominated if and the locations of  2 and  3 are taken from ( 24)- (25).Now, the position of  1 along the median from  1 is given by where . Thus the positions of   ,  = 1, 2, 3, are determined by taking With the present choice, it is noteworthy to observe the self-adjustment of subgrid nodes as the problem evolves from the convection-dominated regime into the reactiondominated one, so that the pseudo bubbles contribute to the stability of the numerical method at their most (see Figure 5).
Finally, we recall that the pseudo bubble functions  *  ( = 1, 2, 3) are approximations to   on the subgrid specified above, through (19), and they are used in place of   to represent   .The approximate representation of   by bubble functions  *  ( = 1, 2, 3) is eventually used to solve (13) for its linear part.

Numerical Results
The methods discussed above and Link-Cutting Bubble (LCB) strategy proposed in [16] are compared for oneand two-dimensional convection-diffusion-reaction problems.This comparison is carried out by computing the numerical solutions for a wide range of problem parameters especially in the interesting case of small diffusion which corresponds to the convection-dominated or reactiondominated regimes.We also remark that the LCB method is employed only for one-dimensional problems as it is not designed for higher dimensions.

Numerical Results for 1D
Convection-Diffusion-Reaction Problems 6.1.1.Experiment 1.We first consider an equation with constant coefficients and compute the approximate solutions on uniform mesh.The uniform mesh is generated by dividing the unit interval [0, 1] into ten elements; that is, the mesh size ℎ = 1/10.In Figure 6, we present the linear parts of the numerical solutions together with the exact solution for   various intensities of reaction ( = 0.1, 1, 10, 20, 50, 100) when  = 10 −5 ,  = 1, and () = 1.The numerical solutions obtained with the PRFB and LCB methods are in excellent agreement with the exact solution for the whole problem parameters while the approximations generated by the SUPG and GLS methods possess spurious oscillations in layer regions, especially in reaction-dominated cases.We also note that there is no major difference between SUPG and GLS methods; however, the instabilities introduced by Galerkin are a little more amplified in GLS compared with SUPG.LCB methods are able to capture both the internal and boundary layers and produce approximations consistent with the physical solution while the approximations generated by the SUPG and GLS methods exhibit spurious oscillations in the reaction-dominated regime.
6.1.3.Experiment 3. Now, we consider a variable coefficient case with a turning point.We set  = 10 −5 ,  = −2(2−1), and () = 4(2 − 1) and decompose the domain into uniform discretization of 40 elements.Since there is no major difference between SUPG/GLS methods and PRFB/LCB methods, respectively, we only display the results with the SUPG and PRFB methods in Figure 8.The comparative study shows that the PRFB method works well even for turning point problems while the SUPG method exhibits slight oscillations around both the internal and the boundary layer.

Numerical Results for 2D
Convection-Diffusion-Reaction Problems 6.2.1.Experiment 1.We consider a test problem taken from [26].This problem is a benchmark problem, which exhibits numerical oscillations in both interior and boundary layer regions.The discontinuity in the boundary data causes interior layer that is propagated across unit square domain and there will also be layers at the outflow boundaries.Boundary conditions are displayed in Figure 9.We first take a set of uniform triangular meshes which is made up of  = 10 elements in  and  directions, respectively.We set  = 10 −4 and  = 72 ∘ and present elevation and contour plots of the solutions obtained with the PRFB and SUPG methods for various values of reaction ( =  = 0.001, 1, 10, 20, 50, 10 2 , 10 3 , 10 4 , 10 5 ) in Figures 10 and 11.The performances of the SUPG, GLS, SGS, and the PRFB methods are displayed in Figures 12 and 13 when  = 20.The stabilized methods deliver similar results unless the reaction term is dominant.We also remark that, in reaction-dominated regimes, there is no significant difference between SUPG and GLS methods; in contrast, the approximations generated by the SGS method are a little more oscillatory.Although the PRFB method could not eliminate the oscillations at all, it captures the layers well like in the RFB method (see [11,12,26,27] and references therein) in any regimes.
6.2.2.Experiment 2. Now, we consider a problem with convection skew to the discretization.We display the boundary conditions and the triangular elements used in discretization of the problem domain in Figure 14.In Figure 15, we present the solutions obtained with the SUPG, GLS, SGS, and PRFB methods when  = 10 −4 ,  = (0.15, 0.1), and  = 0 for various values of reaction ( = 10, 10 2 , 10 3 ).The numerical solutions in Figure 15 show that the PRFB method is more robust as both results are consistent with the physical configuration of the problem even on nonuniform meshes while the classical stabilization techniques, SUPG, GLS, and SGS, present similar qualitative behavior.

Experiment 3:
Thermal Boundary Layer.Finally, we consider the problem statement shown in Figure 16 with the following boundary conditions: The flow is taken as  = (2, 0).This problem may be viewed as the simulation of the development of a thermal boundary layer on a fully developed flow between two parallel plates, where the top plate is moving with velocity equal to one and the bottom plate is fixed.We take a set of uniform triangular meshes which is made up of  = 10 elements in  and  directions, respectively (see Figure 16).We set the  diffusion  = 10 −5 and the reaction term  = 10 −3 , 10 −1 , 5 × 10 −1 , 1, 5, 10, 10 2 , 10 3 , 10 4 .The numerical results and the corresponding contour plots are reported in Figures 17 and  18.The plots show that the PRFB solution captures the characteristic features of the exact solution even on coarse meshes.The SUPG method is again more oscillatory than the PRFB method.The results for SGS and GLS methods present a similar qualitative behavior as in Experiment 1 and thus are not shown.

Conclusion
In this paper, several benchmark problems exhibiting boundary/internal layers are used to test the performance and the  robustness of the stabilized FEMs for the numerical solution of the convection-diffusion-reaction problems.We observe that the classical SUPG method is very similar to the GLS method, whereas the PRFB method has a stabilization effect similar to LCB for one-dimensional convection-diffusionreaction problems.We also note that the numerical results obtained with the SUPG, GLS, SGS, and the PRFB methods are comparable for the convection-dominated regimes.However, as the problem turns into reaction-dominated case, the PRFB approach achieves a better performance than the SUPG, GLS, and SGS methods.Therefore the PRFB strategy can be viewed as a promising tool to solve both      one-and two-dimensional convection-diffusion-reaction problems.

a
different stabilized method and some popular examples to the choice of L appearing in the literature are SUPG : LV =  ⋅ ∇V, GLS : LV = −ΔV +  ⋅ ∇V + V, SGS : LV = ΔV +  ⋅ ∇V − V.

Figure 1 :
Figure 1: Splitting  into three subregions by   (figure reproduced from [15] [under the Creative Commons Attribution License/public domain]).

Figure 2 :Figure 3 :
Figure 2: Configuration of internal nodes for convection-dominated regime: one inflow edge (a) and two inflow edges (b) (figure reproduced from [15] [under the Creative Commons Attribution License/public domain]).

Figure 4 :Figure 5 :
Figure 4: Configuration of internal nodes for reaction-dominated regime: one inflow edge (a) and two inflow edges (b) (figure reproduced from [15] [under the Creative Commons Attribution License/public domain]).

Figure 6 :
Figure 6: The numerical solutions and the exact solution for several values of  when () = 1 and  = 10 −5 .

1Figure 14 :
Figure 14: (a) Configuration of Experiment 2. (b) Triangular elements used in discretization of the problem domain.

Figure 16 :Figure 17 :
Figure 16: (a) Configuration of Experiment 3. (b) Triangular elements used in the discretization of the problem domain.