## 1 Introduction

In this paper we consider the following Stokes equations: Find the velocity and the pressure such that

(1.1) |

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 [29]. 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 [2].

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 [17] and Wang [46] 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 [34], the nonconforming Crouzeix-Raviart element [35, 5], the conforming element [36], the hybrid discontinuous Galerkin method [20], low and high order Taylor-Hood and mini elements [33], hybrid high-order (HHO) method [42] and weak Galerkin (WG) method [39].

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 [38], HHO method [19, 18], WG methods [47, 48], generalized barycentric coordinates methods [44], 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 [56]. 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 [49], 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.

### 2.1 Preliminaries

Letting , the weak formulation for (1.1) reads as follows: Find such that

(2.1) |

where

The weak formulation (2.1) is well posed thanks to the coercivity of the bilinear form and the inf–sup stability of the bilinear , see [25].

We introduce an auxiliary unknown . Then, the model problem (1.1) can be recast into the following first order system of equations

(2.2a) | ||||

(2.2b) | ||||

(2.2c) | ||||

(2.2d) |

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 [7]. Let be a polygon with vertices that are arranged counterclockwise, and let

be the outward unit normal vector on edge

that 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 asThen the Wachspress coordinates are defined as (cf. [45])

where , .

As discusses in [23], 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

###### Lemma 2.1.

Let

be the nodal value interpolation operator based on the generalized barycentric coordinates such as Wachspress coordinates (cf.

[45]). For any , it holds (cf. [7])(2.3) |

### 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 Figure

3)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

(2.4) |

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:

with norm

Now let us define the conforming interpolation operator on polygonal mesh in the spirit of [7]. For any , we define restricted to by

where

Note that any function is continuous over the edge , which by definition yields a function that belongs to .

Then following [56], the discrete formulation for the Stokes equations (2.2) reads as follows: Find such that

(2.5a) | ||||

(2.5b) | ||||

(2.5c) |

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

(2.6) |

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

(2.7) |

Finally, the following inf-sup conditions hold (cf. [12, 30]):

(2.8) | ||||

(2.9) |

###### Remark 2.1.

(implementation).

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 [56].

To facilitate later analysis, we define three interpolation operators , and , which are given explicitly as

It follows immediately that

(2.10) | ||||

(2.11) |

In addition, standard interpolation error estimates (cf. [16, 12]) imply

(2.12) |

Next we will present some properties which play an important role for later analysis. To this end, we decompose as (cf. [1])

(2.13) |

where is the curl of a function in whose tangent trace vanishes on . is such that and .

###### Lemma 2.2.

For and , it holds

(2.14) |

Thus, we have velocity invariance property, i.e.,

(2.15) |

In addition, is divergence free.

###### Proof.

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

This gives

which implies in . On the other hand, is continuous over , thus , which implies that is divergence free.

∎

###### Remark 2.2.

(comparison to existing methods).

We compare our method with weak Galerkin (WG) method proposed in [39]. First of all, our method shares the same degrees of freedom for pressure with WG proposed in [39], 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.

###### Lemma 2.3.

(approximation properties). We have for

###### Proof.

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

(2.17) |

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 [56]. 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.

###### Lemma 3.1.

(uniform a priori bound on the discrete velocity). The following estimate holds

(3.1) |

###### Proof.

Invoking (2.13) and (2.15), we can rewrite (2.5) as

(3.2) | ||||

(3.3) | ||||

(3.4) |

Taking in (3.2), (3.3) and (3.4), and sum up, we can obtain

(3.5) |

We can infer from (2.6), (2.8) and (3.3) that

thereby, we can infer from Lemma 2.3 and (3.5)

which yields the desired estimate by dividing both sides by .

∎

###### Remark 3.1.

Contrary to the classical estimate that can be obtained from [30] and [56], 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 [15], where smallness assumption for is required.

###### Theorem 3.1.

###### Proof.

First, we define such that

(3.6) |

Invoking (2.5b) and (3.6), we can obtain

(3.7) |

We can infer from (2.10) and (3.6) that

which means is the -orthogonal projection of onto . Thus,

An appeal to (2.6), (2.8) and (3.7) yields

(3.8) |

On the other hand, we can infer from the first equation of (1.1) and (2.16) that

Integration by parts yields

Thereby, we can obtain

(3.9) |

On the other hand, we have from the discrete formulation (2.5a)

which can be combined with (3.9) yielding

Hence

(3.10) |

Taking in (3.10), then it leads to

Since , we have

(3.11) |

It follows from (2.6), (2.11), (3.7) and (3.11) that

therefore, we have from Lemma 2.3 and (3.8)

The triangle inequality yields

Thus, we can conclude that

Comments

There are no comments yet.