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

Authors: Ali Sendur
Journal: Journal of Applied Mathematics
 (2018)

---

## Abstract

The disproportionality in the problem parameters of the convection-diffusion-reaction equation may lead to the formation of layer structures in some parts of the problem domain which are difficult to resolve by the standard numerical algorithms. Therefore the use of a stabilized numerical method is inevitable. In this work, we employ and compare three classical stabilized finite element formulations, namely, the Streamline-Upwind Petrov-Galerkin (SUPG), Galerkin/Least-Squares (GLS), and Subgrid Scale (SGS) methods, and a recent Link-Cutting Bubble (LCB) strategy proposed by Brezzi and his coworkers for the numerical solution of the convection-diffusion-reaction equation, especially in the case of small diffusion. On the other hand, we also consider the pseudo residual-free bubble (PRFB) method as another alternative that is based on enlarging the finite element space by a set of appropriate enriching functions. We compare the performances of these stabilized methods on several benchmark problems. Numerical experiments show that the proposed methods are comparable and display good performance, especially in the convection-dominated regime. However, as the problem turns into reaction-dominated case, the PRFB method is slightly better than the other well-known and extensively used stabilized finite element formulations as they start to exhibit oscillations.


---

## Body

## 1. 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-diffusion-reaction (CDR) problems [1–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–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–15]. In particular, their location is determined by minimizing the residual of a local differential problem defining the bubbles, with respect to the L1-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 convection-diffusion-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 Section2, 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.


## 2. Statement of the Problem

The convection-diffusion-reaction problem in two space dimensions consists of finding the scalar-valued functionu in a polygonal domain Ω with a smooth boundary Γ=∂Ω and Dirichlet boundary conditions, such that (1)Lu=-ϵΔu+β·∇u+σu=finΩ,u=0onΓ,where ϵ>0 is the diffusion coefficient, β is the velocity field, σ>0 is the reaction coefficient, and f is a smooth source function in L2(Ω). 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 u∈H01(Ω) such that (2)au,v=f,v∀v∈H01Ω,where (3)au,v=ϵ∫Ω∇u·∇v+∫Ωβ·∇uv+∫Ωσuvis a continuous and coercive bilinear form on the Hilbert space H01(Ω). Let us consider a partition Th of the computational domain Ω into triangular elements K that do not overlap and let hK=diam(K) with h=maxK∈ThhK. If we restrict ourselves to the case of continuous, piecewise linear elements which we denote by (4)VL=v∈H01Ω:vKis  linear  forK∈Th,then the standard Galerkin finite element method reads as follows: Find uL∈VL such that (5)auL,vL=f,vL∀vL∈VL.It is well-known that if the problem is convection or reaction-dominated, then, unless h 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.


## 3. 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 uL∈VL such that (6)auL,vL+SuL,vL=f,vL∀vL∈VL,where S(·,·) indicates a consistent stabilizing term added to the standard Galerkin formulation (5). The general form of the operator S(·,·) is given by (7)SuL,vL=∑KτKLuL-f,L~vLK,where τK is the stabilization parameter that needs to be tuned up by user and (·,·)K denotes the restriction of the inner product over K. Each choice of the operator L~ gives rise to a different stabilized method and some popular examples to the choice of L~ appearing in the literature are (8)SUPG:L~v=β·∇v,GLS:L~v=-ϵΔv+β·∇v+σv,SGS:L~v=ϵΔv+β·∇v-σv.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–23].One of the important points in the implementation of the stabilized methods is the choice of the stabilization parameterτK. The parameter τK 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 τK in the literature [21, 24, 25] are as follows: (9)τKc=14ϵ/h2+2β/h+σ,τKs=94ϵh22+2βh2+σ2-1/2,τKa=112ϵ/h2+2β/h+2σ,τKfv=16ϵ/hK2ζPe2+σζPe1,where (10)ζx=1,0≤x≤1;x,x≥1,Pe1=6ϵσhK2Pe2=βhK3ϵ.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 τKc 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.


## 4. Stabilization through Augmented Spaces

Another successful class of the stabilized finite element methods is based on enlarging the polynomial finite element spaceVL 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 Vh=VL⊕VB where the enriching space VB is defined by (11)VB=⨁K∈ThBKwithBK=H01Kand solving problem (5) on Vh. In this approach, it is required that the bubble component uB of uh satisfies the original differential equation in each K strongly; that is, (12)LuB=-LuL+finKsubject  to  the  boundary  conditionuB=0on∂K.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: Find uh=uL+uB in Vh such that (13)auL,vL+auB,vL=f,vL∀vL∈VL.The term a(uB,vL) is responsible for the stabilization of the numerical method and the bubble component uB should be computed before we solve (13) for its linear part uL. However, since uB is identified by the linear part uL and the source function f 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 touB 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.


## 5. 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 ofK by Vi, the edges of K by ei opposite to Vi, the length of ei by ei, the midpoint of edge ei by Mi, the outward unit normal to ei by ni, νi=eini, βνi=(β,νi), and the area of K by K (see Figure 1). The actual numbering of the vertices will be chosen according to the direction of β.Figure 1
SplittingK into three subregions by Pi (figure reproduced from [15] [under the Creative Commons Attribution License/public domain]).Now let us assume that the location ofPi is along the median from Vi; that is, (14)Pi=1-tiMi+tiVi,0<ti<1,i=1,2,3.To specify the problem regimes, we further introduce (15)ϵi∗=2K-3βνi+σK9e12+e22+e32i=1,2,3.Consider the bubble functions Bi, i=1,2,3, defined by (16)LBi=-LψiinK,Bi=0on⁡∂K,i=1,2,3,where ψi are the restrictions of the piecewise linear basis functions for VL to K and(17)LBf=finK,Bf=0on∂K.Now if (18)uLK=∑i=13ciψithenuBK=∑i=13ciBi+Bf,where  ci∈R,  i=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 Bi∗, that have the same qualitative behavior with their continuous counterpart. To start with, let bi be a piecewise linear function with bi(Vj)=0 and biPi=1∀i,j=1,2,3  and Bi∗(x)=αibi(x) be the classical Galerkin approximation of Bi through (16); that is, (19)aBi∗,biK=-Lψi,biK,i=1,2,3.Here the choice of the internal point Pi is critical and the main criteria used to determine its location are to minimize L1 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 Pi such that (20)Ji=∫KLBi∗+Lψi,i=1,2,3is minimum. The criteria (20) are employed to determine the locations of subgrid nodes Pi in convection and reaction-dominated cases separately.


### 5.1. Convection-Dominated Flows

In order to present the explicit descriptions ofPi, 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).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]).
(a) (b)


#### 5.1.1. One Inflow Edge

The problem is assumed to be convection-dominated if(21)ϵ≤ϵ1∗withσK<3min⁡βν2,βν3.From a previous calculation [15], the explicit location of P1 along the median from V1 is given by (22)t1∗=1--ρ1+ρ12+λ12σK2,where ρ1=-2βν1K+3ϵe2-e32 and λ1=24ϵσK2(e22+e32). The choice of other points P2 and P3 should be consistent with the physics of the problem. Thus we take t1=t1∗,t2=t3=min{1/3,1-t1∗}.


#### 5.1.2. Two Inflow Edges

The problem is assumed to be convection-dominated if(23)ϵ≤min⁡ϵ2∗,ϵ3∗withσK<3βν1.From [15], the position of P2 along the median from V2 is given by (24)t2∗=1--ρ2+ρ22+λ22σK2,where ρ2=-2βν2K+3ϵe1-e32 and λ2=24ϵσK2(e12+e32), and the position of P3 along the median from V3 is given by (25)t3∗=1--ρ3+ρ32+λ32σK2,where ρ3=-2βν3K+3ϵe1-e22 and λ3=24ϵσK2(e12+e22). The choice of P1 should be consistent with the physics of the problem. Thus we take t2=t2∗,t3=t3∗,t1=min{1/3,1-t2∗,1-t3∗}.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 convection-dominated one, so that the pseudo bubbles contribute to the stability of the numerical method at their most (see Figure3).Figure 3
Bubble functions,b1, b2, b3, in a typical element having (a) one inflow edge and (b) two inflow edges, when β=(sin⁡72°,cos⁡72°), σ=0.001, ϵ=10-2,10-3,10-4, and N=10 (figure reproduced from [15] [under the Creative Commons Attribution License/public domain]).
(a) (b)


### 5.2. Reaction-Dominated Flows

In order to present the explicit descriptions ofPi 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).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]).
(a) (b)


#### 5.2.1. One Inflow Edge

The problem is assumed to be reaction-dominated if(26)ϵ≤ϵ1∗withσK>3min⁡βν2,βν3.The location of P1 is taken from (22) and the position of P2 along the median from V2 is given by (27)t2∗∗=1-ρ2+ρ22+λ22σK2,where ρ2=2βν2K-3ϵe1-e32 and λ2=24ϵσK2(e12+e32), and the position of P3 along the median from V3 is given by (28)t3∗∗=1-ρ3+ρ32+λ32σK2,where ρ3=2βν3K-3ϵe1-e22 and λ3=24ϵσK2(e12+e22) (see [15]). Thus the positions of Pi, i=1,2,3, are determined by taking t1=t1∗,t2=max{1-t1,t2∗∗} and t3=max{1-t1,t3∗∗}.


#### 5.2.2. Two Inflow Edges

In this case, the problem is assumed to be reaction-dominated if(29)ϵ≤min⁡ϵ2∗,ϵ3∗withσK>3βν1and the locations of P2 and P3 are taken from (24)-(25). Now, the position of P1 along the median from V1 is given by (30)t1∗∗=1-ρ1+ρ12+λ12σK2,where ρ1=2βν1K-3ϵe2-e32 and λ1=24ϵσK2(e22+e32). Thus the positions of Pi, i=1,2,3, are determined by taking t1=max⁡{min⁡{1/3,1-t2,1-t3},t1∗∗}, t2=t2∗, t3=t3∗.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 reaction-dominated one, so that the pseudo bubbles contribute to the stability of the numerical method at their most (see Figure5).Figure 5
Bubble functions,b1, b2, b3, in a typical element having (a) two inflow edges and (b) one inflow edge, when ϵ=10-3, β=(sin⁡72°,cos⁡72°), σ=10,100,500, and N=10 (figure reproduced from [15] [under the Creative Commons Attribution License/public domain]).
(a) (b)Finally, we recall that the pseudo bubble functionsBi∗ (i=1,2,3) are approximations to Bi on the subgrid specified above, through (19), and they are used in place of Bi to represent uB. The approximate representation of uB by bubble functions Bi∗ (i=1,2,3) is eventually used to solve (13) for its linear part.


## 5.1. Convection-Dominated Flows

In order to present the explicit descriptions ofPi, 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).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]).
(a) (b)


### 5.1.1. One Inflow Edge

The problem is assumed to be convection-dominated if(21)ϵ≤ϵ1∗withσK<3min⁡βν2,βν3.From a previous calculation [15], the explicit location of P1 along the median from V1 is given by (22)t1∗=1--ρ1+ρ12+λ12σK2,where ρ1=-2βν1K+3ϵe2-e32 and λ1=24ϵσK2(e22+e32). The choice of other points P2 and P3 should be consistent with the physics of the problem. Thus we take t1=t1∗,t2=t3=min{1/3,1-t1∗}.


### 5.1.2. Two Inflow Edges

The problem is assumed to be convection-dominated if(23)ϵ≤min⁡ϵ2∗,ϵ3∗withσK<3βν1.From [15], the position of P2 along the median from V2 is given by (24)t2∗=1--ρ2+ρ22+λ22σK2,where ρ2=-2βν2K+3ϵe1-e32 and λ2=24ϵσK2(e12+e32), and the position of P3 along the median from V3 is given by (25)t3∗=1--ρ3+ρ32+λ32σK2,where ρ3=-2βν3K+3ϵe1-e22 and λ3=24ϵσK2(e12+e22). The choice of P1 should be consistent with the physics of the problem. Thus we take t2=t2∗,t3=t3∗,t1=min{1/3,1-t2∗,1-t3∗}.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 convection-dominated one, so that the pseudo bubbles contribute to the stability of the numerical method at their most (see Figure3).Figure 3
Bubble functions,b1, b2, b3, in a typical element having (a) one inflow edge and (b) two inflow edges, when β=(sin⁡72°,cos⁡72°), σ=0.001, ϵ=10-2,10-3,10-4, and N=10 (figure reproduced from [15] [under the Creative Commons Attribution License/public domain]).
(a) (b)


## 5.1.1. One Inflow Edge

The problem is assumed to be convection-dominated if(21)ϵ≤ϵ1∗withσK<3min⁡βν2,βν3.From a previous calculation [15], the explicit location of P1 along the median from V1 is given by (22)t1∗=1--ρ1+ρ12+λ12σK2,where ρ1=-2βν1K+3ϵe2-e32 and λ1=24ϵσK2(e22+e32). The choice of other points P2 and P3 should be consistent with the physics of the problem. Thus we take t1=t1∗,t2=t3=min{1/3,1-t1∗}.


## 5.1.2. Two Inflow Edges

The problem is assumed to be convection-dominated if(23)ϵ≤min⁡ϵ2∗,ϵ3∗withσK<3βν1.From [15], the position of P2 along the median from V2 is given by (24)t2∗=1--ρ2+ρ22+λ22σK2,where ρ2=-2βν2K+3ϵe1-e32 and λ2=24ϵσK2(e12+e32), and the position of P3 along the median from V3 is given by (25)t3∗=1--ρ3+ρ32+λ32σK2,where ρ3=-2βν3K+3ϵe1-e22 and λ3=24ϵσK2(e12+e22). The choice of P1 should be consistent with the physics of the problem. Thus we take t2=t2∗,t3=t3∗,t1=min{1/3,1-t2∗,1-t3∗}.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 convection-dominated one, so that the pseudo bubbles contribute to the stability of the numerical method at their most (see Figure3).Figure 3
Bubble functions,b1, b2, b3, in a typical element having (a) one inflow edge and (b) two inflow edges, when β=(sin⁡72°,cos⁡72°), σ=0.001, ϵ=10-2,10-3,10-4, and N=10 (figure reproduced from [15] [under the Creative Commons Attribution License/public domain]).
(a) (b)


## 5.2. Reaction-Dominated Flows

In order to present the explicit descriptions ofPi 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).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]).
(a) (b)


### 5.2.1. One Inflow Edge

The problem is assumed to be reaction-dominated if(26)ϵ≤ϵ1∗withσK>3min⁡βν2,βν3.The location of P1 is taken from (22) and the position of P2 along the median from V2 is given by (27)t2∗∗=1-ρ2+ρ22+λ22σK2,where ρ2=2βν2K-3ϵe1-e32 and λ2=24ϵσK2(e12+e32), and the position of P3 along the median from V3 is given by (28)t3∗∗=1-ρ3+ρ32+λ32σK2,where ρ3=2βν3K-3ϵe1-e22 and λ3=24ϵσK2(e12+e22) (see [15]). Thus the positions of Pi, i=1,2,3, are determined by taking t1=t1∗,t2=max{1-t1,t2∗∗} and t3=max{1-t1,t3∗∗}.


### 5.2.2. Two Inflow Edges

In this case, the problem is assumed to be reaction-dominated if(29)ϵ≤min⁡ϵ2∗,ϵ3∗withσK>3βν1and the locations of P2 and P3 are taken from (24)-(25). Now, the position of P1 along the median from V1 is given by (30)t1∗∗=1-ρ1+ρ12+λ12σK2,where ρ1=2βν1K-3ϵe2-e32 and λ1=24ϵσK2(e22+e32). Thus the positions of Pi, i=1,2,3, are determined by taking t1=max⁡{min⁡{1/3,1-t2,1-t3},t1∗∗}, t2=t2∗, t3=t3∗.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 reaction-dominated one, so that the pseudo bubbles contribute to the stability of the numerical method at their most (see Figure5).Figure 5
Bubble functions,b1, b2, b3, in a typical element having (a) two inflow edges and (b) one inflow edge, when ϵ=10-3, β=(sin⁡72°,cos⁡72°), σ=10,100,500, and N=10 (figure reproduced from [15] [under the Creative Commons Attribution License/public domain]).
(a) (b)Finally, we recall that the pseudo bubble functionsBi∗ (i=1,2,3) are approximations to Bi on the subgrid specified above, through (19), and they are used in place of Bi to represent uB. The approximate representation of uB by bubble functions Bi∗ (i=1,2,3) is eventually used to solve (13) for its linear part.


## 5.2.1. One Inflow Edge

The problem is assumed to be reaction-dominated if(26)ϵ≤ϵ1∗withσK>3min⁡βν2,βν3.The location of P1 is taken from (22) and the position of P2 along the median from V2 is given by (27)t2∗∗=1-ρ2+ρ22+λ22σK2,where ρ2=2βν2K-3ϵe1-e32 and λ2=24ϵσK2(e12+e32), and the position of P3 along the median from V3 is given by (28)t3∗∗=1-ρ3+ρ32+λ32σK2,where ρ3=2βν3K-3ϵe1-e22 and λ3=24ϵσK2(e12+e22) (see [15]). Thus the positions of Pi, i=1,2,3, are determined by taking t1=t1∗,t2=max{1-t1,t2∗∗} and t3=max{1-t1,t3∗∗}.


## 5.2.2. Two Inflow Edges

In this case, the problem is assumed to be reaction-dominated if(29)ϵ≤min⁡ϵ2∗,ϵ3∗withσK>3βν1and the locations of P2 and P3 are taken from (24)-(25). Now, the position of P1 along the median from V1 is given by (30)t1∗∗=1-ρ1+ρ12+λ12σK2,where ρ1=2βν1K-3ϵe2-e32 and λ1=24ϵσK2(e22+e32). Thus the positions of Pi, i=1,2,3, are determined by taking t1=max⁡{min⁡{1/3,1-t2,1-t3},t1∗∗}, t2=t2∗, t3=t3∗.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 reaction-dominated one, so that the pseudo bubbles contribute to the stability of the numerical method at their most (see Figure5).Figure 5
Bubble functions,b1, b2, b3, in a typical element having (a) two inflow edges and (b) one inflow edge, when ϵ=10-3, β=(sin⁡72°,cos⁡72°), σ=10,100,500, and N=10 (figure reproduced from [15] [under the Creative Commons Attribution License/public domain]).
(a) (b)Finally, we recall that the pseudo bubble functionsBi∗ (i=1,2,3) are approximations to Bi on the subgrid specified above, through (19), and they are used in place of Bi to represent uB. The approximate representation of uB by bubble functions Bi∗ (i=1,2,3) is eventually used to solve (13) for its linear part.


## 6. Numerical Results

The methods discussed above and Link-Cutting Bubble (LCB) strategy proposed in [16] are compared for one- and 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 reaction-dominated regimes. We also remark that the LCB method is employed only for one-dimensional problems as it is not designed for higher dimensions.


### 6.1. 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 h=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 f(x)=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.Figure 6
The numerical solutions and the exact solution for several values ofσ when f(x)=1 and ϵ=10-5.
(a)
σ = 0.1 (b)
σ = 1 (c)
σ = 10 (d)
σ = 20 (e)
σ = 50 (f)
σ = 100


#### 6.1.2. Experiment  2

Next, we consider a test problem which exhibits both internal and boundary layers. External force is defined as(31)fx=50,x≤0.5,-50,x>0.5.In Figure 7, we present the linear parts of the numerical solutions for various intensities of reaction (σ=1,10,25,50) when ϵ=10-5, β=1, and h=1/40. The PRFB and 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.Figure 7
The numerical solutions for several values ofσ when ϵ=10-5 and β=1.
(a)
σ = 1 (b)
σ = 10 (c)
σ = 25 (d)
σ = 50


#### 6.1.3. Experiment  3

Now, we consider a variable coefficient case with a turning point. We setϵ=10-5, β=-2(2x-1), and f(x)=4(2x-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.Figure 8
The numerical solutions whenf(x)=4(2x-1), ϵ=10-5, and β=-2(2x-1).
(a)
σ = 0.1 (b)
σ = 1 (c)
σ = 10 (d)
σ = 20 (e)
σ = 50 (f)
σ = 500


### 6.2. 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 N=10 elements in x and y 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 (σ=f=0.001,1,10,20,50,102,103,104,105) in Figures 10 and 11. The performances of the SUPG, GLS, SGS, and the PRFB methods are displayed in Figures 12 and 13 when N=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.Figure 9
Configuration of Experiment  1.Figure 10
Numerical solutions forθ=72°, ϵ=10-4, and N=10.Figure 11
Corresponding contour plots.Figure 12
Numerical solutions forθ=72°, ϵ=10-4, and N=20.Figure 13
Corresponding contour plots.


#### 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 Figure14. In Figure 15, we present the solutions obtained with the SUPG, GLS, SGS, and PRFB methods when ϵ=10-4, β=(0.15,0.1), and f=0 for various values of reaction (σ=10,102,103). 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.Figure 14
(a) Configuration of Experiment  2. (b) Triangular elements used in discretization of the problem domain.
(a) (b)Figure 15
Numerical solution forϵ=10-4 (top to bottom: SUPG, P-RFB, GLS, and SGS; left to right: σ=10,102,103).


#### 6.2.3. Experiment  3: Thermal Boundary Layer

Finally, we consider the problem statement shown in Figure16 with the following boundary conditions: (32)ux,y=1,ifx=0,0≤y≤0.5,1,ify=0.5,0≤x≤1,0,ify=0,0≤x≤1,2y,ifx=1,0≤y≤0.5.Figure 16
(a) Configuration of Experiment  3. (b) Triangular elements used in the discretization of the problem domain.
(a) (b)The flow is taken asβ=(2y,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 N=10 elements in x and y 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,102,103,104. 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.Figure 17
Numerical approximations for Experiment  3 whenϵ=10-5, β=(2y,0), and N=10.Figure 18
The corresponding contour plots.


## 6.1. 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 h=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 f(x)=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.Figure 6
The numerical solutions and the exact solution for several values ofσ when f(x)=1 and ϵ=10-5.
(a)
σ = 0.1 (b)
σ = 1 (c)
σ = 10 (d)
σ = 20 (e)
σ = 50 (f)
σ = 100


### 6.1.2. Experiment  2

Next, we consider a test problem which exhibits both internal and boundary layers. External force is defined as(31)fx=50,x≤0.5,-50,x>0.5.In Figure 7, we present the linear parts of the numerical solutions for various intensities of reaction (σ=1,10,25,50) when ϵ=10-5, β=1, and h=1/40. The PRFB and 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.Figure 7
The numerical solutions for several values ofσ when ϵ=10-5 and β=1.
(a)
σ = 1 (b)
σ = 10 (c)
σ = 25 (d)
σ = 50


### 6.1.3. Experiment  3

Now, we consider a variable coefficient case with a turning point. We setϵ=10-5, β=-2(2x-1), and f(x)=4(2x-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.Figure 8
The numerical solutions whenf(x)=4(2x-1), ϵ=10-5, and β=-2(2x-1).
(a)
σ = 0.1 (b)
σ = 1 (c)
σ = 10 (d)
σ = 20 (e)
σ = 50 (f)
σ = 500


## 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 h=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 f(x)=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.Figure 6
The numerical solutions and the exact solution for several values ofσ when f(x)=1 and ϵ=10-5.
(a)
σ = 0.1 (b)
σ = 1 (c)
σ = 10 (d)
σ = 20 (e)
σ = 50 (f)
σ = 100


## 6.1.2. Experiment  2

Next, we consider a test problem which exhibits both internal and boundary layers. External force is defined as(31)fx=50,x≤0.5,-50,x>0.5.In Figure 7, we present the linear parts of the numerical solutions for various intensities of reaction (σ=1,10,25,50) when ϵ=10-5, β=1, and h=1/40. The PRFB and 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.Figure 7
The numerical solutions for several values ofσ when ϵ=10-5 and β=1.
(a)
σ = 1 (b)
σ = 10 (c)
σ = 25 (d)
σ = 50


## 6.1.3. Experiment  3

Now, we consider a variable coefficient case with a turning point. We setϵ=10-5, β=-2(2x-1), and f(x)=4(2x-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.Figure 8
The numerical solutions whenf(x)=4(2x-1), ϵ=10-5, and β=-2(2x-1).
(a)
σ = 0.1 (b)
σ = 1 (c)
σ = 10 (d)
σ = 20 (e)
σ = 50 (f)
σ = 500


## 6.2. 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 N=10 elements in x and y 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 (σ=f=0.001,1,10,20,50,102,103,104,105) in Figures 10 and 11. The performances of the SUPG, GLS, SGS, and the PRFB methods are displayed in Figures 12 and 13 when N=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.Figure 9
Configuration of Experiment  1.Figure 10
Numerical solutions forθ=72°, ϵ=10-4, and N=10.Figure 11
Corresponding contour plots.Figure 12
Numerical solutions forθ=72°, ϵ=10-4, and N=20.Figure 13
Corresponding contour plots.


### 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 Figure14. In Figure 15, we present the solutions obtained with the SUPG, GLS, SGS, and PRFB methods when ϵ=10-4, β=(0.15,0.1), and f=0 for various values of reaction (σ=10,102,103). 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.Figure 14
(a) Configuration of Experiment  2. (b) Triangular elements used in discretization of the problem domain.
(a) (b)Figure 15
Numerical solution forϵ=10-4 (top to bottom: SUPG, P-RFB, GLS, and SGS; left to right: σ=10,102,103).


### 6.2.3. Experiment  3: Thermal Boundary Layer

Finally, we consider the problem statement shown in Figure16 with the following boundary conditions: (32)ux,y=1,ifx=0,0≤y≤0.5,1,ify=0.5,0≤x≤1,0,ify=0,0≤x≤1,2y,ifx=1,0≤y≤0.5.Figure 16
(a) Configuration of Experiment  3. (b) Triangular elements used in the discretization of the problem domain.
(a) (b)The flow is taken asβ=(2y,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 N=10 elements in x and y 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,102,103,104. 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.Figure 17
Numerical approximations for Experiment  3 whenϵ=10-5, β=(2y,0), and N=10.Figure 18
The corresponding contour plots.


## 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 N=10 elements in x and y 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 (σ=f=0.001,1,10,20,50,102,103,104,105) in Figures 10 and 11. The performances of the SUPG, GLS, SGS, and the PRFB methods are displayed in Figures 12 and 13 when N=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.Figure 9
Configuration of Experiment  1.Figure 10
Numerical solutions forθ=72°, ϵ=10-4, and N=10.Figure 11
Corresponding contour plots.Figure 12
Numerical solutions forθ=72°, ϵ=10-4, and N=20.Figure 13
Corresponding contour plots.


## 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 Figure14. In Figure 15, we present the solutions obtained with the SUPG, GLS, SGS, and PRFB methods when ϵ=10-4, β=(0.15,0.1), and f=0 for various values of reaction (σ=10,102,103). 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.Figure 14
(a) Configuration of Experiment  2. (b) Triangular elements used in discretization of the problem domain.
(a) (b)Figure 15
Numerical solution forϵ=10-4 (top to bottom: SUPG, P-RFB, GLS, and SGS; left to right: σ=10,102,103).


## 6.2.3. Experiment  3: Thermal Boundary Layer

Finally, we consider the problem statement shown in Figure16 with the following boundary conditions: (32)ux,y=1,ifx=0,0≤y≤0.5,1,ify=0.5,0≤x≤1,0,ify=0,0≤x≤1,2y,ifx=1,0≤y≤0.5.Figure 16
(a) Configuration of Experiment  3. (b) Triangular elements used in the discretization of the problem domain.
(a) (b)The flow is taken asβ=(2y,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 N=10 elements in x and y 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,102,103,104. 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.Figure 17
Numerical approximations for Experiment  3 whenϵ=10-5, β=(2y,0), and N=10.Figure 18
The corresponding contour plots.


## 7. 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-diffusion-reaction 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.



---

*Source: 4259634-2018-03-26.xml*