# Numerical Study on Several Stabilized Finite Element Methods for the Steady Incompressible Flow Problem with Damping

Authors: Jilian Wu; Pengzhan Huang; Xinlong Feng
Journal: Journal of Applied Mathematics
 (2013)

---

## Abstract

We discuss several stabilized finite element methods, which are penalty, regular, multiscale enrichment, and local Gauss integration method, for the steady incompressible flow problem with damping based on the lowest equal-order finite element space pair. Then we give the numerical comparisons between them in three numerical examples which show that the local Gauss integration method has good stability, efficiency, and accuracy properties and it is better than the others for the steady incompressible flow problem with damping on the whole. However, to our surprise, the regular method spends less CPU-time and has better accuracy properties by using Crout solver.


---

## Body

## 1. Introduction

In this paper, we will consider the following steady incompressible flow problem with damping: seek(u,p) such that
(1)-νΔu+α|u|r-2u+∇p=finΩ,∇·u=0inΩ,u=0on∂Ω,
where Ω⊂ℝ2 is a convex polygonal domain with a Lipschitz continuous boundary ∂Ω and the symbols Δ, ∇, and ∇· denote the Laplacian, gradient, and divergence operators, respectively; u=(u1(x,y),u2(x,y)), p=p(x,y), and f=(f1(x,y),f2(x,y)) represent the velocity vector, the pressure, and the prescribed body force, respectively. Further, α|u|r-2u represents the damping term with two constants α>0 and r>2. In addition, |u|=u12+u22, and we linearize the nonlinear term α|u|r-2u by allowing it to lag one step behind.These equations describe various physical situations such as porous media flow, drag or friction effects, and some dissipative mechanisms from the resistance to the motion of the flow. If the damping system is different, energy dissipation is different. So the application of these equations is very extensive in the daily life (see [1–6] and references therein).It is well known that it is very difficult to compute some PDEs directly while numerical method plays an important role in these problems, so developing an efficient and effective computational method for solving the incompressible flow problem has practical significance and has drawn the attention of many researchers (see [4, 7–12] and the references cited therein). During this time, mixed finite element methods [13] are a natural choice for solving fluid mechanics equations because these equations naturally appear in mixed form in terms of velocity and pressure.In the analysis and practice of employing mixed finite element methods in solving the incompressible flow problems, the inf-sup condition has played an important role because it ensures stability and accuracy of the underlying numerical schemes. Pairs of finite element spaces that are used to approximate the velocity and the pressure unknown are said to be stable if they satisfy the inf-sup condition. Intuitively speaking, the inf-sup condition is something that enforces a certain correlation between two finite element spaces so that they both have the required properties when employed for the incompressible flow problems. However, due to computational convenience and efficiency in practice, some mixed finite element pairs which do not satisfy the inf-sup condition are also popular. Thus, much attention has been paid to the study of the stabilized method for the Stokes problem [4, 7].In the present years, studies have focused on stabilization techniques [14, 15], which include penalty method [16, 17], regular method [15], multiscale enrichment method [18], and local Gauss integration method [19]. There exist a lot of theoretical results for the stabilized mixed finite element methods for the Stokes equations and the comparisons between them are also given (see [4, 7, 15, 16, 18–21] and the references cited therein). It is pointed out that authors considered the performance of several stabilized methods for the Stokes equations based on the lowest equal-order pairs in [21]. And authors studied several stabilized methods for the Stokes eigenvalue problem by using different conforming and nonconforming lower-order pairs in [7]. In this report, we will adopt four kinds of stabilized finite element methods but will mainly focus on the incompressible flow problem with damping based on the lowest equal-order pairs. Moreover, we present the comparisons between these methods of the considered problem.A brief outline of the rest of our paper is organized as follows: in Section2, we introduce some notations and present some preliminary materials and some well-known results of the steady incompressible flow problem with damping to be used in our subsequent sections; then in Section 3, we review several stabilized mixed finite element methods and recall their key stabilization techniques; in Section 4, comparisons between these stabilized methods are performed numerically; finally, we end with some short conclusions in Section 5.


## 2. Notations and Preliminaries

We will use the following Hilbert spaces:(2)X=H01(Ω)2∩Lr(Ω)2,Y=L2(Ω)2,M=L02(Ω)={q∈L2(Ω):∫Ω‍qdx=0}.
Here and in what follows, the spaces L2(Ω)m (m=1,2) are equipped with the L2-scalar product (·,·) and L2-norm ∥·∥L2 or ∥·∥0, respectively. Further, we will consider the standard definitions for Sobolev spaces Wm,p(Ω) which, for any integer m>0 and any number p≥1, are equipped with the norm
(3)∥u∥m,p,Ω=∥u∥m,p.
Notice that
(4)Wm,2(Ω)=Hm(Ω),∥·∥m,2=∥·∥m.
Then the variational formulation of problem (1) is to seek (u,p)∈X×M such that
(5)B((u,p);(v,q))+α(|u|r-2u,v)=(f,v),ffff∀(v,q)∈X×M,
where
(6)B((u,p);(v,q))=a(u,v)+d(v,p)+d(u,q),g1g∀(u,p),(v,q)∈X×M,a(u,v)=ν(∇u,∇v),∀u,v∈X,d(v,p)=-(p,divv),∀(v,p)∈X×M.Now, for convenience, we introduce the finite subspacesXh×Mh⊂X×M, assumed to be uniformly regular in the usual sense. Suppose that Kh is a triangular decomposition of the domain Ω and h is the maximum mesh size of the partition. Therefore, we define
(7)Xh={u∈C0(Ω-)2∩X:u|K∈P1(K)2,∀K∈Kh},Mh={q∈C0(Ω-)∩M:q|K∈P1(K),∀K∈Kh},
where P1(K) represents the space of linear functions on K. And we assume the following basis functions:
(8)Xh=Xh×Xh,Xh=span{ϕ1h,…,ϕnh},Mh=span{φ1h,…,φmh},
where n and m are the dimensions of Xh and Mh, respectively. And the bilinear forms a(vh,vh) and d(vh,qh) satisfy the following conditions [7].(i)
There is a constantc1>0 such that
(9)a(vh,vh)≥c1∥vh∥1,Ω2,∀vh∈Zh,where Zh={vh∈Xh∣d(vh,qh)=0,∀qh∈Mh}.(ii)
There is a constantc2>0 independent of h such that
(10)supvh∈Xh|d(vh,qh)|∥vh∥0≥c2∥qh∥,∀(vh,qh)∈Xh×Mh.Then we have a unique solution (uh,ph) of (5) satisfying
(11)∥u-uh∥1,Ω+∥p-ph∥0,Ω≤c3(infvh∈Xh∥u-vh∥1,Ω+infqh∈Mh∥p-qh∥0,Ω),
where c3 is a positive constant.


## 3. Stabilized Mixed Finite Element Methods

In this section, we will give several stabilized mixed finite element algorithms to show their different aspects and several ways have been used to stabilize the lowest equal-order finite element space pair as follows (see [7, 21]). First we introduce the following classical Uzawa iterative algorithm.LetH1 and H2 be finite dimensional spaces. We consider the following saddle point problems:
(12)[ABTB-C][XY]=[fg],
where X∈H1 and Y∈H2 are unknown variables; f∈H1, g∈H1; A:H1→H1 is symmetric positive definite operator, B:H1→H2 is linear map, BT:H2→H1 is the transpose operator of B, and C:H2→H2 is symmetric semidefinite operator. If the initial values X0∈H1 and Y0∈H2 are given, then Xi and Yi(i=1,2,…) are defined by
(13)Xi=Xi-1+A-1(f-AXi-1-BTYi-1),Yi=Yi-1+τ(BXi-CYi-1-g),
where τ is a given real number.LeteiY=Y-Yi be the iteration error generated by the above method. It is easy to show that
(14)eiY=(I-τ(BA-1BT+C))ei-1Y.
Let λ1 denote the largest eigenvalue of matrix BA-1BT+C. Then Yi converges to Y if τ is chosen such that 0<τ≤2/λ1. In this case, Xi and Yi converge, respectively, to X and Y with a rate of convergence bounded by the absolute value of 1-τλ1. For more details about the saddle point problems, please see [22–24] and the references cited therein.Next we present four kinds of stabilized finite element method for the steady incompressible flow problem with damping.Remark 1 (nonconforming finite element space).
For convenience, we leth be a positive parameter and Kh={Kj} a regular triangulation of Ω. Denote by Γj=∂Ω∩∂Kj the boundary edge and by Γjk=Γkj=∂Kj∩∂Kk the interior boundary. Set the centers of Γj and Γjk by ζj and ζjk, respectively. The nonconforming finite element space can be defined as
(15)NCh={v:vj=v|Kj∈P1(Kj)2,vj(ζjk)=vk(ζkj),fv(ζj)=0,Kj∈Kh,∀j,kP1(Kj)2},
where P1(Kj) is the set of all polynomials on Kj of degree less than 1. Note that NCh is not a subspace of X. However, in this nonconforming case, the pair of finite element spaces is NCh×Mh; that is, the conforming space is still used for pressure.So the corresponding discrete variational formulation of (5) for the Stokes equations with damping reads as follows. Seek (uh,ph)∈Xh×Mh such that
(16)B((uh,ph);(v,q))+α(|uh|r-2uh,v)=(f,v),gggg∀(v,q)∈Xh×Mh.
Then we can get the following equations from (16):
(17)A1{u1,h}+αA2{u1,h}+BxT{ph}={∫Ω‍ϕif1},A1{u2,h}+αA2{u2,h}+ByT{ph}={∫Ω‍ϕif2},Bx{u1,h}+By{u2,h}=0.
Next, let Uh and Ph be the array of the velocity and pressure, respectively. Then it is easy to see that (17) can be written in matrix form:
(18)[A1+αA2B1TB10][UhPh]=[Fh0],
where  A1=diag(A1,A1), A1=(A1i,j), and A1i,j=ν(∇ϕjh,∇ϕih), i,j=1,…,n; A2=diag(A2,A2), A2=(A2i,j), and A2i,j=𝒰(ϕjh,ϕih), i,j=1,2,…,n; B1=(Bx,By), Bx=(Bxij), By=(Byij), Bxij=-∫Ω‍(∂ϕjh/∂x)φih, and Byij=-∫Ω‍(∂ϕjh/∂y)φih, i,j=1,…,n; Uh=[{u1,h},{u2,h}]T; Fh=[{∫Ω‍ϕihf1},{∫Ω‍ϕihf2}]T; and Ph=[p1,…,pm]T.Remark 2.
Here,𝒰=({u1,h}2+{u2,h}2)r-2. Because we linearized the nonlinear term α|u|r-2u, then 𝒰 is a known-term in the process of solving (18). And for (16), the process of linearizing is as follows.
Given(uhn,phn), seek (uhn+1,phn+1)∈Xh×Mh such that
(19)B((uhn+1,phn+1);(v,q))+α(|uhn|r-2uhn+1,v)=(f,v).
In the nonconforming case, the discrete nonconforming formulation for the steady incompressible flow problem with damping is to seek(uh,ph)∈NCh×Mh such that
(20)B′((uh,ph);(v,q))+α∑K∈Kh‍(|uh|r-2uh,v)K=∑K∈Kh‍(f,v)K,
where
(21)B′((u,p);(v,q))=a′(u,v)+d′(v,p)+d′(u,q),a′(u,v)=ν∑K∈Kh‍(∇u,∇v)K,d′(v,p)=-∑K∈Kh‍(divv,p)K.
Note that (17) is a saddle point problem and the lowest equal-order pair does not satisfy the discrete inf-sup condition
(22)supvh∈Xhd(vh,ph)∥∇vh∥0≥c4∥ph∥0orsupvh∈NChd′(vh,ph)∥∇vh∥0,h≥c4∥ph∥0∀qh∈Mh,
where the constant c4>0 is independent of h and ∥∇vh∥0,h=(∑j‍|v|1,kj2)2, for all v∈NCh.Algorithm  I (Penalty method). The penalty method compensates for the inf-sup condition deficiency by adding the penalty term as follows.Seek(uh,ph)∈Xh×Mh such that
(23)B((uh,ph);(v,q))+α(|uh|r-2uh,v)+εν(ph,q)=(f,v),gg∀(v,q)∈Xh×Mh,
where ε>0 is a penalty parameter. The performance of this method obviously depends on the choice of the penalty parameter ε. Then the matrix form of (23) can be expressed as
(24)[A1+αA2B1TB1ενD1][UhPh]=[Fh0],
where the matrixes A1 and B1 are presented and D1 is deduced from (ph,q), using the base for Mh in the usual manner; that is,
(25)D1=(D1i,j),D1i,j=∫Ω‍φihφjhdxdy,i,j=1,2,…,m.
Let eiPh=Ph-Pih(i=1,2,…); we use the above Uzawa iterative algorithm and get
(26)ei+1Ph=(I-τ(B1(A1+αA2)-1B1T-ενD1))eiPh.Algorithm  II (Regular method). This method uses a simple way to stabilize the mixed finite element approximation without a loss of accuracy, that is, to seek (uh,ph)∈Xh×Mh such that
(27)B((uh,ph);(v,q))+α(|uh|r-2uh,v)-δ∑K∈Kh‍(∇ph-f,∇q)K=(f,v),
for all (v,q)∈Xh×Mh, where δ=h2/(βν) is a stabilization parameter and β>0. The matrix form of the above stabilized version can be expressed as
(28)[A1+αA2B1TB1δD2][UhPh]=[FhC1h],
where additional blocks D2 and C1h correspond to
(29)-∑K∈Kh‍(∇ph,∇q)K,-δ∑K∈Kh‍(f,∇q)K,
respectively; that is,
(30)D2=(D2i,j),D2i,j=-∑K∈Kh‍∫K‍∇φih∇φjhdxdy,ggggfgggi,j=1,…,m,
and C1h=[-δ{∑‍∫K‍(∂φih/∂x)f1},-δ{∑‍∫K‍(∂φih/∂y)f2}]T, where
(31){∑‍∫K‍∂φih∂xf1}=∑‍[∫K‍∂φ1h∂xf1,…,∫K‍∂φmh∂xf1]T,{∑‍∫K‍∂φih∂yf2}=∑‍[∫K‍∂φ1h∂yf2,…,∫K‍∂φmh∂yf2]T.
Similar to (26), we also have
(32)ei+1Ph=(I-τ(B1(A1+αA2)-1B1T-δD2))eiPh.Algorithm  III (Multiscale enrichment method). Another stabilized way is the multiscale enrichment approach which includes the usual Galerkin least squares stabilized terms on each finite element and positive jump terms at interelement boundaries. Namely, seek (uh,ph)∈Xh×Mh such that
(33)B((uh,ph);(v,q))+α(|uh|r-2uh,v)-δ1∑K∈Kh‍(∇ph-f,∇q)K+δ2∑Γjk‍([ν∂nuh],[ν∂nv])Γjk=(f,v),fffffffffffffffffff∀(v,q)∈Xh×Mh,
where δ1=h2/(β1ν) and δ2=h/(β2ν) are the positive stabilization parameters, n is the normal outward vector, ∂n is normal derivative operator, and [v] denotes the jump of v across e. Moreover, a direct algebraic manipulation leads to the matrix form
(34)[A1+αA2B1TB1+δ2D3δ1D2][UhPh]=[FhC2h],
where the matrix D3 is deduced from the term δ2∑Γjk‍([ν∂nuh],[ν∂nv])Γjk.Similar to (26), we have
(35)ei+1Ph=(I-τ((B1+δ2D3)(A1+αA2)-1B1T-δ1D2))eiPh.Algorithm  IV (Local Gauss integration method). The local Gauss integration method is to add two Gauss integrals rather than any stabilization parameter to the original discrete formulation (16) as follows. Seek (uh,ph)∈Xh×Mh such that
(36)B((uh,ph);(v,q))+α(|uh|r-2uh,v)-G(ph,q)=(f,v),(v,q)∈Xh×Mh,
where G(ph,q) is defined by
(37)G(ph,q)=δ3∑K∈Kh‍{∫K,2‍phqdx-∫K,1‍phqdx},∀ph,q∈Mh,δ3>0, and ∫K,i‍g(x)dx indicates a local Gauss integral over K that is exact for polynomials of degree i(i=1,2). Then the corresponding matrix form of this stabilized method is
(38)[A1+αA2B1TB1δ3G][UhPh]=[Fh0],
where
(39)G=-∑K∈Kh‍(G2-G1),G2=(∫K,2‍φihφjh),G1=(∫K,1‍φihφjh),i,j=1,…,m.
As (26), we get
(40)ei+1Ph=(I-τ(B1(A1+αA2)-1B1T-δ3G))eiPh.It is well known that, ifτ is well chosen, then Ui and Pi converge, respectively, to U and P with a rate of convergence based on (26), (32), (35), and (40). From these equations, we can find that the penalty method converges faster than the local Gauss integration method. We can obtain that coefficient matrices of penalty algorithm, regular algorithm, and local Gauss integration algorithm are all symmetric; however, the multiscale enrichment method’s coefficient matrices are not symmetric from (10), (28), (34), and (38). What is more, it is easy to see that the matrix calculations of multiscale enrichment algorithm are more complex and cumbersome, so this method maybe costs more time.By using the regularity assumptions and well-established techniques for velocity and pressure [4, 7], the theoretical convergence rates should be of order O(h) for the velocity in the L2-norm and of order O(h2) for the pressure in the H1-norm, respectively, by using all these stabilized methods.


## 4. Numerical Experiments

In this section, we will give three numerical tests to confirm the numerical theory developed in the previous section. In the given experiments, the pressure and velocity are approximated by the lowest equal-order finite element pairs defined with respect to the same uniform triangulation; that is, the mesh consists of triangular elements that are obtained by dividingΩ into subsquares of equal size and then drawing the diagonal in each subsquare.


### 4.1. Numerical Test 1

In this example, we consider the exact solution problem firstly. Let the domainΩ be the unit square Ω=(0,1)×(0,1)⊂ℝ2. The exact solution for the velocity u=(u1,u2) and pressure p is given as follows:
(41)p(x,y)=cos(πx)cos(πy),u1(x,y)=2πsin2(πx)sin(πy)cos(πy),u2(x,y)=-2πsin(πx)cos(πx)sin2(πy),
and the right-hand side f=(f1(x,y),f2(x,y)) is determined by the original problem (1). Our goal in this test is to compare CPU-time, the L2-error of the pressure, and H1-error of the velocity; the experimental rates of convergence for these methods with different values of h are tabulated in Tables 1, 2, 3, 4, and 5. What is more, the rates of convergence are calculated by the formula log(Ei/Ei+1)/log(hi/hi+1), where Ei and Ei+1 are the relative errors corresponding to the meshes of sizes hi and hi+1, respectively.Table 1
Numerical results for the penalty method withν=1.0e-4, α=1.0e-4, and r=3.
1
/
h
CPU-time
∥
u
-
u
h
∥
1
/
∥
u
∥
1
∥
p
-
p
h
∥
0
/
∥
p
∥
0
u
H
1-rate
p
L
2-rate
12
0.375
2.3745
E
-
1
1.3649
E
-
1
—
—
24
3.265
1.1800
E
-
1
6.8136
E
-
2
1.0088
1.0023
36
10.39
7.8528
E
-
2
4.5442
E
-
2
1.0045
0.9990
48
27.562
5.8844
E
-
2
3.4113
E
-
2
1.0031
0.9968
60
61.25
4.7051
E
-
2
2.7326
E
-
2
1.0023
0.9942Table 2
Results got from the regular method withν=1.0e-4, α=1.0e-4, and r=3.
1
/
h
CPU-time
∥
u
-
u
h
∥
1
/
∥
u
∥
1
∥
p
-
p
h
∥
0
/
∥
p
∥
0
u
H
1-rate
p
L
2-rate
12
0.282
1.1945
E
-
0
1.5179
E
-
2
—
—
24
1.641
4.4196
E
-
1
4.7747
E
-
3
1.4344
1.6686
36
4.657
2.4423
E
-
1
2.5361
E
-
3
1.4628
1.5605
48
10.313
1.6070
E
-
1
1.6313
E
-
3
1.4550
1.5337
60
21.094
1.1646
E
-
1
1.1619
E
-
3
1.4431
1.5207Table 3
Results got from the multiscale enrichment method withν=1.0e-4, α=1.0e-4, and r=3.
1
/
h
CPU-time
∥
u
-
u
h
∥
1
/
∥
u
∥
1
∥
p
-
p
h
∥
0
/
∥
p
∥
0
u
H
1-rate
p
L
2-rate
12
2.496
0.8924
1.1544
—
—
24
12.746
0.6523
0.2823
0.4520
2.0319
36
37.175
0.4853
0.2094
0.7293
0.7370
48
83.648
0.3737
0.1713
0.9091
0.6977
60
161.32
0.2996
0.1445
0.9900
0.7624Table 4
Results got from the local Gauss integration method withν=1.0e-4, α=1.0e-4, and r=3.
1
/
h
CPU-time
∥
u
-
u
h
∥
1
/
∥
u
∥
1
∥
p
-
p
h
∥
0
/
∥
p
∥
0
u
H
1-rate
p
L
2-rate
12
0.406
2.2348
E
-
1
8.9420
E
-
3
—
—
24
4.578
1.1458
E
-
1
2.8981
E
-
3
1.0353
1.6255
36
8.5
7.4890
E
-
2
1.8494
E
-
3
1.0488
1.1078
48
16.328
5.5377
E
-
2
1.3620
E
-
3
1.0493
1.0635
60
27.39
4.3894
E
-
2
1.0468
E
-
3
1.0414
1.1797Table 5
Results got from the local Gauss integration method with the nonconforming element withν=1.0e-4, α=1.0e-4, and r=3.
1
/
h
CPU-time
∥
u
-
u
h
∥
1
/
∥
u
∥
1
∥
p
-
p
h
∥
0
/
∥
p
∥
0
u
H
1-rate
p
L
2-rate
12
0.297
3.6007
5.8335
E
-
2
—
—
24
2.500
1.8851
1.6524
E
-
2
0.9336
1.8198
36
5.735
1.2682
7.5810
E
-
3
0.9776
1.9217
48
11.469
0.9543
4.3219
E
-
3
0.9887
1.9534
60
21.719
0.7646
2.7858
E
-
3
0.9932
1.9681From this experiment, we can learn the following several points.(1) For penalty method, the result of this method is well in which parameter “r” can be 2, 3, and 4 and we can use UMFPACK or default solver in the process of calculation. (2) Results got from the penalty, regular, multiscale enrichment, and local Gauss integration methods by using conforming and nonconforming elements are presented in Tables 1–5, respectively. Here, we choose ε=1.0e-6, α=1.0e-4, r=3, Re=10000, and β=160, β1=160, β2=100, because they can deal with the considered problem well. (3) For regular method, it works well only we choose Crout solver; however, for other methods, several solvers can be used and their difference is not big. What is more, the value of “r” has little influence on the results within a certain range; in this paper, we give the results of r=3 presented in tables.From Tables1–5, we can see that these methods work well and keep the convergence rates just like the theoretical analysis except for the multiscale enrichment method. Meanwhile, it can be seen that the penalty method requires the least CPU-time, which validates the analysis in Section 3. As expected, we have an interesting observation that the error of the nonconforming local Gauss integration method is better than that of the conforming version, which is not surprising since the degree of freedom of the nonconforming method is nearly three times than that of the conforming one on uniform mesh. Hence, it is natural that the nonconforming local Gauss integration method is more accurate and costs more CPU-time. The penalty method should use less time than the regular method and local Gauss integration method theoretically but in fact not; maybe this is caused by the different solver.Besides, from the convergence results on this example, we can see that regular and multiscale enrichment methods are not better than other methods. And the nonconforming local Gauss integration method shows the best numerical stability.


### 4.2. Numerical Test 2

In this test, we test a popular benchmark problem, the lid-driven flow. Let the computation be carried out in the regionΩ={(x,y)∣0<x,y<1}. We assume the normal component of the velocity to be zero on ∂Ω and the tangential component to be zero except along y=1, where it is set to one.In this example, we simulate the referred physics phenomena. In Figures1, 2, 3, 4, and 5, we present the velocity streamlines and pressure level lines for ε=10e-6, α=100, h=1/30, and r=4 based on these five methods. From Figures 1(b)–5(b), only Gauss methods can obtain resolved pressure. For the velocity, from Figures 1(a)–5(a), we can see that Gauss methods can capture this model better than the other methods.Velocity streamlines (a) and pressure level lines (b) for the penalty method withε=1.0e-6, α=100, and r=4, respectively.
(a)
(b)Velocity streamlines (a) and pressure level lines (b) for the regular method withε=1.0e-6, α=100, and r=4, respectively.
(a)
(b)Velocity streamlines (a) and pressure level lines (b) for the multiscale enrichment method withε=1.0e-6, α=100, and r=4, respectively.
(a)
(b)Velocity streamlines (a) and pressure level lines (b) for the local Gauss integration method withε=1.0e-6, α=100, and r=4, respectively.
(a)
(b)Velocity streamlines (a) and pressure level lines (b) for the local Gauss integration method with the nonconforming element withε=1.0e-6, α=100, and r=4, respectively.
(a)
(b)


### 4.3. Numerical Test 3

In this example, we choose the exact solution for the velocity and pressure in the unit square as follows:(42)p(x,y)=10(2x-1)(2y-1),u1(x,y)=10x2(x-1)2y(y-1)(2y-1),u2(x,y)=-10x(x-1)(2x-1)y2(y-1)2,
and the right-hand side is determined just as Test 1. In this example, parameters we choose are ν=1.0e-4, α=1.0e-4, and r=3. For the penalty and local Gauss integration methods, CPU-time, the L2-error of the pressure, H1-error of the velocity, and the experimental rates of convergence for these methods with different values of h are tabulated in Tables 6 and 7. What is more, the numerical solutions are given in Figure 6. From this test, we can get conclusion that is similar to Test 1. However, the result got from multiscale method is not better.Table 6
Numerical results for the penalty method withν=1.0e-4, α=1.0e-4, and r=3.
1
/
h
CPU-time
∥
u
-
u
h
∥
1
/
∥
u
∥
1
∥
p
-
p
h
∥
0
/
∥
p
∥
0
u
H
1-rate
p
L
2-rate
12
0.219
2.8476
E
-
1
1.8917
E
-
1
—
—
24
1.625
1.4126
E
-
1
9.8071
E
-
2
1.0114
0.9478
36
5.719
9.3881
E
-
2
6.6229
E
-
2
1.0077
0.9682
48
15.64
7.0298
E
-
2
5.0001
E
-
2
1.0056
0.9771
60
32.281
5.6183
E
-
2
4.0161
E
-
2
1.0044
0.9820Table 7
Results got from the local Gauss integration method withν=1.0e-4, α=1.0e-4, and r=3.
1
/
h
CPU-time
∥
u
-
u
h
∥
1
/
∥
u
∥
1
∥
p
-
p
h
∥
0
/
∥
p
∥
0
u
H
1-rate
p
L
2-rate
12
0.219
1.78153
E
-
1
2.01165
E
-
2
—
—
24
3.109
7.64659
E
-
1
4.97672
E
-
3
1.22023
2.01511
36
5.437
4.78738
E
-
1
2.20800
E
-
3
1.15491
2.00432
48
9.578
3.40798
E
-
1
1.24253
E
-
3
1.18138
1.99853
60
16.062
2.58872
E
-
1
7.96681
E
-
4
1.23219
1.99177Numerical solutions for the penalty method (a) and pressure level lines for the local Gauss integration method (b) withε=1.0e-4, α=1.0e-4, and r=4, respectively.
(a)
(b)


## 4.1. Numerical Test 1

In this example, we consider the exact solution problem firstly. Let the domainΩ be the unit square Ω=(0,1)×(0,1)⊂ℝ2. The exact solution for the velocity u=(u1,u2) and pressure p is given as follows:
(41)p(x,y)=cos(πx)cos(πy),u1(x,y)=2πsin2(πx)sin(πy)cos(πy),u2(x,y)=-2πsin(πx)cos(πx)sin2(πy),
and the right-hand side f=(f1(x,y),f2(x,y)) is determined by the original problem (1). Our goal in this test is to compare CPU-time, the L2-error of the pressure, and H1-error of the velocity; the experimental rates of convergence for these methods with different values of h are tabulated in Tables 1, 2, 3, 4, and 5. What is more, the rates of convergence are calculated by the formula log(Ei/Ei+1)/log(hi/hi+1), where Ei and Ei+1 are the relative errors corresponding to the meshes of sizes hi and hi+1, respectively.Table 1
Numerical results for the penalty method withν=1.0e-4, α=1.0e-4, and r=3.
1
/
h
CPU-time
∥
u
-
u
h
∥
1
/
∥
u
∥
1
∥
p
-
p
h
∥
0
/
∥
p
∥
0
u
H
1-rate
p
L
2-rate
12
0.375
2.3745
E
-
1
1.3649
E
-
1
—
—
24
3.265
1.1800
E
-
1
6.8136
E
-
2
1.0088
1.0023
36
10.39
7.8528
E
-
2
4.5442
E
-
2
1.0045
0.9990
48
27.562
5.8844
E
-
2
3.4113
E
-
2
1.0031
0.9968
60
61.25
4.7051
E
-
2
2.7326
E
-
2
1.0023
0.9942Table 2
Results got from the regular method withν=1.0e-4, α=1.0e-4, and r=3.
1
/
h
CPU-time
∥
u
-
u
h
∥
1
/
∥
u
∥
1
∥
p
-
p
h
∥
0
/
∥
p
∥
0
u
H
1-rate
p
L
2-rate
12
0.282
1.1945
E
-
0
1.5179
E
-
2
—
—
24
1.641
4.4196
E
-
1
4.7747
E
-
3
1.4344
1.6686
36
4.657
2.4423
E
-
1
2.5361
E
-
3
1.4628
1.5605
48
10.313
1.6070
E
-
1
1.6313
E
-
3
1.4550
1.5337
60
21.094
1.1646
E
-
1
1.1619
E
-
3
1.4431
1.5207Table 3
Results got from the multiscale enrichment method withν=1.0e-4, α=1.0e-4, and r=3.
1
/
h
CPU-time
∥
u
-
u
h
∥
1
/
∥
u
∥
1
∥
p
-
p
h
∥
0
/
∥
p
∥
0
u
H
1-rate
p
L
2-rate
12
2.496
0.8924
1.1544
—
—
24
12.746
0.6523
0.2823
0.4520
2.0319
36
37.175
0.4853
0.2094
0.7293
0.7370
48
83.648
0.3737
0.1713
0.9091
0.6977
60
161.32
0.2996
0.1445
0.9900
0.7624Table 4
Results got from the local Gauss integration method withν=1.0e-4, α=1.0e-4, and r=3.
1
/
h
CPU-time
∥
u
-
u
h
∥
1
/
∥
u
∥
1
∥
p
-
p
h
∥
0
/
∥
p
∥
0
u
H
1-rate
p
L
2-rate
12
0.406
2.2348
E
-
1
8.9420
E
-
3
—
—
24
4.578
1.1458
E
-
1
2.8981
E
-
3
1.0353
1.6255
36
8.5
7.4890
E
-
2
1.8494
E
-
3
1.0488
1.1078
48
16.328
5.5377
E
-
2
1.3620
E
-
3
1.0493
1.0635
60
27.39
4.3894
E
-
2
1.0468
E
-
3
1.0414
1.1797Table 5
Results got from the local Gauss integration method with the nonconforming element withν=1.0e-4, α=1.0e-4, and r=3.
1
/
h
CPU-time
∥
u
-
u
h
∥
1
/
∥
u
∥
1
∥
p
-
p
h
∥
0
/
∥
p
∥
0
u
H
1-rate
p
L
2-rate
12
0.297
3.6007
5.8335
E
-
2
—
—
24
2.500
1.8851
1.6524
E
-
2
0.9336
1.8198
36
5.735
1.2682
7.5810
E
-
3
0.9776
1.9217
48
11.469
0.9543
4.3219
E
-
3
0.9887
1.9534
60
21.719
0.7646
2.7858
E
-
3
0.9932
1.9681From this experiment, we can learn the following several points.(1) For penalty method, the result of this method is well in which parameter “r” can be 2, 3, and 4 and we can use UMFPACK or default solver in the process of calculation. (2) Results got from the penalty, regular, multiscale enrichment, and local Gauss integration methods by using conforming and nonconforming elements are presented in Tables 1–5, respectively. Here, we choose ε=1.0e-6, α=1.0e-4, r=3, Re=10000, and β=160, β1=160, β2=100, because they can deal with the considered problem well. (3) For regular method, it works well only we choose Crout solver; however, for other methods, several solvers can be used and their difference is not big. What is more, the value of “r” has little influence on the results within a certain range; in this paper, we give the results of r=3 presented in tables.From Tables1–5, we can see that these methods work well and keep the convergence rates just like the theoretical analysis except for the multiscale enrichment method. Meanwhile, it can be seen that the penalty method requires the least CPU-time, which validates the analysis in Section 3. As expected, we have an interesting observation that the error of the nonconforming local Gauss integration method is better than that of the conforming version, which is not surprising since the degree of freedom of the nonconforming method is nearly three times than that of the conforming one on uniform mesh. Hence, it is natural that the nonconforming local Gauss integration method is more accurate and costs more CPU-time. The penalty method should use less time than the regular method and local Gauss integration method theoretically but in fact not; maybe this is caused by the different solver.Besides, from the convergence results on this example, we can see that regular and multiscale enrichment methods are not better than other methods. And the nonconforming local Gauss integration method shows the best numerical stability.


## 4.2. Numerical Test 2

In this test, we test a popular benchmark problem, the lid-driven flow. Let the computation be carried out in the regionΩ={(x,y)∣0<x,y<1}. We assume the normal component of the velocity to be zero on ∂Ω and the tangential component to be zero except along y=1, where it is set to one.In this example, we simulate the referred physics phenomena. In Figures1, 2, 3, 4, and 5, we present the velocity streamlines and pressure level lines for ε=10e-6, α=100, h=1/30, and r=4 based on these five methods. From Figures 1(b)–5(b), only Gauss methods can obtain resolved pressure. For the velocity, from Figures 1(a)–5(a), we can see that Gauss methods can capture this model better than the other methods.Velocity streamlines (a) and pressure level lines (b) for the penalty method withε=1.0e-6, α=100, and r=4, respectively.
(a)
(b)Velocity streamlines (a) and pressure level lines (b) for the regular method withε=1.0e-6, α=100, and r=4, respectively.
(a)
(b)Velocity streamlines (a) and pressure level lines (b) for the multiscale enrichment method withε=1.0e-6, α=100, and r=4, respectively.
(a)
(b)Velocity streamlines (a) and pressure level lines (b) for the local Gauss integration method withε=1.0e-6, α=100, and r=4, respectively.
(a)
(b)Velocity streamlines (a) and pressure level lines (b) for the local Gauss integration method with the nonconforming element withε=1.0e-6, α=100, and r=4, respectively.
(a)
(b)


## 4.3. Numerical Test 3

In this example, we choose the exact solution for the velocity and pressure in the unit square as follows:(42)p(x,y)=10(2x-1)(2y-1),u1(x,y)=10x2(x-1)2y(y-1)(2y-1),u2(x,y)=-10x(x-1)(2x-1)y2(y-1)2,
and the right-hand side is determined just as Test 1. In this example, parameters we choose are ν=1.0e-4, α=1.0e-4, and r=3. For the penalty and local Gauss integration methods, CPU-time, the L2-error of the pressure, H1-error of the velocity, and the experimental rates of convergence for these methods with different values of h are tabulated in Tables 6 and 7. What is more, the numerical solutions are given in Figure 6. From this test, we can get conclusion that is similar to Test 1. However, the result got from multiscale method is not better.Table 6
Numerical results for the penalty method withν=1.0e-4, α=1.0e-4, and r=3.
1
/
h
CPU-time
∥
u
-
u
h
∥
1
/
∥
u
∥
1
∥
p
-
p
h
∥
0
/
∥
p
∥
0
u
H
1-rate
p
L
2-rate
12
0.219
2.8476
E
-
1
1.8917
E
-
1
—
—
24
1.625
1.4126
E
-
1
9.8071
E
-
2
1.0114
0.9478
36
5.719
9.3881
E
-
2
6.6229
E
-
2
1.0077
0.9682
48
15.64
7.0298
E
-
2
5.0001
E
-
2
1.0056
0.9771
60
32.281
5.6183
E
-
2
4.0161
E
-
2
1.0044
0.9820Table 7
Results got from the local Gauss integration method withν=1.0e-4, α=1.0e-4, and r=3.
1
/
h
CPU-time
∥
u
-
u
h
∥
1
/
∥
u
∥
1
∥
p
-
p
h
∥
0
/
∥
p
∥
0
u
H
1-rate
p
L
2-rate
12
0.219
1.78153
E
-
1
2.01165
E
-
2
—
—
24
3.109
7.64659
E
-
1
4.97672
E
-
3
1.22023
2.01511
36
5.437
4.78738
E
-
1
2.20800
E
-
3
1.15491
2.00432
48
9.578
3.40798
E
-
1
1.24253
E
-
3
1.18138
1.99853
60
16.062
2.58872
E
-
1
7.96681
E
-
4
1.23219
1.99177Numerical solutions for the penalty method (a) and pressure level lines for the local Gauss integration method (b) withε=1.0e-4, α=1.0e-4, and r=4, respectively.
(a)
(b)


## 5. Conclusions

We have used several stabilized mixed finite element methods in solving the steady incompressible flow problem with damping based on the lowest equal-order pairs in this paper. We give some conclusions by comparing numerically as follows.All of those methods’ stability and efficiency depends on their parameter values. In terms of the penalty method, the smaller its parameter value, the more stable the method. However, it can not be too small, otherwise, the condition number of the system matrix arising from this method will become too large to be solve. For the regular and multiscale enrichment methods whose performance heavily depends on the choice of the stabilization parameters, however, it is difficult to choose fine parameters in fact. What is more, a poor choice of these stabilization parameters can also lead to serious deterioration in the convergence rates. The local Gauss integration method is free of stabilization parameters and shows numerically the best performance among the methods considered for the given problem.



---

*Source: 985864-2013-11-26.xml*