A Weak Galerkin Finite Element Scheme for solving the stationary Stokes Equations
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 gradientdivergence operators. The WG method is highly flexible by allowing the use of discontinuous functions on arbitrary polygons or polyhedra with certain shape regularity. Optimalorder 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
\hb@xt@.01(1.1)  
\hb@xt@.01(1.2)  
u  \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 tensorvalued 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
where is the unit outward normal vector on the domain boundary .
This problem mainly arises from approximations of
lowReynoldsnumber flows. The finite element methods for Stokes and
NavierStokes problems enforce the
divergencefree property in
finite element spaces, which satisfy the infsup (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: , , ,
and
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
\hb@xt@.01(1.4)  
\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
It is easy to know that the weak form (LABEL:VF1)(LABEL:VF2) is also equivalent to the following variational problem: find such that
\hb@xt@.01(1.6)  
\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 hpadaptivity. 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 highcontrast parameter dependent combination of the Darcy and Stokes models. The numerical method of [17] is based on the traditional gradientdivergence variational form for the Brinkman equations. In [24], we presented a new WG scheme based on the gradientgradient 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 BrinkmanStokes 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 gradientgradient 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. Optimalorder 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
with notations
The Sobolev norm is defined as
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 vectorvalued function and the scalarvalued function, also we will introduce the weak gradients and the corresponding discrete forms.
3.1 Weak gradient for weak vectorvalued function
Let be a polygonal or polyhedral domain with boundary . A weak vectorvalued function on the domain is defined as such that and . Let
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
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
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
\hb@xt@.01(3.1) 
3.2 Weak gradient for weak scalarvalued function
We define a weak scalarvalued function on the domain as such that and . Let
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
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
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
\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
It should be noticed that is single valued on each edge . For pressure variable, we define
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,
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:
Lemma 4.1
([19]) For any and the following equations hold true
\hb@xt@.01(4.1)  
\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
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
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.
\hb@xt@.01(4.3)  
\hb@xt@.01(4.4) 
Now we introduce four bilinear forms as follows:
\hb@xt@.01(4.5)  
\hb@xt@.01(4.6)  
\hb@xt@.01(4.7)  
\hb@xt@.01(4.8) 
Using these bilinear forms we define the following two norms. For any and ,
\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
\hb@xt@.01(4.11)  
\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
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 ,
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 ,
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
\hb@xt@.01(5.1)  
\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
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
\hb@xt@.01(5.3)  
\hb@xt@.01(5.4) 
Proof. First, from (LABEL:ifw1) and the property (LABEL:qcfv) we obtain
Summing over all elements , we have
\hb@xt@.01(5.5) 
From the commutative property (LABEL:qcfp) we arrive at
\hb@xt@.01(5.6) 
It follows from (LABEL:ifw2) that
Summing over all yields
\hb@xt@.01(5.7)  
Next, using in to test (LABEL:ose1), we have
Integrating by parts, we obtain
\hb@xt@.01(5.8) 
where we have used the fact that . Using in to test (LABEL:ose2), we arrive at
Using the fact that one has
\hb@xt@.01(5.9)  
Finally combining the equations (LABEL:eefa1) and (LABEL:eefb1) with (LABEL:eefa2) yields
\hb@xt@.01(5.10)  
Substituting it into (LABEL:wf1), then we would have
Combining the equations (LABEL:eefb2) with (LABEL:eefc1) we arrive at
\hb@xt@.01(5.11)  
Substituting (LABEL:eefl1) into (LABEL:wf2) yields the following error equation
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)  
\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
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
\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
\hb@xt@.01(6.5) 
where C is a constant.
From (LABEL:qcfv), we have
It follows from the definition of , (LABEL:czi1), and (LABEL:ti) that