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.
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 -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 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.
2. 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 is the diffusion coefficient, is the velocity field, is the reaction coefficient, and is a smooth source function in . 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 such that where is a continuous and coercive bilinear form on the Hilbert space . Let us consider a partition of the computational domain into triangular elements that do not overlap and let with . 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 reaction-dominated, 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.
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 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 gives rise to a different stabilized method and some popular examples to the choice of appearing in the literature are
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 . 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.
4. 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: Find in such that The term 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.
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 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 , , defined by where are the restrictions of the piecewise linear basis functions for to andNow if where , . 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 and 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 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 reaction-dominated cases separately.
5.1. 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).
(a)
(b)
5.1.1. One Inflow Edge
The problem is assumed to be convection-dominated if From a previous calculation [15], the explicit location of along the median from is given by where and . The choice of other points and should be consistent with the physics of the problem. Thus we take
5.1.2. Two Inflow Edges
The problem is assumed to be convection-dominated if From [15], the position of along the median from is given by where and , and the position of along the median from is given by where and . The choice of 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 convection-dominated one, so that the pseudo bubbles contribute to the stability of the numerical method at their most (see Figure 3).
(a)
(b)
5.2. 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).
(a)
(b)
5.2.1. One Inflow Edge
The problem is assumed to be reaction-dominated if The location of is taken from (22) and the position of along the median from is given by where and , and the position of along the median from is given by where and (see [15]). Thus the positions of , , are determined by taking and
5.2.2. Two Inflow Edges
In this case, the problem is assumed to be reaction-dominated if and the locations of and are taken from (24)-(25). Now, the position of along the median from is given by where and Thus the positions of , , 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 reaction-dominated one, so that the pseudo bubbles contribute to the stability of the numerical method at their most (see Figure 5).
(a)
(b)
Finally, we recall that the pseudo bubble functions () 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 () 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 into ten elements; that is, the mesh size . In Figure 6, we present the linear parts of the numerical solutions together with the exact solution for various intensities of reaction () when , , and . 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.
(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 asIn Figure 7, we present the linear parts of the numerical solutions for various intensities of reaction () when , , and . 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.
(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 , , and 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.
(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 elements in and directions, respectively. We set and and present elevation and contour plots of the solutions obtained with the PRFB and SUPG methods for various values of reaction () in Figures 10 and 11. The performances of the SUPG, GLS, SGS, and the PRFB methods are displayed in Figures 12 and 13 when . 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 , , and for various values of reaction (). 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.
(a)
(b)
6.2.3. Experiment 3: Thermal Boundary Layer
Finally, we consider the problem statement shown in Figure 16 with the following boundary conditions:
(a)
(b)
The flow is taken as . 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 elements in and directions, respectively (see Figure 16). We set the diffusion and the reaction term . 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.
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.
Conflicts of Interest
The author declares that there are no conflicts of interest regarding the publication of this paper.