In this paper we consider the following Stokes equations: Find the velocity and the pressure such that
where is the computational domain in , is a constant viscosity parameter and is the body force.
A large amount of work has been dedicated to solving (1.1), see [25, 32] and the references therein. However, classical finite element methods can not deliver exactly divergence free solutions in the sense of and the divergence constraint is relaxed in order to guarantee inf-sup stability . As such, these methods suffer from a lack of pressure robustness, specifically, their velocity error estimates are pressure dependent and possibly deteriorates unboundedly for , which shows some kind of locking phenomena in the sense of .
To remedy this issue, various approaches have been proposed. Exactly divergence-free or conforming methods of order is developed in [43, 50, 22]. An alternative strategy to achieve divergence-free property is to enrich the -conforming elements with divergence free rational functions, which provides the correct flexibility to enforce (strong) continuity (cf. [26, 27]). Moreover, Cockburn  and Wang  modify the variational formulation and introduce the tangential penalty and thus obtain discontinuous Galerkin divergence-free schemes. Another approach that can deliver pressure independent error estimates is to add grad-div stabilization [41, 40]. Recently, a method based on modifying the right hand side is shown to be able to retrieve pressure robustness, and this strategy has been successfully applied to finite volume method , the nonconforming Crouzeix-Raviart element [35, 5], the conforming element , the hybrid discontinuous Galerkin method , low and high order Taylor-Hood and mini elements , hybrid high-order (HHO) method  and weak Galerkin (WG) method .
Recently, polygonal finite element methods have become a hot topic due to their great flexibility in handling complicated geometries. Various approaches in the framework of polygonal grids have been proposed to solve partial differential equations arising from practical applications, such as virtual element method (VEM)[3, 4], mimetic finite difference methods , HHO method [19, 18], WG methods [47, 48], generalized barycentric coordinates methods , staggered discontinuous Galerkin (SDG) method [52, 56], etc. SDG as the new generation discretisation methods for PDEs based on discrete unknowns that enjoys staggered continuity properties is initially introduced to solve wave propagation problem [12, 13]. Its close connection to a hybridizable DG method is unveiled in [10, 11]. Note that SDG method differs from other DG methods in the sense that the basis functions are locally conforming and penalty terms are not required. This method has many desirable features and has been successfully applied to a wide range of partial differential equations [9, 21, 52, 56, 53, 54, 51, 55, 31].
The purpose of this paper is to develop a pressure robust staggered discontinuous Galerkin (PR-SDG) method on general polygonal meshes for the Stokes equations with piecewise constant approximations. Our approach is based on the framework of lowest order SDG method introduced in . The novelty here is to exploit velocity reconstruction operator in the discretization of the source term in the spirit of [34, 39]. The principal idea behind this is that discrete divergence-free velocity test functions are mapped to exact divergence-free ones by velocity reconstruction operator. The crux of SDG method for the Stokes equations is to generate a sub-triangulation by connecting an interior point of the polygonal grid to all the vertices of the polygon. Then finite element spaces that enjoy staggered continuity properties for velocity gradient, velocity and pressure are developed, which naturally lead to inf-sup stable pairs. In particular, the velocity space is continuous over the edges of the polygon, which enables us to establish velocity reconstruction operator based on the polygonal grid. A rigorous convergence analysis for errors of velocity gradient, velocity and pressure is investigated. The main difficulty lies in the proof of superconvergence and traditional techniques no long work. To attack this issue, we exploit Aubin-Nitsche duality argument with nonstandard incorporation of continuous and discrete formulation of the dual problem; in addition, velocity reconstruction operator is employed in the discretization of the dual problem. We keep track explicitly of the dependence on the viscosity in the analysis so as to address the practically important issue of body forces with large irrotational part. Note that the lowest order SDG method can be viewed as finite volume method, thus our work provides new perspectives for the understanding of finite volume method. Indeed, our primal and dual partitions are exactly the same as the finite volume method proposed in , whereas our method is based on the first order system and piecewise constant functions are exploited for all the involved unknowns. Importantly, the continuity of all the unknowns are staggered on the interelement boundaries, which naturally yields a stable numerical scheme without resorting to further stabilization. We emphasize that the development of pressure robust method with velocity reconstruction operator on polygonal mesh is still in its infancy [8, 39], the approach developed in this paper will definitely inspire more works in this direction.
The rest of the paper is organized as follows. In the next section, we formulate the PR-SDG method for the Stokes equations based on velocity reconstruction operator. Then in section 3, we prove the optimal convergence for errors of velocity gradient, velocity and pressure as well as superconvergence of velocity. Several numerical experiments are carried out in section 4 to verify the proposed theories. Finally, a conclusion is given.
2 Description of PR-SDG method
In this section we begin with introducing the construction of minimal degree conforming function on convex polygon. Then we propose the PR-SDG method by employing conforming velocity reconstruction in the discretization of the body force. Finally, some fundamental properties inherited from the proposed method are provided.
Letting , the weak formulation for (1.1) reads as follows: Find such that
We introduce an auxiliary unknown . Then, the model problem (1.1) can be recast into the following first order system of equations
Before closing this subsection we introduce some notations that will be employed throughout the paper. Let , we adopt the standard notations for the Sobolev spaces and their associated norms , and semi-norms for . The space coincides with , for which the norm is denoted as . We use to denote the inner product for . If , the subscript will be dropped unless otherwise mentioned. In the sequel, we use to denote a generic positive constant which may have different values at different occurrences.
2.2 conforming function on polygon
In this subsection, we briefly introduce the construction of conforming function on convex polygon by following . Let be a polygon with vertices that are arranged counterclockwise, and let
be the outward unit normal vector on edgethat connects vertices and , where we conveniently denote when the subscript is not in the range of , see Figure 1 for an illustration. Let be an interior point of , its distance to edge and a scaled normal vector are defined as
Then the Wachspress coordinates are defined as (cf. )
where , .
As discusses in , one introduces an auxiliary ratio function
By the quotient rule for differentiation, one is led to
Then, a counterclockwise rotation of the gradient is defined by
Now we are ready to present the construction of minimal degree conforming finite element function on convex polygon. Denote as the area and the length of edge for . Let be an arbitrary point inside polygon , and denote by the area of the triangle formed by . Then for each , we define by
where , and . Here is the Kronecker symbol.
Furthermore, for these basis functions, we have the following properties
2.3 Construction of PR-SDG method
In this subsection we will present the construction of PR-SDG method. To begin with, we will describe the construction of the staggered mesh needed in our method by following [12, 13, 14, 30]. We let be the initial partition of the domain into non-overlapping simple convex polygonal elements (primal mesh). We use to represent the set of all the edges in this partition (primal edges) and to represent the subset of all interior edges, that is, . For each polygon in the initial partition , we select an interior point and create new edges by connecting to the vertices of polygon. This process will divide into the union of triangles, where the triangle is denoted as . Moreover, we will use to denote the set of all the new edges generated by this subdivision process (dual edges) and use to denote the resulting triangulation (simplicial submeshes), on which our basis functions are defined. In addition, we define and . This construction is illustrated in Figure 2, where solid lines are edges in and dashed lines are edges in . For each triangle , we let be the diameter of and . Here we assume that the primal partition satisfies the standard mesh regularity assumptions: 1) Every primal element is star-shaped with respect to a ball of radius , where is a positive constant. 2) For every primal element and every edge , it satisfies , where is a positive constant. Note that these assumptions can guarantee that the resulting triangulation is shape regular.
For each interior edge , we use to denote the dual mesh, which is the union of the two triangles in sharing the edge , and for each boundary edge , we use to denote the triangle in having the edge , see Figure 2. We define a unit normal vector on each edge as follows: If is a boundary edge, then we define as the unit normal vector of pointing towards outside of . If is an interior edge, then we fix as one of the two possible unit normal vectors on . We will use instead of to simplify the notation when there is no confusion.
Let be the order of approximation. For every and , we define and as the spaces of polynomials of degree less than or equal to on and , respectively. Now we are ready to describe the finite element spaces that will be used to define our numerical scheme. First, the locally conforming SDG space for velocity is defined as:
The degrees of freedom for this space can be described as (see Figure3)
The discrete norm and norm for the space are given by
Here and denotes the jump on an interior edge that is shared by two triangles and belonging to , and we simply take for .
Next, the locally conforming SDG space for the velocity gradient approximation is defined as:
which is equipped by
Invoking scaling arguments, we have
We define the following degrees of freedom for and it is illustrated in Figure 3
In the above definition, the jump over an edge is defined as
where , is the common edge of the two triangles and that belong to , and is a unit normal to the edge .
Finally, locally conforming finite element space for pressure is defined as:
Now let us define the conforming interpolation operator on polygonal mesh in the spirit of . For any , we define restricted to by
Note that any function is continuous over the edge , which by definition yields a function that belongs to .
where the bilinear forms and are defined as
and the bilinear forms and are defined as
Here, we list some important properties that will be employed later. First, integration by parts yields the following adjoint property
Next, we notice that for , , we have from integration by parts and the definition of
thereby summing up over all the elements yields the following adjoint property
Since the bilinear forms and are adjoint, we only need to compute in the actual implementation. This bypasses the computation of the global velocity reconstruction in the bilinear form, thus our implementation only modifies the right hand side assembling compared to .
To facilitate later analysis, we define three interpolation operators , and , which are given explicitly as
It follows immediately that
Next we will present some properties which play an important role for later analysis. To this end, we decompose as (cf. )
where is the curl of a function in whose tangent trace vanishes on . is such that and .
For and , it holds
Thus, we have velocity invariance property, i.e.,
In addition, is divergence free.
For , we have from integration by parts
which gives (2.14). Thereby, we can obtain
Finally, we will show that is divergence free. For , we have from integration by parts and the definition of
In addition, we can infer from (2.5c) that
which implies in . On the other hand, is continuous over , thus , which implies that is divergence free.
(comparison to existing methods).
We compare our method with weak Galerkin (WG) method proposed in . First of all, our method shares the same degrees of freedom for pressure with WG proposed in , see Figures 3 and 4. Our method has less degrees of freedom for velocity since our velocity space only consists of edge degrees of freedom while WG consists of edge and interior degrees of freedom (cf. Figure 4). On the other hand, our approach is based on the first order system, thus velocity gradient can be calculated simultaneously. However, WG is based on the primal formulation, straightforward calculation of velocity gradient is not available.
(approximation properties). We have for
For any , notice that is a constant over the sub-triangle of and we use to denote the value of restricted to each sub-triangle . Let and we use to represent restricted to , then we have
where we use the fact that if the right side vanishes, then the left side also vanishes. An application of scaling arguments yields the desired estimate.
Since for , we can proceed analogously to (2.17) to obtain
where Lemma 2.1 is employed in the last inequality.
Finally, the triangle inequality yields
Summing up over all the primal elements yields the desired estimate.
3 A priori error analysis
In this section we will present the convergence estimates for all the variables involved. In particular we will prove that the error estimates for velocity are independent of . The main contribution of this section is to prove the superconvergence of velocity via duality argument, which is non-trivial.
The unique solvability of the discrete formulation (2.5) can be proved in line with . The proof will not be repeated here for the sake of conciseness. We will give the uniform a priori bound for velocity in the next lemma.
(uniform a priori bound on the discrete velocity). The following estimate holds
which yields the desired estimate by dividing both sides by .
Contrary to the classical estimate that can be obtained from  and , the a priori bound (3.1) persists in the limit . This bound can be incorporated into Navier-Stokes equations to establish error estimates under smallness assumption that only concerns the solenoidal part of the body force. This can improve the existing bound given in , where smallness assumption for is required.
First, we define such that
which means is the -orthogonal projection of onto . Thus,
Integration by parts yields
Thereby, we can obtain
On the other hand, we have from the discrete formulation (2.5a)
which can be combined with (3.9) yielding
Taking in (3.10), then it leads to
Since , we have
The triangle inequality yields
Thus, we can conclude that