Stabilised hybrid discontinuous Galerkin methods for the Stokes problem with non-standard boundary conditions

In several studies it has been observed that, when using stabilised $\mathbb{P}_k^{}\times\mathbb{P}_k^{}$ elements for both velocity and pressure, the error for the pressure is smaller, or even of a higher order in some cases, than the one obtained when using inf-sup stable $\mathbb{P}_k^{}\times\mathbb{P}_{k-1}^{}$ (although no formal proof of either of these facts has been given). This increase in polynomial order requires the introduction of stabilising terms, since the finite element pairs used do not stability the inf-sup condition. With this motivation, we apply the stabilisation approach to the hybrid discontinuous Galerkin discretisation for the Stokes problem with non-standard boundary conditions.


Introduction
The interest of this paper is to discretise the Stokes problem with non-standard boundary conditions. In [1], a hybrid discontinuous Galerkin (hdG) method was proposed and analysed for this problem. The finite element method used was the combination of BDM elements of order k for the velocity, and discontinuous elements of order k − 1 for the pressure. In this paper we increase the order of the pressure space to k, while keeping the order for the velocity space fixed as k. Since this pair does not satisfy the inf-sup condition, a stabilisation term needs to be added.
The stabilisation term referred to above can be built using a diversity of approaches, but, roughly speaking, the stabilisation can be residual or non-residual. In [8] the authors added a mesh-dependent term penalising the gradient of the pressure to the formulation. Later, in [14] this method was restricted and reinterpreted as a Petrov-Galerkin scheme leading to the first consistent stabilised method, and further developments were presented in the works [7] and [13]. For a review of different residual stabilised finite element methods for the Stokes problem, see the review paper [2]. Now, due to their nature, residual methods include unphysical couplings to the formulation, and modify all the entries of the stiffness matrix. Hence, non-residual methods where only a positive semi-definite term penalising the pressure is added have also being proposed. Examples of this type of methods are the pressure gradient projection [9] and local pressure gradient stabilisation [3]. The methods just mentioned typically use two nested meshes in order to build the method. Thus, to avoid this complication, the local pressure gradient stabilisation has been also presented on the same mesh in [12]. Additionally, methods that use fluctuations of the pressure gradient are not effective when the finite element space for pressure is the piecewise constant space. The usual way to overcome this is to add pressure jumps to the formulation, as it has been done, e.g., in [16]. These have been shown to be very effective, but they do somehow temper with the data structure of the code. To avoid this, the authors in [10] present an approach that is based on polynomial-pressure-projection. This method works for low order of polynomials as was shown in [4], and preserves symmetry of the original equation.
In the light of the discussion of the previous paragraphs, in this work we propose a stabilised hdG method for the Stokes problem with non-standard boundary conditions. The method is reminiscent of the Dorhmann-Bochev method (from [10]), but uses the same velocity space used in the hdG method from [1].

Notations and model problem
Let Ω be an open polygonal domain in R 2 with Lipschitz boundary Γ := ∂Ω. We use boldface font for tensor or vector variables e.g. u is a velocity vector field. The scalar variables will be italic e.g. p denotes pressure scalar value. We define the stress tensor σ := ν∇u − pI (where ν > 0 is the fluid viscosity and I is the identity matrix) and the flux as σn := σ n. In addition, we denote normal and tangential components as follows un := u · n, ut := u · t, σnn := σn · n, where n is the outward unit normal vector to the boundary Γ and t is a vector tangential to Γ such that n · t = 0. For D ⊂ Ω, we use the standard L 2 (D) space with the following norm Let us define, for m ∈ N, the following Sobolev spaces where, for α = (α1, α2) ∈ N 2 , |α| = α1 + α2, and ∂ α = . In addition, we will use the standard semi-norm and norm for the Sobolev space H m (D) In this work, we consider the two dimensional Stokes problem with tangential-velocity and normal-flux (TVNF) boundary conditions where u :Ω → R 2 is the unknown velocity field, p :Ω → R the pressure, ν > 0 the viscosity, which is considered to be constant, and f ∈ [L 2 (Ω)] 2 , g ∈ L 2 (Γ) are given functions. The restriction to homogeneous Dirichlet conditions on ut is made only to simplify the presentation. Let {T h } h>0 be a regular family of triangulations ofΩ made of triangles. For each triangulation T h , E h denotes the set of its edges. In addition, for each of element K ∈ T h , hK := diam(K), and we denote h := maxK∈T h hK . We define following Sobolev spaces on the triangulation T h and the set of all edges in E h with the corresponding broken norms. Now we will introduce the finite element spaces that discretise the above spaces. Let k ≥ 1. We start by introducing the velocity and pressure spaces. To discretise the velocity u we use the Brezzi-Douglas-Marini space (see [5, Section 2.3.1]) of order k ≥ 1 defined by Associated to this space, we introduce the BDM projection Π k : The pressure is discretised using the following space Associated to this space we define the local , and we define the continuous projection Ψ k |K = Ψ k K for all K ∈ T h . The last ingredient needed in the method described below is a finite element space associated to a family of Lagrange multipliers associated to the edges of the triangulation. These multipliers will be denoted byũ and are meant to approximate the tangential trace of the velocity u on the edges of the triangulation. For this, and in order to propose a discretisation with fewer degrees of freedom, we discretise the Lagrange multiplierũ using the space

The stabilised method
Our approach is to write the discrete problem with the same degree of polynomials for velocity and pressure spaces. In other words, denoting Vh := BDM k h × M k−1 h,0 , we want to use the space Vh × Q k h , instead of Vh × Q k−1 h as it was done in [1]. To do this, we need the proper stabilisation term, because this choice of spaces does not guarantee inf-sup stability.
The first ingredient in the definition of the stabilised method for (1.1) we use the same bilinear forms as in [1], this is where ε ∈ {−1, 1} and τ > 0 is a stabilisation parameter. In addition, to compensate for the non-inf-sup stability of the finite element spaces we have chosen, we introduce the bilinear form With these ingredients we can now present the finite element method analysed in this work:

Well-posedness of the discrete problem
Let us consider the following norm on Vh (see [1, Lemma 3.2] for a proof that this is actually a norm in Vh) The first step towards proving the stability of Method (2.2) is the following weak inf-sup condition for b. Lemma 1. There exist constants C1, C2 > 0, independent of hK and ν, such that Proof. We consider an arbitrary q h ∈ Q k h . LetΩ be a convex, open, Lipschitz set such that Ω ⊂Ω, and let us consider following extension Let now φ be the unique weak solution of the problem SinceΩ is convex, then φ ∈ H 2 (Ω). Then w := ∇φ|Ω belongs to [H 1 (Ω)] 2 , and forw := wt, In addition, applying standard regularity results, see [ Let (wh,w h ) := Π (w), then thanks to (2.6), (2.4) and the continuity of b (see [1,Lemma 3 Using the approximation properties of the BDM interpolation operator (see [ where, in the last estimate we have used the stability of the Fortin operator Π in the ||| · ||| norm (2.7). This proves the result with C1 = Before showing an inf-sup condition, we prove the continuity of bilinear form A. Lemma 2. There exists a constant C > 0 such that, for all (wh,w h ) , (vh,ṽ h ) ∈ Vh and r h , q h ∈ Q k h , we have Proof. We use the continuity of the bilinear forms (see [1,Lemma 3.3]) and the fact that the projection is a bounded operator.
The final step towards stability is proving the inf-sup condition for bilinear form A.
Lemma 3. There exists β > 0 independent of hK such that for all (wh,w h , r h ) ∈ Vh × Q k h the following holds As a consequence, Problem (2.2) is well-posed.
Proof. Let (wh,w h , r h ) ∈ Vh × Q k h . The idea of the proof is to construct an appropriate (vh,ṽ h , q h ) such that To achieve that we use coercivity of a (see [1,Lemma 3.4]), continuity of a (see [1,Lemma 3.3]) and Lemma 2. For details see [6].

Error analysis
In this section we present the error estimates for the method. The addition of the stabilising bilinear form s(·, ·) introduced a consistency error. However according to [4], this should not be viewed as a serious flaw, as this consistency error can be bounded in an optimal way. The following result is the first step towards that goal.
(Ω) be a solution of (1.1),ũ = ut on all edges in E h , and (uh,ũ h , p h ) ∈ Vh × Q k h solves the discrete problem (2.2). Then there exists C > 0, independent of h and ν, such that Proof. It is a combination of Lemmas 1, 2 and 3. For details see [6].
be a solution of (1.1),ũ = ut on all edges in E h , and (uh,ũ h , p h ) ∈ Vh × Q k h solves the discrete problem (2.2). Then there exists C > 0, independent of h and ν, such that Proof. It is a combination of [1, Lemmas 3.8] and Lemma 5 with the local L 2 -projection approximation [11,Theorem 1.103].

Numerical experiments
The computational domain is the unit square Ω = (0, 1) 2 . We present the results for k = 1, that is the discrete space is given by BDM 1 h ×M 0 h,0 ×Q 1 h . We test both the symmetric method (ε = −1) and the non-symmetric method (ε = 1). We have followed the recommendation given in [15,Section 2.5.2] and taken τ = 6. We choose the right hand side f and the boundary condition g such that the exact solution is given by p = tan(xy). In Figure 1a and 1b we depict the errors for both the symmetric and non-symmetric cases, respectively. We can see that they not only validate the theory from Section 2.2, but also perform an optimal h 2 convergence rate for u − uh Ω. Furthermore, we observe an increased order of convergence for p − p h Ω. In fact, the error seems to decrease with O(h 3/2 ), rather than the O(h) predicted by the theory.  To stress the last point made in the previous paragraph, in Table 1 we compare the L 2 error of the pressure (||p − p h ||Ω) for hdG method introduced in [1] and stabilised hdG method from Section 2. Columns p h ∈ Q 0 h are associated with hdG method and p h ∈ Q 1 h with stabilised hdG ones. There, we confirm that the pressure error for the stabilised version is much smaller than the one for the inf-sup stable case, in addition to having an increased order of convergence.

Conclusion
In this work we have applied the idea introduced in [10] to stabilise the hdG method proposed in [1] for the Stokes problem with TVNF boundary conditions. The method adds a simple, symmetric, term to the formulation, and allowed us to use a higher order pressure space, which, in turn, improved the pressure convergence (although a proof of this fact is, in general, not available). This approach was also applied to NVTF boundary conditions (see [6]) and can be used for other discontinuous Galerkin methods that deal with Stokes or nearly incompressible elasticity problems. Future testing using higher order discretisations is needed to assess whether this approach provides an increase of the convergence rate for the pressure. Thus, the numerical tests with higher order of polynomials for discontinuous finite methods is interest for further research to look for the improvement of the convergence.