Journal of Applied Mathematics

Volume 2018, Article ID 4259634, 16 pages

https://doi.org/10.1155/2018/4259634

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

Alanya Alaaddin Keykubat University, 07490 Antalya, Turkey

Correspondence should be addressed to Ali Sendur; rt.ude.aynala@rudnes.ila

Received 12 November 2017; Revised 31 January 2018; Accepted 14 February 2018; Published 26 March 2018

Academic Editor: Srinivasan Natesan

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

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