A Weak Galerkin Finite Element Scheme for solving the stationary Stokes Equations

# A Weak Galerkin Finite Element Scheme for solving the stationary Stokes Equations

Ruishu Wang Department of Mathematics, Jilin University, Changchun, China (ruishu@email.jlu.edu.cn).    Xiaoshen Wang Department of Mathematics, University of Arkansas at Little Rock, Little Rock, AR 72204, United States(xxwang@ualr.edu).    Qilong Zhai Department of Mathematics, Jilin University, Changchun, China (diql13@mails.jlu.edu.cn).    Ran Zhang Department of Mathematics, Jilin University, Changchun, China (zhangran@mail.jlu.edu.cn). The research of Zhang was supported in part by China Natural National Science Foundation(11271157, 11371171, 11471141), and by the Program for New Century Excellent Talents in University of Ministry of Education of China.
###### Abstract

A weak Galerkin (WG) finite element method for solving the stationary Stokes equations in two- or three- dimensional spaces by using discontinuous piecewise polynomials is developed and analyzed. The variational form we considered is based on two gradient operators which is different from the usual gradient-divergence operators. The WG method is highly flexible by allowing the use of discontinuous functions on arbitrary polygons or polyhedra with certain shape regularity. Optimal-order error estimates are established for the corresponding WG finite element solutions in various norms. Numerical results are presented to illustrate the theoretical analysis of the new WG finite element scheme for Stokes problems.

Key words. weak Galerkin finite element methods, weak gradient, Stokes equations, polytopal meshes.

AMS subject classifications. Primary, 65N30, 65N15, 65N12, 74N20; Secondary, 35B45, 35J50, 35J35

## 1 Introduction

The aim of this paper is to present a novel weak Galerkin finite element method for solving the stationary Stokes equations. Let be a polygonal or polyhedral domain in . As a model for the flow of an incompressible viscous fluid confined in , we consider the following equations

 −μΔu+∇p = f,in Ω, \hb@xt@.01(1.1) ∇⋅u = 0,in Ω, \hb@xt@.01(1.2) u = g,on ∂Ω, \hb@xt@.01(1.3)

for unknown velocity function and pressure function (we require that has zero average in order to guarantee the uniqueness of the pressure). Bold symbols are used to denote vector- or tensor-valued functions or spaces of such functions. Here is a body source term, is the kinematic viscosity and is a boundary condition that satisfies the compatibility condition

 ∫∂Ωg⋅n ds=0,

where is the unit outward normal vector on the domain boundary .

This problem mainly arises from approximations of low-Reynolds-number flows. The finite element methods for Stokes and NavierStokes problems enforce the
divergence-free property in finite element spaces, which satisfy the inf-sup (LBB) condition, in order for them to be numerically stable [2, 1, 10, 11, 8]. The Stokes problem has been studied with various different new numerical methods: [4, 12, 13, 22, 23].

Throughout this paper, we would follow the standard definitions for Lebesgue and Sobolev spaces: , , ,

 [H10(Ω)]d={v∈[H1(Ω)]d:v=0 on ∂Ω}

and

 L20(Ω):={q∈L2(Ω):∫Ωqdx=0}

are the natural spaces for the weak form of the Stokes problem [10, 7]. Denote for inner products in the corresponding spaces.

Next we assume that and . Then one of the variational formulations for the Stokes problem (LABEL:ose1)-(LABEL:ose3) is to find and such that

 (∇u,∇v)−(∇⋅v,p) = (f,v), \hb@xt@.01(1.4) (∇⋅u,q) = 0, \hb@xt@.01(1.5)

for all and . Here denotes the velocity gradient tensor . It is well known that under our assumptions on the domain and the data, problem (LABEL:VF1)-(LABEL:VF2) has a unique solution .

For any , define a functional such that

 ⟨∇p,v⟩=−(∇⋅v,p),∀v∈[H10(Ω)]d.

It is easy to know that the weak form (LABEL:VF1)-(LABEL:VF2) is also equivalent to the following variational problem: find such that

 (∇u,∇v)+⟨∇p,v⟩ = (f,v), \hb@xt@.01(1.6) ⟨∇q,u⟩ = 0, \hb@xt@.01(1.7)

for all and . The unique solvability of (LABEL:VF3)-(LABEL:VF4) follows directly from that of the (LABEL:VF1)-(LABEL:VF2).

The WG method refers to a general finite element technique for partial differential equations where differential operators are approximated as distributions for generalized functions. This method was first proposed in [20, 21, 15] for second order elliptic problem, then extended to other partial differential equations [14, 16, 18, 17, 25, 26]. Weak functions and weak derivatives can be approximated by polynomials with various degrees. The WG method uses weak functions and their weak derivatives which are defined as distributions. The most prominent features of it are:

• The usual derivatives are replaced by distributions or discrete approximations of distributions.

• The approximating functions are discontinuous. The flexibility of discontinuous functions gives WG methods many advantages, such as high order of accuracy, high parallelizability, localizability, and easy handling of complicated geometries.

The above features motivate the use of WG methods for the Stokes equations. It can easily handle meshes with hanging nodes, elements of general shapes with certain shape regularity and ideally suited for hp-adaptivity. In [19], Wang et. al. considered WG methods for the Stokes equations (LABEL:VF1)-(LABEL:VF2). Similarly, in [17], they presented WG methods for the Brinkman equations, which is a model with a high-contrast parameter dependent combination of the Darcy and Stokes models. The numerical method of [17] is based on the traditional gradient-divergence variational form for the Brinkman equations. In [24], we presented a new WG scheme based on the gradient-gradient variational form. It is shown that this scheme is suit for the mixed formulation of Darcy which would present a better approximation for this case. In fact, for complex porous media with interface conditions, people often use Brinkman-Stokes interface model to describe this problem, which is an ongoing work for us now. In order to present a more efficient WG scheme, we prefer to utilize this gradient-gradient weak form to approximate the model. In order to unify the weak form of this interface problem, we need the numerical analysis results of this form for Stokes problem. However, to the best of our knowledge, the numerical analysis of methods based on the variational form (LABEL:VF3)-(LABEL:VF4) has never been done before. Therefore in this paper, we propose a WG method based on the weak form (LABEL:VF3)-(LABEL:VF4) of the primary problem. In addition, if we choose high order polynomials to approximate the model and use Schur complement to reduce the interior DOF of the velocity and pressure by the boundary DOF, the total DOF of this new method is less than the scheme of [19].

The rest of this paper is organized as follows. In Section 2 we shall introduce some preliminaries and notations for Sobolev spaces. Section 3 is devoted to the definitions of weak functions and weak derivatives. The WG finite element schemes for variational form of the Stokes equation (LABEL:VF3)-(LABEL:VF4) are presented in Section 4. This section also contains some local projection operators and then derives some approximation properties which are useful in a convergence analysis. In Section 5, we derive an error equation for the WG finite element approximation. Optimal-order error estimates for the WG finite element approximations are derived in Section 6 in an -equivalent norm for the velocity, and norm for both the velocity and the pressure. In Section 7, we present some numerical results which confirm the theory developed in earlier sections. Finally, we present some technical estimates in the appendix for quantities related to the local projections into various finite element spaces.

## 2 Preliminaries and Notations

Let be an open bounded domain with Lipschitz continuous boundary in . We shall use standard definitions of the Sobolev spaces and inner products , their norms , and seminorms , for any . For instance, for any integer , the seminorm is defined as

 |v|s,K=⎛⎝∑|α|=s∫K|∂αv|2dK⎞⎠12,

with notations

 α=(α1,α2,⋯,αd),|α|=α1+α2+⋯+αd,∂α=d∏j=1∂αjxj.

The Sobolev norm is defined as

 ∥v∥m,K=(m∑j=0|v|2j,K)12.

The space is same as , whose norm and inner product are denoted by and , respectively. If , we would drop the subscript in the notations of the norm and the inner product.

## 3 Weak Differential Operators

In this section we will define weak functions for both the vector-valued function and the scalar-valued function, also we will introduce the weak gradients and the corresponding discrete forms.

### 3.1 Weak gradient for weak vector-valued function

Let be a polygonal or polyhedral domain with boundary . A weak vector-valued function on the domain is defined as such that and . Let

 V(T)={v={v0,vb};v0∈[L2(T)]d,vb∈[L2(∂T)]d},

where is not necessarily the trace of .

###### Definition 3.1

([19]) For any , the weak gradient of , denoted by , is defined as a linear functional in the dual space of whose action on each is given by

 (∇wv,τ)T=−(v0,∇⋅τ)T+⟨vb,τ⋅n⟩∂T,   ∀τ∈[H1(T)]d×d,

where is the outer unit normal vector to , is the inner product of and , and is the inner product of and in .

Consider the inclusion map defined below

 iV(ϕ)={ϕ|T,ϕ|∂T},    ϕ∈[H1(T)]d.

By this map the Sobolev space can be embedded into the space . With the help of map , the Sobolev space can be considered as a subspace of by identifying each with .

Let be the set of polynomials on T with degree no more than .

###### Definition 3.2

([19]) The discrete weak gradient operator is defined as follows: for each , is the unique element such that

 (∇w,r,Tv,τ)T=−(v0,∇⋅τ)T+⟨vb,τ⋅n⟩∂T,    ∀τ∈[Pr(T)]d×d. \hb@xt@.01(3.1)

### 3.2 Weak gradient for weak scalar-valued function

We define a weak scalar-valued function on the domain as such that and . Let

 W(T)={q={q0,qb};q0∈L2(T),qb∈L2(∂T)},

where is not necessarily the trace of .

###### Definition 3.3

([20]) For any , the weak gradient of , denote by , is defined as a linear functional in the dual space of whose action on each is given by

 (˜∇wq,w)T=−(q0,∇⋅w)T+⟨qb,w⋅n⟩∂T,   ∀w∈[H1(T)]2,

where is the outer unit normal vector to , is the inner product of and , and is the inner product of and in .

Consider the inclusion map defined as follows

 iW(ϕ)={ϕ|T,ϕ|∂T},    ϕ∈H1(T).

By which the Sobolev space is embedded into the space . With the help of map , the Sobolev space can be considered as a subspace of by identifying each with .

###### Definition 3.4

([20]) The discrete weak gradient operator is defined as follows: for each , is the unique element such that

 (˜∇w,r,Tq,w)T=−(q0,∇⋅w)T+⟨qb,w⋅n⟩∂T,    ∀w∈[Pr(T)]d. \hb@xt@.01(3.2)

## 4 A Weak Galerkin Finite Element Scheme

Let be a partition of the domain into polygons in 2D or polyhedral in 3D. Assume that is shape regular in the sense as defined in [18]. Denote by the set of all edges or flat faces in , and let be the set of all interior edges or flat faces. Denote by the diameter of and the meshsize for the partition .

For any interger , we define weak Galerkin finite element spaces as follows: for velocity variable, let

 Vh={v={v0,vb};{v0,vb}|T∈[Pk(T)]d×[Pk(e)]d,e⊂∂T,vb=0 on ∂Ω}.

It should be noticed that is single valued on each edge . For pressure variable, we define

 Wh=⎧⎨⎩q={q0,qb};∑T∈Th∫Tq0dT=0,{q0,qb}|T∈Pk−1(T)×Pk(e),e⊂∂T⎫⎬⎭.

Also is single valued on each edge .

The discrete weak gradients and on the spaces and can be computed by the equations (LABEL:dwvg) and (LABEL:dwqg) on each element respectively, that is,

 (∇w,k−1v)|T = ∇w,k−1,T(v|T),∀v∈Vh, (˜∇w,kq)|T = ˜∇w,k,T(q|T),   ∀q∈Wh.

For the sake of simplicity, we shall drop the subscripts and of and in the rest of the paper.

We use the inner product to denote the sum of inner products on each of the elements as follows:

 (∇wv,∇ww) = ∑T∈Th(∇wv,∇ww)T, (˜∇wq,v) = ∑T∈Th(˜∇wq,v)T.
###### Lemma 4.1

([19]) For any and the following equations hold true

 (∇wv,τ)T = (∇v0,τ)T−⟨v0−vb,τ⋅n⟩∂T,    ∀τ∈[Pk−1(T)]d×d, \hb@xt@.01(4.1) (˜∇wp,w)T = (∇p0,w)T−⟨p0−pb,w⋅n⟩∂T,   ∀w∈[Pk(T)]d. \hb@xt@.01(4.2)

For each element , denote by the projection operator from onto . For each edge or face , denote by the projection from onto .We shall combine with as a projection onto , such that on each element

 Qhu={Q0u,Qbu}.

On each element , denote by the projection onto . Denote by the projection operator from onto . For each edge or face , denote by the projection from onto . We shall combine with as a projection onto space , such that on each element

 ˜Qhq={˜Q0q,˜Qbq}.

Then we shall present a useful property which indicates the discrete weak gradient operators are good approximation to the gradient operators in the classical sense.

###### Lemma 4.2

([19]) The following equations hold true.

 ∇wQhv = Qh∇v,   ∀v∈[H1(Ω)]d, \hb@xt@.01(4.3) ˜∇w˜Qhp = Q0∇p,    ∀p∈H1(Ω). \hb@xt@.01(4.4)

Now we introduce four bilinear forms as follows:

 s(w,v) = ∑T∈Thh−1T⟨w0−wb,v0−vb⟩∂T, \hb@xt@.01(4.5) a(w,v) = (∇ww,∇wv)+s(w,v), \hb@xt@.01(4.6) b(w,q) = (w0,˜∇wq), \hb@xt@.01(4.7) c(ρ,q) = ∑T∈ThhT⟨ρ0−ρb,q0−qb⟩∂T. \hb@xt@.01(4.8)

Using these bilinear forms we define the following two norms. For any and ,

 |||v|||2 = a(v,v)=(∇wv,∇wv)+∑T∈Thh−1T⟨v0−vb,v0−vb⟩∂T, \hb@xt@.01(4.9)

and

 = \hb@xt@.01(4.10)

where is a seminorm.

It is easy to verify that and are norms in and , respectively,

###### Weak Galerkin Algorithm 1

A numerical approximation for (LABEL:ose1)-(LABEL:ose3) can be obtained by seeking and such that

 a(uh,v)+b(v,ph) = (f,v0), \hb@xt@.01(4.11) b(uh,q)−c(ph,q) = 0, \hb@xt@.01(4.12)

for all and .

Next we shall show that the weak Galerkin finite element algorithm (LABEL:wf1)-(LABEL:wf2) has only one solution. Since the system is linear, it suffices to show that if , the only solution is .

###### Lemma 4.3

The WG finite element scheme (LABEL:wf1)-(LABEL:wf2) has a unique solution.

Proof. Let , we shall show that the solution of (LABEL:wf1)-(LABEL:wf2) is trivial. To this end, taking and and subtracting (LABEL:wf2) from (LABEL:wf1) we arrive at

 a(uh,uh)+c(ph,ph)=0.

By the definition of and , we know on each , , and on each . Thus and are continuous.

By (LABEL:ifw1) and the fact that on we have, for any ,

 0 = (∇wuh,τ)T = (∇u0,τ)T−⟨u0−ub,τ⋅n⟩∂T = (∇u0,τ)T,

which implies on each and thus is a constant. Since on each and on , we arrive at in . It follows from (LABEL:wf1), , and that for any ,

 0 = b(v,ph) = (v0,˜∇wph) = ∑T∈Th(v0,∇p0)T−∑T∈Th⟨v0⋅n,p0−pb⟩∂T = ∑T∈Th(v0,∇p0)T.

Hence we have on each . Thus is a constant in . From , we would obtain in . Since on each , .

This completes the proof of the lemma.

## 5 Error Equation

In this section, we shall derive the error equations for the WG finite element solution we get from (LABEL:wf1)-(LABEL:wf2). This error equation is essential for the following analysis.

Now we define two bilinear forms

 l1(w,v) = ∑T∈Th⟨v0−vb,(∇w−Qh∇w)⋅n⟩∂T, \hb@xt@.01(5.1) l2(w,q) = ∑T∈Th⟨q0−qb,(w−Q0w)⋅n⟩∂T, \hb@xt@.01(5.2)

for all and .

Let be the exact solution of (LABEL:ose1)-(LABEL:ose3), and be the solution of (LABEL:wf1)-(LABEL:wf2).

Define

 eh=Qhu−uh,   εh=˜Qhp−ph.

We shall derive the error equations that and satisfy.

###### Lemma 5.1

Let and be the solution of the numerical scheme (LABEL:wf1)-(LABEL:wf2), and be the exact solution of (LABEL:ose1)-(LABEL:ose3). Then, for any and we have

 a(eh,v)+b(v,εh) = s(Qhu,v)+l1(u,v), \hb@xt@.01(5.3) b(eh,q)−c(εh,q) = l2(u,q)−c(˜Qhp,q). \hb@xt@.01(5.4)

Proof. First, from (LABEL:ifw1) and the property (LABEL:qcfv) we obtain

 (∇wQhu,∇wv)T= (Qh∇u,∇wv)T = (∇v0,Qh∇u)T−⟨v0−vb,(Qh∇u)⋅n⟩∂T = (∇v0,∇u)T−⟨v0−vb,(Qh∇u)⋅n⟩∂T.

Summing over all elements , we have

 (∇wQhu,∇wv)=(∇u,∇v0)−∑T∈Th⟨v0−vb,(Qh∇u)⋅n⟩∂T. \hb@xt@.01(5.5)

From the commutative property (LABEL:qcfp) we arrive at

 b(v,˜Qhp)=(v0,˜∇w˜Qhp)=(v0,∇p). \hb@xt@.01(5.6)

It follows from (LABEL:ifw2) that

 (Q0u,˜∇wq)T= (∇q0,Q0u)T−⟨q0−qb,Q0u⋅n⟩∂T =

Summing over all yields

 b(Qhu,q) =(Q0u,˜∇wq) \hb@xt@.01(5.7) =(u,∇q0)−∑T∈Th⟨q0−qb,Q0u⋅n⟩∂T.

Next, using in to test (LABEL:ose1), we have

 (−Δu,v0)+(v0,∇p) = (f,v0).

Integrating by parts, we obtain

 (∇u,∇v0)+(v0,∇p) = (f,v0)+∑T∈Th⟨∇u⋅n,v0−vb⟩∂T, \hb@xt@.01(5.8)

where we have used the fact that . Using in to test (LABEL:ose2), we arrive at

 (∇⋅u,q0)=0.

Using the fact that one has

 0= (∇⋅u,q0) \hb@xt@.01(5.9) = −(u,∇q0)+∑T∈Th⟨u⋅n,q0⟩∂T = −(u,∇q0)+∑T∈Th⟨u⋅n,q0−qb⟩∂T.

Finally combining the equations (LABEL:eefa1) and (LABEL:eefb1) with (LABEL:eefa2) yields

 a(Qhu,v)+b(v,˜Qhp)= (∇wQhu,∇wv)+s(Qhu,v)+(v0,˜∇w˜Qhp) \hb@xt@.01(5.10) = (f,v0)+∑T∈Th⟨v0−vb,(∇u−Qh∇u)⋅n⟩∂T +s(Qhu,v).

Substituting it into (LABEL:wf1), then we would have

 a(eh,v)+b(v,εh)=s(Qhu,v)+l1(u,v).

Combining the equations (LABEL:eefb2) with (LABEL:eefc1) we arrive at

 b(Qhu,q)−c(˜Qhp,q)= (∇q0,u)−∑T∈Th⟨q0−qb,Q0u⋅n⟩∂T−c(˜Qhp,q) \hb@xt@.01(5.11) = ∑T∈Th⟨u⋅n,q0−qb⟩∂T−∑T∈Th⟨q0−qb,Q0u⋅n⟩∂T −c(˜Qhp,q) = ∑T∈Th⟨q0−qb,(u−Q0u)⋅n⟩∂T−c(˜Qhp,q).

Substituting (LABEL:eefl1) into (LABEL:wf2) yields the following error equation

 b(eh,q)−c(εh,q)=l2(u,q)−c(˜Qhp,q),

for all , which completes the proof of (LABEL:ee2).

## 6 Error Estimates

In this section we shall present the error estimates between the exact solution of (LABEL:ose1)-(LABEL:ose3) and the numerical solution of WG finite element method (LABEL:wf1)-(LABEL:wf2). The two norms and are essentially norm and norm on and respectively. In this section we always assume is shape regular ([18]).

###### Theorem 6.1

Let be the exact solution of (LABEL:ose1)-(LABEL:ose3), be the numerical solution of (LABEL:wf1)-(LABEL:wf2), then the following error estimates hold true

 \hb@xt@.01(6.1) ∥ε0∥≤Chk(∥u∥k+1+∥p∥k), \hb@xt@.01(6.2)

and consequently, one has

 \hb@xt@.01(6.3)

Proof. Letting in (LABEL:ee1) and in (LABEL:ee2), we would obtain

 s(Qhu,eh)+l1(u,eh) −l2(u,εh)+c(˜Qhp,εh).

Then from (LABEL:iqfs)-(LABEL:iqfc) we arrive at

from which we would have

For any given , it follows from [19, 6, 7, 5, 10, 11] that there is a such that

 (∇⋅˜v,ρ)∥˜v∥1≥C∥ρ∥, \hb@xt@.01(6.4)

where C is a positive constant which is dependent only on . Let , we claim that the following inequality holds true

 |||v|||≤C∥˜v∥1, \hb@xt@.01(6.5)

where C is a constant.

From (LABEL:qcfv), we have

 ∑T∈Th∥∇wv∥2T=∑T∈Th∥∇w(Qh˜v)∥2T=∑T∈Th∥Qh(∇˜v)∥2T≤∥∇˜v∥2.

It follows from the definition of , (LABEL:czi1), and (LABEL:ti) that

 ∑T∈Thh−1T∥v0−vb∥2∂T =∑T∈Thh−1T∥Q0˜v−Qb˜v∥2∂T =∑T∈Thh−1T∥Qb(Q0˜v)−Qb˜v∥2∂T ≤∑T∈Thh−1T∥Q0˜v−˜v∥2∂T ≤C∥∇˜v∥