# A Posteriori Error Estimates for Hughes Stabilized SUPG Technique and Adaptive Refinement for a Convection-Diffusion Problem

Authors: Brehmit Kaur; Vivek Sangwan
Journal: Advances in Mathematical Physics
 (2020)

---

## Abstract

The motive of the present work is to propose an adaptive numerical technique for singularly perturbed convection-diffusion problem in two dimensions. It has been observed that for small singular perturbation parameter, the problem under consideration displays sharp interior or boundary layers in the solution which cannot be captured by standard numerical techniques. In the present work, Hughes stabilization strategy along with the streamline upwind/Petrov-Galerkin (SUPG) method has been proposed to capture these boundary layers. Reliable a posteriori error estimates in energy norm on anisotropic meshes have been developed for the proposed scheme. But these estimates prove to be dependent on the singular perturbation parameter. Therefore, to overcome the difficulty of oscillations in the solution, an efficient adaptive mesh refinement algorithm has been proposed. Numerical experiments have been performed to test the efficiency of the proposed algorithm.


---

## Body

## 1. Introduction

Singularly perturbed problems occur frequently in various branches of applied science and engineering, e.g., fluid dynamics, aerodynamics, oceanography, quantum mechanics, chemical reactor theory, reaction-diffusion processes, and radiating flows. In general, it has been observed that singularly perturbed problems exhibit singularities as the singular perturbation parameterε⟶0. Therefore, it becomes essential to implement some robust numerical technique to capture these singularities. In literature, there exist various numerical techniques to handle these singularities. Adaptive mesh refinement techniques are one of such techniques. Very few researchers have proposed adaptive refinement strategies for singularly perturbed convection-diffusion problems. Generally, adaptive refinement techniques are based on two types of error estimates, namely, a priori error estimates and a posteriori error estimates. Nicaise [1] developed a posteriori residual error estimates for convection-diffusion-reaction problems using some cell-centered finite volume methods. Based on a posteriori error estimates, the author proposed an adaptive algorithm. John [2, 3] did numerical study of various a posteriori error estimates and indicators for convection-diffusion problems. On the basis of error estimates, the author proposed numerical solution of singularly perturbed convection-dominated problems on adaptive refined grid. Repin and Nicaise [4] derived a posteriori error estimates for linear convection-diffusion-reaction problems using functional arguments. Verfurth [5] derived a posteriori error estimates for convection-dominated stationary convection-diffusion equation using locally refined isotropic meshes. Zhao et al. [6] proposed adaptive numerical technique for convection-diffusion equations based on semirobust residual a posteriori error estimates for lower order nonconforming finite element approximations of streamline diffusion method.From literature, we know that the classical finite element methods [7] fail to provide satisfactory results for small values of singular perturbation parameter, i.e., when ε⟶0. It has also been observed that streamline upwind/Petrov-Galerkin (SUPG) method provides good approximate solution in the region where there is no sharp change in the solution but fails in the small subregions of sharp boundary layers. It has been observed that occurrence of these nonphysical oscillations in the region of sharp boundary layers in the discrete solution of SUPG method is based on the fact that this scheme is not monotonicity preserving. To overcome this difficulty, in the present work, we have proposed Hughes stabilization strategy [8] along with the SUPG method. It involves suitable addition of one more term which is multiple of a function in the direction where spurious oscillations were seen in approximate SUPG solution. This additional term is added on the left-hand side of SUPG discretization of convection-diffusion problem. The a posteriori error estimates have been derived for the proposed scheme. Based on these estimates, an anisotropic mesh refinement strategy has been proposed for singularly perturbed problems.The outline of the paper is as follows.In Section2, the continuous problem under consideration and its streamline upwind Petrov-Galerkin finite element approximation have been presented. In Section 3, some auxiliary tools which are required for deriving reliable error bounds have been presented. In Section 4, we have discussed residual-based a posteriori error estimates and derived error bounds on anisotropic meshes. An adaptive refinement algorithm based on derived a posteriori error estimates has been proposed in Section 5. Section 6 deals with some numerical experiments which have been performed to analyze the robustness and efficiency of the proposed adaptive refinement strategy. In the last section, conclusion has been presented.


## 2. Continuous Problem

Consider the following convection-diffusion equation in two dimensions:(1a)−∇·ε∇u+a·∇u+bu=finΩ,(1b)u=0on∂ΩD,(1c)ε∂u∂n=gon∂ΩN,where Ω⊂ℝ2 is a bounded domain with Lipschitz-continuous boundary ∂Ω=∂ΩD∪∂ΩN, ε (0<ε≪1) is singular perturbation parameter, and a, b, and f are sufficiently smooth. Here, ∂ΩD and ∂ΩN denote the Dirichlet and Neumann boundaries of the domain, respectively.W1,∞Ω and L∞Ω represent the usual Sobolev and Lebesgue spaces, respectively. The notation (.,.) has been used for inner product (.,.)Ω.Throughout the paper, we assume that−1/2∇·a+b≥b0>0.For any open bounded subsetK⊂Ω¯, let H1K be the standard Sobolev space. Further, we define
(2)V0=v∈H1Ω,v=0on∂ΩD.Let(3)vK2=ε∇vK2+b0vK2be energy norm on bounded subset K⊂Ω¯. The weak formulation of equations (1a), (1b) and (1c) is given by the following.Findu∈H1Ω such that
(4)Bu,v=F,vwhere
(5)Bu,v=ε∇u,∇v+a·∇u,v+bu,v,(6)F,v=f,v+g,v∂ΩN∀v∈V0.The existence and uniqueness of the solution of the above weak formulation (4) are guaranteed using Lax Milgram Lemma together with condition (1c). Let Γh be the admissible and shape-regular triangulation of domain Ω¯ consisting of triangles. Let L be any two-dimensional element with edge E. Let nL,E=nx,ny be unit outward normal vector to L along E (see Figure 1). Fixing one of the two normal vectors, let nE be the normal vector for each edge E.Figure 1
Orthogonality condition.It has been observed that the solution of singularly perturbed problem displays boundary layers if the Peclet number, as discussed below, is large. Define local mesh Peclet number as(7)PeK=a∞,KhminK2ε,where hminK is the minimal length of element K as defined in the next section.LetVh=vh∈H1Ω:vhK∈P1K, where P1K is the space of all linear polynomials over the element K and V0h=vh∈Vh:vh∂ΩD=0. Next, we discuss the SUPG method along with the Hughes stabilization technique for approximating the solution of problem (1a), (1b) and (1c). The SUPG method [9] for problem (1a), (1b) and (1c) is defined as follows.Finduh∈Vh such that
(8)Bρuh,vh=F,vh∀vh∈V0h,where Bρuh,vh=Buh,vh+Rhuh,ρa·∇hvh, Rhu=−εΔhu+a·∇hu+bu−f, ρ is nonnegative stabilization parameter, and Bu,v and F,v are defined in (5).


### 2.1. Hughes Stabilization Technique

It has been observed that the streamline upwind/Petrov-Galerkin (SUPG) method provides good approximate solution in the region where there is no sharp change in the solution but fails badly in the small subregions of sharp boundary layers. To overcome this difficulty, we use Hughes stabilization technique [10] to SUPG method. It involves an additional term Rhuh,σah·∇hvh in the left-hand side of SUPG finite element discretization of convection-diffusion equation where
(9)ah=a⋅∇uh∇uh∇uh2,if∇uh≠0,0,if∇uh=0,and σ is a nonnegative stabilization parameter. This additional term increases the robustness of SUPG method in the boundary layer region by controlling oscillations. Using Hughes stabilization technique to SUPG finite element method, equation (1a), (1b) and (1c) is discretized as follows.Finduh∈Vh such that
(10)Bρ,σuh,vh=F,vh∀vh∈V0h,where Bρ,σuh,vh=Buh,vh+Rhuh,ρa·∇hvh+Rhuh,σah·∇hvh and Rhu=−εΔhu+a·∇hu+bu−f.LetρK be the stabilization parameter over each element K. Ross et al. [7] showed that the approximate solution uh obtained using SUPG finite element discretization exists and is unique provided stabilization parameter ρK is small and satisfies
(11)0≤ρK≤12minb0b∞,K−2,hminK2ε−1ν−2,where hminK is the minimal length of element K and the constant ν satisfies the inequality
(12)∇·∇vhK≤νhminK−1∇vhK∀vh∈V0h.From inequality (12), it can easily be observed that ν=0 for piecewise linear functions in V0h. Therefore, the above bounds reduce to 0≤ρK≤b0/2b∞,K−2. In order to simplify the calculations, we introduce the notation c≲d which means that there exists a positive constant A independent of c, d, Γh, and ε such that c≤Ad. Further, we assume that
(13)ρK≲hminKa∞,K−1∀K∈Γh.Also, for any mesh functionvh∈V0h, using (12) and scaling arguments, we can get
(14)∇vhK≲hminK−1vhK.Using energy norm def. (3),
(15)vhK≤b0−1/2vhK.Thus, we have(16)∇vhK≲hminK−1b0−12vhK.Again, from energy norm, we have(17)∇vhK≤ε−1/2vhK.Using (16) and (17), we get
(18)∇vhK≲minhminK−1b0−1/2,ε−1/2vhK.


## 2.1. Hughes Stabilization Technique

It has been observed that the streamline upwind/Petrov-Galerkin (SUPG) method provides good approximate solution in the region where there is no sharp change in the solution but fails badly in the small subregions of sharp boundary layers. To overcome this difficulty, we use Hughes stabilization technique [10] to SUPG method. It involves an additional term Rhuh,σah·∇hvh in the left-hand side of SUPG finite element discretization of convection-diffusion equation where
(9)ah=a⋅∇uh∇uh∇uh2,if∇uh≠0,0,if∇uh=0,and σ is a nonnegative stabilization parameter. This additional term increases the robustness of SUPG method in the boundary layer region by controlling oscillations. Using Hughes stabilization technique to SUPG finite element method, equation (1a), (1b) and (1c) is discretized as follows.Finduh∈Vh such that
(10)Bρ,σuh,vh=F,vh∀vh∈V0h,where Bρ,σuh,vh=Buh,vh+Rhuh,ρa·∇hvh+Rhuh,σah·∇hvh and Rhu=−εΔhu+a·∇hu+bu−f.LetρK be the stabilization parameter over each element K. Ross et al. [7] showed that the approximate solution uh obtained using SUPG finite element discretization exists and is unique provided stabilization parameter ρK is small and satisfies
(11)0≤ρK≤12minb0b∞,K−2,hminK2ε−1ν−2,where hminK is the minimal length of element K and the constant ν satisfies the inequality
(12)∇·∇vhK≤νhminK−1∇vhK∀vh∈V0h.From inequality (12), it can easily be observed that ν=0 for piecewise linear functions in V0h. Therefore, the above bounds reduce to 0≤ρK≤b0/2b∞,K−2. In order to simplify the calculations, we introduce the notation c≲d which means that there exists a positive constant A independent of c, d, Γh, and ε such that c≤Ad. Further, we assume that
(13)ρK≲hminKa∞,K−1∀K∈Γh.Also, for any mesh functionvh∈V0h, using (12) and scaling arguments, we can get
(14)∇vhK≲hminK−1vhK.Using energy norm def. (3),
(15)vhK≤b0−1/2vhK.Thus, we have(16)∇vhK≲hminK−1b0−12vhK.Again, from energy norm, we have(17)∇vhK≤ε−1/2vhK.Using (16) and (17), we get
(18)∇vhK≲minhminK−1b0−1/2,ε−1/2vhK.


## 3. Some Important Notations and Tools

Since singularly perturbed convection-diffusion problems exhibit sharp boundary layers when Peclet number becomes large or the singular perturbation parameter becomes smaller, in such situations, elements with large aspect ratio (anisotropic meshes) are recommended. In this section, we will discuss some important results on anisotropic meshes.


### 3.1. Notations

Consider an arbitrary triangleK∈Γh with Q0Q1 as the longest edge (see Figure 2). Denote two orthogonal vectors qi,K with length hi,K=qi,K,i=1,2, where q1,K is taken along the largest edge Q0Q1. From Figure 2, it can be verified that h1,K≥h2,K. Define hminK=h2,K. These qi,K’s correspond to two anisotropic directions. Further, define an orthogonal matrix CK=q1,K,q2,K∈ℝ2×2. Let αK be the scaling factor defined as
(19)αK=minb0−1/2,ε−1/2·hminK.Figure 2
TriangleK.We represent triangles byK or K′ or Ki and its edges by E. Further define its height over edge E as
(20)hE,K=2·KE,where K represents the area of triangle K. Let wE be the bounded domain formed by using two triangles having common edge E. Further, define wK to be the domain consisting of triangle K and its edge neighboring triangles. Let
(21)PewK=maxK′⊂wKPeK′be mesh Peclet number on the domain wK where PeK is defined in (7). For an interior edge E=K1∩K2, define parameters hE=hE,K1+hE,K2/2, hminE=hminK1+hminK2/2, and αE=αK1+αK2/2.For boundary edgeE⊂∂K∩∂Ω, we define hE=hE,K, hminE=hminK, and αE=αK.Since the mesh considered is assumed to be shape-regular and admissible, along with these requirements, we take(22)hi,K~hi,K′∀K,K′withK∩K′≠∅,i=1,2,and the number of triangles with node yj is bounded uniformly.


### 3.2. Interpolation

In order to obtain reliable error upper bounds, we define a suitable matching function [11, 12] to measure the alignment of anisotropic mesh Γh and anisotropic function.Definition 1 (matching function).
Letu∈H1Ω and Γh∈F be the triangulation of Ω. We define M1:H1Ω×F⟶ℝ by
(23)M1u,Γh≔∑K∈ΓhhminK−2·CKT∇uK21/2∇u,where CK∈ℝ2×2 as defined earlier.We can easily verify thatM1u,Γh~1 for isotropic meshes. Similarly, it can easily be observed that M1u,Γh~1 for anisotropic meshes suitably aligned with anisotropic function u. Therefore, M1u,Γh≈C for anisotropic meshes.To propose reliable error estimates in energy norm, we will use Clément interpolation operatorRC [13] for u∈H1Ω as standard Lagrange interpolation cannot be defined for these functions.Lemma 2.
Letu∈H01Ω and αK be the scaling factor defined by (19). Then, the Clément interpolation operator RC:H01Ω↦V0h satisfies
(24)∑K∈ΓhαK−2·u−RcuK2≲C2u2,ε1/2∑E⊂Ω\∂ΩDαE−1·u−RcuE2≲C2u2.Proof.
The proof is discussed in [14].


## 3.1. Notations

Consider an arbitrary triangleK∈Γh with Q0Q1 as the longest edge (see Figure 2). Denote two orthogonal vectors qi,K with length hi,K=qi,K,i=1,2, where q1,K is taken along the largest edge Q0Q1. From Figure 2, it can be verified that h1,K≥h2,K. Define hminK=h2,K. These qi,K’s correspond to two anisotropic directions. Further, define an orthogonal matrix CK=q1,K,q2,K∈ℝ2×2. Let αK be the scaling factor defined as
(19)αK=minb0−1/2,ε−1/2·hminK.Figure 2
TriangleK.We represent triangles byK or K′ or Ki and its edges by E. Further define its height over edge E as
(20)hE,K=2·KE,where K represents the area of triangle K. Let wE be the bounded domain formed by using two triangles having common edge E. Further, define wK to be the domain consisting of triangle K and its edge neighboring triangles. Let
(21)PewK=maxK′⊂wKPeK′be mesh Peclet number on the domain wK where PeK is defined in (7). For an interior edge E=K1∩K2, define parameters hE=hE,K1+hE,K2/2, hminE=hminK1+hminK2/2, and αE=αK1+αK2/2.For boundary edgeE⊂∂K∩∂Ω, we define hE=hE,K, hminE=hminK, and αE=αK.Since the mesh considered is assumed to be shape-regular and admissible, along with these requirements, we take(22)hi,K~hi,K′∀K,K′withK∩K′≠∅,i=1,2,and the number of triangles with node yj is bounded uniformly.


## 3.2. Interpolation

In order to obtain reliable error upper bounds, we define a suitable matching function [11, 12] to measure the alignment of anisotropic mesh Γh and anisotropic function.Definition 1 (matching function).
Letu∈H1Ω and Γh∈F be the triangulation of Ω. We define M1:H1Ω×F⟶ℝ by
(23)M1u,Γh≔∑K∈ΓhhminK−2·CKT∇uK21/2∇u,where CK∈ℝ2×2 as defined earlier.We can easily verify thatM1u,Γh~1 for isotropic meshes. Similarly, it can easily be observed that M1u,Γh~1 for anisotropic meshes suitably aligned with anisotropic function u. Therefore, M1u,Γh≈C for anisotropic meshes.To propose reliable error estimates in energy norm, we will use Clément interpolation operatorRC [13] for u∈H1Ω as standard Lagrange interpolation cannot be defined for these functions.Lemma 2.
Letu∈H01Ω and αK be the scaling factor defined by (19). Then, the Clément interpolation operator RC:H01Ω↦V0h satisfies
(24)∑K∈ΓhαK−2·u−RcuK2≲C2u2,ε1/2∑E⊂Ω\∂ΩDαE−1·u−RcuE2≲C2u2.Proof.
The proof is discussed in [14].


## 4. Residual Error Estimates

In this section, firstly we discuss exact and approximate residuals. Further, we will develop reliable error upper bounds for Hughes stabilized SUPG finite element solution on anisotropic meshes. It is shown that the error bounds obtained depend on anisotropic interpolation.


### 4.1. Exact Residuals

We define exact element residualRK and exact edge residual RE as
(25)RK=f−−εΔvh+a·∇vh+bvhonK,REx=ε·lims→+0∂nEvhx+snE−∂nEvhx−snEifE⊂Ω\∂Ω,g−ε∂nvhifE⊂∂ΩN,0ifE⊂∂ΩD,where nE⊥E⊂Ω\∂Ω is the unitary normal vector and n⊥E⊂∂ΩN is the outer unitary normal vector.


### 4.2. Approximate Residuals

LetQ be the approximation operator used to approximate the element residual and the face residual, i.e.,
(26)rK=QRK∈P0K∀K∈Γh,rE=QRE∈P0E∀E,where we have denoted (approximate) element residual by rK and the (approximate) face residual by rE. Since the finite element solution vh is linear, thus
(27)rE=RE∀E⊂Ω\∂ΩN.


### 4.3. Residual Error Estimator

Residual error estimatorηK and the approximation term ζK over triangle K are given by
(28)ηK2=aK2·rKK2+ε−1/2·αK·∑E⊂∂K\∂ΩDrEE2,ζK2=aK2·rK−RKwK2+ε−1/2·αK·∑E⊂∂K∩∂ΩNrE−REE2,where αK is the scaling factor defined earlier and aK2=3αK2. Further, we define global error estimators as
(29)η2=∑K∈ΓhηK2,ζ2=∑K∈ΓhζK2.Next, we derive reliable upper error bounds on anisotropic meshes.Theorem 1 (residual error estimation).
Letv∈H01Ω be the exact solution and vh∈V0h be the finite element solution obtained by the proposed scheme. Then, the error in energy norm is bounded above globally by
(30)v−vh≲C∑K∈Γh3αK2rK−RKK2+rKK2+∑E⊂∂K\∂ΩDε−1/2αErE−REE2+rEE21/2.Proof.
We know thatBv,v≥v2∀v∈H01Ω.
Using this result, we get(31)v−vh≤Bv−vh,uv−vh=Bv−vh,uu,where u=v−vh. Introducing Clément interpolation operator RC, we can write the bilinear form B(.,.) as
(32)Bv−vh,u=Bv−vh,u−RCu+Bv−vh,RCu.Now, using the error equation and integration by parts, we have(33)Bv−vh,w=∑K∈ΓhRK,wK+∑E⊂Ω\∂ΩDRE,wE∀w∈H01Ω.Using equation (33), the middle term of equation (32) can be written as
(34)Bv−vh,u−RCu=∑K∈ΓhRK,u−RCuK+∑E⊂Ω\∂ΩDRE,u−RCuE.Using Cauchy Schwarz inequality, we get(35)∑K∈ΓhRK,u−RCuK≤∑K∈ΓhαK2RKK21/2·∑K∈ΓhαK−2u−RCuK21/2,∑E⊂Ω\∂ΩDRE,u−RCuE≤∑E⊂Ω\∂ΩDε−1/2αEREE21/2·∑E⊂Ω\∂ΩDε1/2αE−1u−RCuE21/2.Further, using Lemma2, we get
(36)∑K∈ΓhRK,u−RCuK≲∑K∈ΓhαK2RKK21/2Cu,∑E⊂Ω\∂ΩDRE,u−RCuE≲∑E⊂Ω\∂ΩDε−1/2αEREE21/2Cu.Therefore, the termBv−vh,u−RCu is bounded above by
(37)Bv−vh,u−RCu≲∑K∈ΓhαK2RKK2+∑E⊂Ω\∂ΩDε−1/2αEREE21/2Cu.From (18) and (19), we have
(38)∇uhK≲minhminK−1b0−1/2,ε−1/2uhK=hminK−1αKuhKforuh∈V0h.Next, we will find the bounds on the second termBv−vh,RCu of equation (32).
(39)Bv−vh,RCu=ε∇v,∇RCu+a·∇v,RCu+bv,RCu−ε∇vh,∇RCu+a·∇vh,RCu+bvh,RCu+∑KρKRK,a·∇RCu+∑KσKRK,ah·∇RCu.Using the standard Galerkin orthogonality condition and standard scaling results, the above equation reduces to(40)Bv−vh,RCu=−∑KρKRK,a·∇RCu−∑KσKRK,ah·∇RCu≤∑K∈ΓhρKRKKa∞,K∇RCuK+∑K∈ΓhσKRKKah∞,K∇RCuK≤∑K∈ΓhρKRKKa∞,KhminK−1αKRCuK+∑K∈ΓhσKRKKah∞,KhminK−1αKRCuK.We know that for Clément operator [11]
(41)RCu≲Cu∀u∈H01Ω.Thus, we have(42)Bv−vh,RCu≤∑K∈ΓhρK2RKK2a∞,K2hminK−2αK21/2Cu+∑K∈ΓhσK2RKK2ah∞,K2hminK−2αK21/2Cu.It may be noted that the effect of nonlinear term in theL∞ norm will be bounded by that of the term a∞ as shown below, i.e.,
(43)ah=a·∇uh∇uh∇uh2,∇uh≠0,ah≤a∇uh∇uh∇uh2UsingCauchy–Schwarzinequality,ah≤a∇uh∇uh∇uh2,≤a.Sincea, the convection coefficient, is assumed to be smooth in the domain under consideration, it is bounded above. Hence, the nonlinear term ah∞,K is taken as bounded above by some constant and is absorbed in the constant term.We know that(44)ρK≲hminK/a∞,K∀K∈Γh,σK≲hminK/ah∞,K∀K∈Γh.Therefore,(45)Bv−vh,RCu≲∑K∈ΓhαK2RKK2+∑K∈ΓhαK2RKK21/2·Cu·M1u,Γh≈C.Since(46)v−vh≤Bv−vh,uu,using equation (45) and equation (37) in equation (32), we get
(47)v−vh≲∑K∈Γh3αK2RKK2+∑E⊂Ω\∂ΩDε−1/2αEREE21/2C.Using triangle inequalities,(48)RKK2≤rK−RKK2+rKK2,REE2≤rE−REE2+rEE2,we get
(49)v−vh≲C∑K∈Γh3αK2rK−RKK2+rKK2+∑E⊂∂K\∂ΩDε−1/2αErE−REE2+rEE21/2.


## 4.1. Exact Residuals

We define exact element residualRK and exact edge residual RE as
(25)RK=f−−εΔvh+a·∇vh+bvhonK,REx=ε·lims→+0∂nEvhx+snE−∂nEvhx−snEifE⊂Ω\∂Ω,g−ε∂nvhifE⊂∂ΩN,0ifE⊂∂ΩD,where nE⊥E⊂Ω\∂Ω is the unitary normal vector and n⊥E⊂∂ΩN is the outer unitary normal vector.


## 4.2. Approximate Residuals

LetQ be the approximation operator used to approximate the element residual and the face residual, i.e.,
(26)rK=QRK∈P0K∀K∈Γh,rE=QRE∈P0E∀E,where we have denoted (approximate) element residual by rK and the (approximate) face residual by rE. Since the finite element solution vh is linear, thus
(27)rE=RE∀E⊂Ω\∂ΩN.


## 4.3. Residual Error Estimator

Residual error estimatorηK and the approximation term ζK over triangle K are given by
(28)ηK2=aK2·rKK2+ε−1/2·αK·∑E⊂∂K\∂ΩDrEE2,ζK2=aK2·rK−RKwK2+ε−1/2·αK·∑E⊂∂K∩∂ΩNrE−REE2,where αK is the scaling factor defined earlier and aK2=3αK2. Further, we define global error estimators as
(29)η2=∑K∈ΓhηK2,ζ2=∑K∈ΓhζK2.Next, we derive reliable upper error bounds on anisotropic meshes.Theorem 1 (residual error estimation).
Letv∈H01Ω be the exact solution and vh∈V0h be the finite element solution obtained by the proposed scheme. Then, the error in energy norm is bounded above globally by
(30)v−vh≲C∑K∈Γh3αK2rK−RKK2+rKK2+∑E⊂∂K\∂ΩDε−1/2αErE−REE2+rEE21/2.Proof.
We know thatBv,v≥v2∀v∈H01Ω.
Using this result, we get(31)v−vh≤Bv−vh,uv−vh=Bv−vh,uu,where u=v−vh. Introducing Clément interpolation operator RC, we can write the bilinear form B(.,.) as
(32)Bv−vh,u=Bv−vh,u−RCu+Bv−vh,RCu.Now, using the error equation and integration by parts, we have(33)Bv−vh,w=∑K∈ΓhRK,wK+∑E⊂Ω\∂ΩDRE,wE∀w∈H01Ω.Using equation (33), the middle term of equation (32) can be written as
(34)Bv−vh,u−RCu=∑K∈ΓhRK,u−RCuK+∑E⊂Ω\∂ΩDRE,u−RCuE.Using Cauchy Schwarz inequality, we get(35)∑K∈ΓhRK,u−RCuK≤∑K∈ΓhαK2RKK21/2·∑K∈ΓhαK−2u−RCuK21/2,∑E⊂Ω\∂ΩDRE,u−RCuE≤∑E⊂Ω\∂ΩDε−1/2αEREE21/2·∑E⊂Ω\∂ΩDε1/2αE−1u−RCuE21/2.Further, using Lemma2, we get
(36)∑K∈ΓhRK,u−RCuK≲∑K∈ΓhαK2RKK21/2Cu,∑E⊂Ω\∂ΩDRE,u−RCuE≲∑E⊂Ω\∂ΩDε−1/2αEREE21/2Cu.Therefore, the termBv−vh,u−RCu is bounded above by
(37)Bv−vh,u−RCu≲∑K∈ΓhαK2RKK2+∑E⊂Ω\∂ΩDε−1/2αEREE21/2Cu.From (18) and (19), we have
(38)∇uhK≲minhminK−1b0−1/2,ε−1/2uhK=hminK−1αKuhKforuh∈V0h.Next, we will find the bounds on the second termBv−vh,RCu of equation (32).
(39)Bv−vh,RCu=ε∇v,∇RCu+a·∇v,RCu+bv,RCu−ε∇vh,∇RCu+a·∇vh,RCu+bvh,RCu+∑KρKRK,a·∇RCu+∑KσKRK,ah·∇RCu.Using the standard Galerkin orthogonality condition and standard scaling results, the above equation reduces to(40)Bv−vh,RCu=−∑KρKRK,a·∇RCu−∑KσKRK,ah·∇RCu≤∑K∈ΓhρKRKKa∞,K∇RCuK+∑K∈ΓhσKRKKah∞,K∇RCuK≤∑K∈ΓhρKRKKa∞,KhminK−1αKRCuK+∑K∈ΓhσKRKKah∞,KhminK−1αKRCuK.We know that for Clément operator [11]
(41)RCu≲Cu∀u∈H01Ω.Thus, we have(42)Bv−vh,RCu≤∑K∈ΓhρK2RKK2a∞,K2hminK−2αK21/2Cu+∑K∈ΓhσK2RKK2ah∞,K2hminK−2αK21/2Cu.It may be noted that the effect of nonlinear term in theL∞ norm will be bounded by that of the term a∞ as shown below, i.e.,
(43)ah=a·∇uh∇uh∇uh2,∇uh≠0,ah≤a∇uh∇uh∇uh2UsingCauchy–Schwarzinequality,ah≤a∇uh∇uh∇uh2,≤a.Sincea, the convection coefficient, is assumed to be smooth in the domain under consideration, it is bounded above. Hence, the nonlinear term ah∞,K is taken as bounded above by some constant and is absorbed in the constant term.We know that(44)ρK≲hminK/a∞,K∀K∈Γh,σK≲hminK/ah∞,K∀K∈Γh.Therefore,(45)Bv−vh,RCu≲∑K∈ΓhαK2RKK2+∑K∈ΓhαK2RKK21/2·Cu·M1u,Γh≈C.Since(46)v−vh≤Bv−vh,uu,using equation (45) and equation (37) in equation (32), we get
(47)v−vh≲∑K∈Γh3αK2RKK2+∑E⊂Ω\∂ΩDε−1/2αEREE21/2C.Using triangle inequalities,(48)RKK2≤rK−RKK2+rKK2,REE2≤rE−REE2+rEE2,we get
(49)v−vh≲C∑K∈Γh3αK2rK−RKK2+rKK2+∑E⊂∂K\∂ΩDε−1/2αErE−REE2+rEE21/2.


## 5. Adaptive Refinement Algorithm

In this section, we propose an adaptive refinement strategy based on the a posteriori error estimates obtained in the last section. We propose the following adaptive refinement algorithm:(1)
Discretize the domain using triangular elements. Triangulation is being carried out using red refinement(2)
Solve the problem using the proposed scheme described in Section2(3)
Over each elementK, the residual error estimates ηK have been calculated as defined in Section 4(4)
Mark the elementsKeiei=1M satisfying ηKei>CmaxL′ηL′, where C is a user chosen constant from 0,1, for refinement(5)
Refine these marked elements using green refinement procedure(6)
Refine all those elements having hanging nodes also to avoid the discontinuity in the solution(7)
Solve the problem again on the new adapted mesh(8)
Repeat the process of grid refinement until the solution has been obtained up to a given desired accuracy


## 6. Numerical Results

In this section, numerical experiments have been carried out in order to test the efficiency and robustness of the proposed adaptive refinement technique based on the derived error estimates.Example 4.
Consider the following singularly perturbed convection-diffusion problem:(50)−∇·ε∇u+2ux+3uy+u=finΩ=0,12,u=0on∂Ω.
The right-hand side functionf is so chosen to satisfy the exact solution
(51)u=sinx1−e−21−x/εy21−e−31−y/ε.The solution of the above problem exhibits exponential boundary layers along the linesx=1 and y=1. For adaptive refinement, an anisotropic triangular mesh has been taken into consideration. In Figures 3(a) and 3(b), we present a portion of adaptive triangular mesh for ε=2−3 with different degrees of freedom. Figures 4(a) and 4(b) present adaptive refined meshes for ε=2−5 with different degrees of freedoms.Figure 3
Portion of adaptive mesh forε=2−3 with DOF=187,224.Figure 4
Portion of adaptive mesh forε=2−5 with DOF=559,587.In Figures5 and 6, the numerical solution obtained using the proposed refinement algorithm for different values of the singular perturbation parameter ε=2−3 and ε=2−5 has been plotted. It can be easily seen that even very sharp boundary layers have been efficiently captured using the proposed refinement algorithm. From the solutions, it can also be observed that the problem is very sensitive to the singular perturbation parameter, i.e., even for ε=2−5, very sharp boundary layers appear in the solution. In Figure 7, energy norm errors for ε=2−5 have been presented. The behavior of effectivity index ψ=v−vh/ErrorEstimator which is used to measure reliability of the estimator is shown in Figure 8.Figure 5
Numerical solution profile forε=2−3.Figure 6
Numerical solution profile forε=2−5.Figure 7
Errorv−vh and error estimator for ε=2−5.Figure 8
ψ (effectivity index) for ε=2−5.


## 7. Conclusion

In the presented work, an adaptive numerical technique has been proposed for singularly perturbed convection-diffusion problems in two dimensions. The singularly perturbed problem under consideration has been solved using Hughes stabilization technique under SUPG finite element framework. Anisotropic meshes have been considered for the domain discretization. Reliable a posteriori error estimates have been developed in energy norm on anisotropic meshes. Based on these a posteriori error estimates, an adaptive mesh refinement strategy has been proposed. It has been depicted through numerical experiments that the proposed adaptive refinement strategy is very much efficient in capturing sharp boundary layers as the singular perturbation parameterε approaches to 0.



---

*Source: 1361498-2020-05-16.xml*