A low order divergence-free H(div)-conforming finite element method for Stokes flows

Xu Li School of Mathematics, Shandong University, Jinan 250100, China xulisdu@126.com and Hongxing Rui School of Mathematics, Shandong University, Jinan 250100, China hxrui@sdu.edu.cn

Abstract.

In this paper, we propose a ${\bm{P}_{1}^{c}}\oplus{RT0}-P0$ discretization of the Stokes equations on general simplicial meshes in two/three dimensions (2D/3D), which yields an exactly divergence-free and pressure-independent velocity approximation with optimal order. Our method has the following features. Firstly, the global number of the degrees of freedom of our method is the same as the low order Bernardi and Raugel ( $B$ - $R$ ) finite element method (Bernardi and Raugel, 1985), while the number of the non-zero entries of the former is about half of the latter in the velocity-velocity region of the coefficient matrix. Secondly, the ${\bm{P}_{1}^{c}}$ component of the velocity, the $RT0$ component of the velocity and the pressure seem to solve a popular ${\bm{P}_{1}^{c}}-{RT0}-P0$ discretization of a poroelastic-type system formally. Finally, our method can be easily transformed into a pressure-robust and stabilized ${\bm{P}_{1}^{c}}-P0$ discretization for the Stokes problem via the static condensation of the $RT0$ component, which has a much smaller number of global degrees of freedom. Numerical experiments illustrating the robustness of our method are also provided.

Key words and phrases:

Stokes problem; divergence-free velocity; pressure-robustness; optimal error estimates.

The second author is the corresponding author.

1. Introduction

The importance of strong mass conservation has been found in a variety of real-world applications (see, e.g., [36, 17, 34]). For the incompressible Stokes problem, pointwise mass conservation requires an exactly divergence-free velocity approximation. Another highly relevant concept is pressure-robustness (see, e.g., [25, 35, 18]), which describes a type of robustness that the velocity error is independent of the pressure. In this paper, we design a low order divergence-free and pressure-robust $H(\operatorname{div})$ -conforming mixed method for the Stokes problem on general shape-regular simplicial meshes.

Due to the fascinating properties of divergence-free mixed finite element methods [16, 25], constructing efficient divergence-free methods has been a hot research topic during recent years and great progress has been made. However, this is still not an easy thing especially in three dimensions. For conforming divergence-free methods, we refer the readers to [46, 2, 50, 21, 23, 9, 15, 16], to name just a few. Among these methods, the lowest order cases on simplicial meshes might be the ones proposed in [9] and [23], where they have the same number of global degrees of freedom (DOFs) as the classical Bernardi and Raugel finite element method [4]. Like the Bernardi and Raugel element, they both use $d+1$ face bubbles to enrich the linear polynomial space on each simplex in $d$ dimensions. The difference lies in the construction of the face bubbles: the bubbles in [9] are piecewise linear on the Powell-Sabin split (e.g., in 3D, a tetrahedron is split into $12$ small tetrahedrons) of a simplex, and the bubbles in [23] are typically piecewise polynomials of degree $d$ on the barycentric refinement of a simplex. Some applications of the former could be found in [8]. Another popular class of mixed methods is based on $H(\operatorname{div})$ -conforming elements, which are a class of nonconforming elements for the Stokes equations. In general, there are two ways to stabilize the nonconformity (tangential discontinuity) of the elements. One is to modify the velocity-velocity bilinear form into a discontinuous Galerkin (DG) framework (see, e.g., [12, 48, 27, 26]), the other is to modify the discrete velocity space to impose tangential continuity in some weak sense such as [22, 39, 49, 47]. Among these methods the lowest order case on simplicial meshes might be the DG method based on $BDM1-P0$ pair [12, 48] ( $BDM1$ : the lowest order Brezzi-Douglas-Marini element [7], $P0$ : the piecewise constant space with zero mean). Recently, there is also a class of pressure-robust discretizations developed on classical non-divergence-free mixed methods via divergence-free reconstruction of test functions (cf. [37, 38, 35, 30]), which can be seen as a class of Petrov-Galerkin methods.

One of the advantages of DG $H(\operatorname{div})$ -conforming mixed methods lies in the simple construction of their velocity spaces. For instance, the lowest order velocity is piecewise linear in arbitrary dimensions. However, in a DG framework we should introduce some integral terms over element interfaces which may increase the cost and enlarge the stencil (sparsity). In addition, the stabilization parameters (the parameters of the jump-penalty term) need to be chosen carefully to guarantee the stability in some cases. The features of conforming and the second type of $H(\operatorname{div})$ -conforming divergence-free mixed methods are exactly opposite to the DG methods mentioned above. They do not need to modify the bilinear form and are always stable. The constructions of their spaces are, however, a little complicated and high order polynomials are usually included. In a word, constructing divergence-free mixed methods is non-trivial.

The aim of this paper is to investigate a way on how to combine some advantages of the conforming and (DG-) $H(\operatorname{div})$ -conforming divergence-free mixed methods. In our method, the velocity space is obtained by enriching the continuous vector-valued piecewise linear polynomial space ( ${\bm{P}_{1}^{c}}$ ) with the lowest order Raviart-Thomas-Nédélec space ( $RT0$ ) [40, 42], and the pressure is approximated by piecewise constants with zero mean. In a sense, our method can be seen as replacing the Bernardi and Raugel bubbles [4] with the lowest order Raviart-Thomas-Nédélec elements (see Fig. 1). This replacement does bring some advantages. For example, the discrete velocity becomes pointwise divergence-free. However, this turns into a type of nonconforming method and some stabilization strategies are needed. To resolve this issue,

Refer to caption — Figure 1. DOFs of $B$ - $R$ (left), $RT0$ (middle) and ${\bm{P}_{1}^{c}}$ (right) element in 2D. Solid circles indicate function evaluation and arrows indicate normal component evaluation.

we propose a conforming-like formulation which does not involve any integrals over element interfaces. The resulting method has several advantages. Firstly, the construction of the space is very simple in any dimensions. Secondly, we can use the lowest order quadrature rule to compute each entry of the coefficient matrix exactly in some cases. Thirdly, our method instead has much better sparsity than the Bernardi and Raugel method. Next, the discrete system is stable as long as the parameters introduced are positive. Finally, we can further obtain a stabilized pressure-robust ${\bm{P}_{1}^{c}}-P0$ discretization for the Stokes equations via static condensation. Compared to a direct ${\bm{P}_{1}^{c}}-P0$ discretization (unstable), this only changes the pressure region of the right-hand side and the pressure-pressure region of the coefficient matrix, which is different from the similar strategy for the Bernardi and Raugel element in [45]. The price to pay for these advantages is a consistency error bounded in an optimal order.

The motivation for studying low order divergence-free elements is threefold: (1) the low order element in this paper is easy to implement; (2) in some practical applications, the available mesh might be fixed due to the complex geometry [24], in which case lower order methods allow a smaller global system; (3) the divergence-free nature usually makes low order elements also give a satisfactory result for some non-trivial Stokes equations such as those with large pressure gradient, small viscosity or long time simulation (unsteady case), cf. [25, 1].

We also want to mention that hybridization is a very efficient way to reduce the computational cost of DG methods, cf. [11]. In [32] and [33], a class of $H(\operatorname{div})$ -conforming hybrid DG (HDG) methods with discontinuous pressure elements was presented. A cheaper version of them with relaxed $H(\operatorname{div})$ -conformity for velocity and pressure-robust reconstructions could be found in [28] and [29], where the lowest order pair (after static condensation) has one DOF per facet and dimension for velocity and one DOF per element for pressure, the same as $CR1-P0$ pair ( $CR1$ : the first order nonconforming Crouzeix-Raviart element, cf. [13]). A total hybridization (simultaneously for velocity and pressure) formulation for the Stokes equations was analyzed in [43], where the velocity is automatically $H(\operatorname{div})$ -conforming if the order of the pressure space is one less than the order of the velocity space. To further lower the cost, an embedded DG (EDG) method with pressure-robust reconstruction and an embedded-hybridized DG (EDG-HDG, i.e., EDG for velocity and HDG for pressure) method were developed in [31] and [44], respectively. EDG methods use continuous facet unknowns which have the same number of DOFs as the corresponding conforming methods, while retaining the versatility of DG methods. In a sense, our method shares a similar topic with EDG methods, i.e., combining some advantages of conforming methods and DG methods. The number of the DOFs (in the reduced version) of the lowest order pair in [31] is equal to the one of ${\bm{P}_{1}^{c}}-P_{1}^{c}$ pair, where $P_{1}^{c}$ is the scalar version of ${\bm{P}_{1}^{c}}$ . This is also equal to that of the lowest order MINI elements with static condensation of the cell bubbles, which can also be pressure-robust with reconstructions, cf. [30]. The two low order pairs have less DOFs than ${\bm{P}_{1}^{c}}-P0$ due to a smaller pressure space and can be easily extended to higher order cases. Compared to them, the advantage of our method lies in its automatically divergence-free nature.

The rest of the paper is organized as follows. Section 2 is devoted to describing the model and its discretizations. In Section 3, we analyze the LBB stability, consistency and error estimates. Some numerical studies are shown in Section 4. Finally we do some conclusions in Section 5.

Throughout the paper we use $C,$ with or without subscript, to denote a generic positive constant. The standard inner product and norm (seminorm) of the Sobolov space $[H^{m}(X)]^{n}$ or $[H^{m}(X)]^{n\times n}$ ( $n\in\mathbb{N}^{+}$ ) are denoted by $(\cdot,\cdot)_{m,X}$ and $\|\cdot\|_{m,X}$ ( $|\cdot|_{m,X}$ ), respectively. When $m=0$ ( $X=\Omega$ ), with the convention that the index $m$ ( $X$ , respectively) is omitted. For any sub-dimensional face $e$ , we have $\langle\cdot,\cdot\rangle_{e}$ and $||\cdot||_{e}$ similarly.

2. Model, notation and method

In this section we consider the incompressible Stokes equations in a bounded domain $\Omega\subset\mathbb{R}^{d}~(d=2,3)$


(2.1a)	$\displaystyle-\nu\Delta\bm{u}+\nabla p$	$\displaystyle=\bm{f}\quad\text{ in }\Omega,$
(2.1b)	$\displaystyle\nabla\cdot\bm{u}$	$\displaystyle=0\,\quad\text{ in }\Omega,$
(2.1c)	$\displaystyle\bm{u}$	$\displaystyle={\bm{0}}\quad\text{ on }\partial\Omega,$

where the unknowns $(\bm{u},p)$ represent the velocity and pressure, respectively; $\nu$ denotes the constant fluid viscosity; $\bm{f}$ is the external body force; Eq. 2.1a is usually called the momentum equation and Eq. 2.1b is usually called the continuity equation. For simplicity, we assume that $\Omega$ has a Lipschitz-continuous polygonal/polyhedral boundary $\Gamma$ .

Let $V=[H_{0}^{1}(\Omega)]^{d}$ and $W=L_{0}^{2}(\Omega)$ . A commonly used variational form of Eq. 2.1 is:


	$\displaystyle{\rm find}~(\bm{u},p)\in V\times W$	$\displaystyle~{\rm such~that}$
(2.2a)		$\displaystyle\nu a(\bm{u},\bm{v})-b(\bm{v},p)$	$\displaystyle=(\bm{f},\bm{v})\quad\forall~\bm{v}\in V,$
(2.2b)		$\displaystyle b(\bm{u},q)$	$\displaystyle=0\qquad\quad\forall~q\in W,$

where

a(\bm{u},\bm{v})=(\nabla\bm{u},\nabla\bm{v}),\quad b(\bm{v},q)=(\nabla\cdot\bm{v},q).

And the spaces $(V,W)$ satisfy the following well-known inf-sup condition (See [5, §4.2])

(2.3)

\sup_{\bm{v}\in V\setminus\{\bm{0}\}}\frac{b(\bm{v},q)}{||\bm{v}||_{1}}\geq\beta\|q\|\quad\forall~q\in W,

where $\beta$ is a positive constant.

2.1. Discretizations

Let $\mathcal{T}_{h}$ be shape-regular triangular or tetrahedral partitions of $\Omega$ [10]. Let $h_{T}$ and $h_{e}$ denote the diameters of elements $T$ and faces $e$ , respectively, and $h=\max_{T\in\mathcal{T}_{h}}h_{T}$ . We denote the set of interior faces of $\mathcal{T}_{h}$ by $\mathcal{E}^{0}$ , the set of boundary faces by $\mathcal{E}^{\partial}$ and $\mathcal{E}=\mathcal{E}^{0}\cup\mathcal{E}^{\partial}$ . We define

H(\operatorname{div};\Omega)=\left\{\bm{v}\in[L^{2}(\Omega)]^{d}:\nabla\cdot\bm{v}\in L^{2}(\Omega)\right\},\quad H_{0}(\operatorname{div};\Omega)=\left\{\bm{v}\in H(\operatorname{div};\Omega):\left.\bm{v}\cdot\bm{n}\right|_{\Gamma}=0\right\}.

Let $P_{k}(T)$ ( $k\geq 0$ ) denote the space of polynomials of degree no more than $k$ on elements $T$ , and the finite element spaces are set as

V_{h}^{1}=\{\bm{v}\in V:\bm{v}|_{T}\in[P_{1}(T)]^{d}\quad\forall\ T\in\mathcal{T}_{h}\},

(2.4)

RT0=\{\bm{v}\in H_{0}(\operatorname{div};\Omega):\bm{v}|_{T}\in[P_{0}(T)]^{d}\oplus\bm{x}P_{0}(T)\quad\forall\ T\in\mathcal{T}_{h}\},

W_{h}=\left\{q\in L_{0}^{2}(\Omega):\left.q\right|_{T}\in P_{0}(T)\quad\forall~T\in\mathcal{T}_{h}\right\}.

We note that the functions in $V_{h}^{1}$ are totally continuous, while the functions in $RT0$ may have discontinuous tangential components across the interelement boundaries. The following lemma depicts a simple but important fact in our method.

Lemma 1.

$V_{h}^{1}\cap RT0=\{\bf 0\}$ .

Proof.

Consider a simplest case in two dimensions (see Fig. 2).

Two axes are respectively labeled as $x$ and $y$ . From 2.4, a function $\bm{f}^{rt}\in RT0$ can be written as follows,

\bm{f}^{rt}|_{T_{1}}=\begin{pmatrix}a_{1}+c_{1}x,b_{1}+c_{1}y\end{pmatrix}^{\top},\quad\bm{f}^{rt}|_{T_{2}}=\begin{pmatrix}a_{2}+c_{2}x,b_{2}+c_{2}y\end{pmatrix}^{\top},\quad a_{i},b_{i},c_{i}\in\mathbb{R},\quad i=1,2.

It suffices to prove that $\bm{f}^{rt}\in V_{h}^{1}\Rightarrow\bm{f}^{rt}=\bm{0}$ . If $\bm{f}^{rt}\in V_{h}^{1}$ , $\bm{f}^{rt}$ is totaly continuous across the element interface $\ell$ . In this case, we have


(2.5a)		$\displaystyle a_{1}+c_{1}x=a_{2}+c_{2}x\quad\text{on}~\ell,$
(2.5b)		$\displaystyle b_{1}+c_{1}y=b_{2}+c_{2}y\quad\text{on}~\ell.$

We note that there is at least one of the variables ( $x$ or $y$ ) varying over $\ell$ . For brevity, we assume that $x$ varies over $\ell$ . From 2.5a we can obtain

a_{1}=a_{2},\quad c_{1}=c_{2}.

Substituting $c_{1}=c_{2}$ into 2.5b gives $b_{1}=b_{2}.$ Then we have

\bm{f}^{rt}|_{\Omega}=\begin{pmatrix}a_{1}+c_{1}x,b_{1}+c_{1}y\end{pmatrix}^{\top}.

If we further demand that $\bm{f}^{rt}|_{\Gamma}=\bm{0}$ , through the similar analysis, one can obtain $a_{1}=b_{1}=c_{1}=0.$ Thus we complete the proof of this simplest case. The proofs of general decompositions and three-dimensional case are a natural extension of this. ∎

Denote by $V_{h}=V_{h}^{1}\oplus RT0$ the discrete velocity space. Let $\mathcal{RT}:=\{\bm{\Phi}_{e},~e\in\mathcal{E}^{0}\}$ denote the set of the commonly used bases in $RT0$ , where $\bm{\Phi}_{e}$ is only supported on two elements $T_{1}$ and $T_{2}$ which share the face $e$ and $\bm{\Phi}_{e}\cdot\bm{n}_{e^{\prime}}$ vanishes on any face $e^{\prime}\neq e$ . Here $\bm{n}_{e^{\prime}}$ denotes the unit normal vector to $e^{\prime}$ .

For any $\bm{v}_{h}\in V_{h}$ , we have the unique decompositions

\bm{v}_{h}=\bm{v}_{h}^{1}+\bm{v}_{h}^{R},\quad\bm{v}_{h}^{R}=\sum_{e\in\mathcal{E}^{0}}v_{e}\bm{\Phi}_{e}\quad\text{ with }\bm{v}_{h}^{1}\in V_{h}^{1},\ \bm{v}_{h}^{R}\in RT0.

Our discretization is presented as:


	$\displaystyle{\rm Find}~(\bm{u}_{h},p_{h})\in V_{h}\times W_{h}~$	$\displaystyle{\rm such~that}$
(2.6a)		$\displaystyle\nu a_{h}(\bm{u}_{h},\bm{v}_{h})-b(\bm{v}_{h},p_{h})$	$\displaystyle=(\bm{f},\bm{v}_{h})\quad\forall~\bm{v}_{h}\in V_{h},$
(2.6b)		$\displaystyle b(\bm{u}_{h},q_{h})$	$\displaystyle=0\qquad\ \,\quad\forall~q_{h}\in W_{h},$

where

(2.7)

a_{h}(\bm{u}_{h},\bm{v}_{h})=a(\bm{u}_{h}^{1},\bm{v}_{h}^{1})+a^{R}(\bm{u}_{h}^{R},\bm{v}_{h}^{R}).

Here $a^{R}:RT0\times RT0\rightarrow R$ has several choices as

(2.8)

a^{R}(\bm{u}_{h}^{R},\bm{v}_{h}^{R})=a^{0}(\bm{u}_{h}^{R},\bm{v}_{h}^{R}):=\sum_{T\in\mathcal{T}_{h}}\alpha_{T}h_{T}^{-2}(\bm{u}_{h}^{R},\bm{v}_{h}^{R})_{T},

(2.9)

a^{R}(\bm{u}_{h}^{R},\bm{v}_{h}^{R})=a^{D}(\bm{u}_{h}^{R},\bm{v}_{h}^{R}):=\sum_{T\in\mathcal{T}_{h}}\sum_{e\in\partial T\cap\mathcal{E}^{0}}\alpha_{T}h_{T}^{-2}u_{e}v_{e}(\bm{\Phi}_{e},\bm{\Phi}_{e})_{T},

(2.10)

a^{R}(\bm{u}_{h}^{R},\bm{v}_{h}^{R})=a^{\operatorname{div}}(\bm{u}_{h}^{R},\bm{v}_{h}^{R}):=\sum_{T\in\mathcal{T}_{h}}\sum_{e\in\partial T\cap\mathcal{E}^{0}}\alpha_{T}u_{e}v_{e}(\nabla\cdot\bm{\Phi}_{e},\nabla\cdot\bm{\Phi}_{e})_{T},

where $\alpha_{T}$ , $T\in\mathcal{T}_{h}$ are positive parameters.

Note that

a^{0}(\bm{u}_{h}^{R},\bm{v}_{h}^{R})=\sum_{T\in\mathcal{T}_{h}}\sum_{e\in\mathcal{E}^{0}\cap\partial T}\sum_{e^{\prime}\in\mathcal{E}^{0}\cap\partial T}\alpha_{T}h_{T}^{-2}u_{e}v_{e^{\prime}}(\bm{\Phi}_{e},\bm{\Phi}_{e^{\prime}})_{T},

and

a^{D}(\bm{u}_{h}^{R},\bm{v}_{h}^{R})=\sum_{T\in\mathcal{T}_{h}}\sum_{e\in\mathcal{E}^{0}\cap\partial T}\sum_{e^{\prime}\in\mathcal{E}^{0}\cap\partial T}{\delta}_{ee^{\prime}}\alpha_{T}h_{T}^{-2}u_{e}v_{e^{\prime}}(\bm{\Phi}_{e},\bm{\Phi}_{e^{\prime}})_{T},

where ${\delta}_{ee^{\prime}}=1$ if $e=e^{\prime}$ , otherwise vanishes. Thus, $a^{D}(\cdot,\cdot)$ can be regarded as a diagonalization version of $a^{0}(\cdot,\cdot)$ . Clearly the matrices from $a^{D}(\cdot,\cdot)$ and $a^{\operatorname{div}}(\cdot,\cdot)$ are both diagonal. In the next section, we will prove that the three choices of $a^{R}(\cdot,\cdot)$ are equivalent and formula 2.6 is stable as long as $\alpha_{T}>0$ for all $T\in\mathcal{T}_{h}$ .

Remark 1 (The motivation of $a^{R}$ ).

The three forms of $a^{R}$ proposed here can also be regarded as some penalty-type stabilizations. From Lemma 1 we can further obtain $RT0\cap V=\{\bm{0}\}$ . In other words, there is no conforming element in $RT0$ . Thus, we penalize the nonconformity in another way. That is, we penalize the whole $RT0$ component, and enforce $\bm{u}_{h}^{R}=\bm{0}$ when $h\rightarrow 0$ .

The formulation 2.6 can be rewritten to the following equivalent three-field form:


	$\displaystyle{\rm Find}~(\bm{u}_{h}^{1},\bm{u}_{h}^{R},p_{h})\in V_{h}^{1}\times RT0\times W_{h}~$	$\displaystyle{\rm such~that}$
(2.11a)		$\displaystyle\nu a(\bm{u}_{h}^{1},\bm{v}_{h}^{1})-b(\bm{v}_{h}^{1},p_{h})$	$\displaystyle=(\bm{f},\bm{v}_{h}^{1})~\quad\forall~\bm{v}^{1}_{h}\in V_{h}^{1},$
(2.11b)		$\displaystyle\nu a^{R}(\bm{u}_{h}^{R},\bm{v}_{h}^{R})-b(\bm{v}_{h}^{R},p_{h})$	$\displaystyle=(\bm{f},\bm{v}_{h}^{R})\quad\forall~\bm{v}^{R}_{h}\in RT0,$
(2.11c)		$\displaystyle b(\bm{u}_{h}^{1},q_{h})+b(\bm{u}_{h}^{R},q_{h})$	$\displaystyle=0\qquad\ \,\quad\forall~q_{h}\in W_{h}.$

Remark 2.

Eq. 2.11a is a ${\bm{P}_{1}^{c}}$ discretization of the momentum equation 2.1a, while 2.11b is very similar to a low order discretization of the Darcy equation with the permeability $\kappa=(\nu\alpha_{T}h_{T}^{-2})^{-1}$ on the element $T$ , when $a^{0}$ is used. Thus, our method is very similar to a ${\bm{P}_{1}^{c}}-RT0-P0$ discretization of a three-field (displacement - (Darcy) velocity - pressure) poroelastic-type system [41].

Since $\nabla\cdot RT0=W_{h}$ (see [5, Proposition 2.3.3]) and $\nabla\cdot V_{h}^{1}\subset W_{h}$ , we have $\nabla\cdot V_{h}=W_{h}$ . Thus, from 2.6b we know

Theorem 2.1 (Mass Conservation).

The finite element solution $\bm{u}_{h}$ of 2.6 features a full satisfaction of the continuity equation, which means

(2.12)

\nabla\cdot\bm{u}_{h}\equiv 0.

In other words, our method conserves mass pointwise.

Proof.

Taking $q_{h}=\nabla\cdot\bm{u}_{h}$ in 2.6b gives $||\nabla\cdot\bm{u}_{h}||^{2}=0.$ Then 2.12 follows. ∎

In the continuous case, there is a fundamental invariance property: Changing the external force by a gradient field changes only the pressure solution, and not the velocity; in symbols,

\bm{f}\rightarrow\bm{f}+\nabla\psi\quad\Longrightarrow\quad(\bm{u},p)\rightarrow(\bm{u},p+\psi),

where $\psi\in H^{1}(\Omega)/\mathbb{R}$ . This phenomenon is highly relevant to pressure-robustness [25]. In our method, we have the following discrete invariance property.

Theorem 2.2.

The following invariance property holds in the discrete scheme 2.6:

(2.13)

\bm{f}\rightarrow\bm{f}+\nabla\psi\quad\Longrightarrow\quad(\bm{u}_{h},p_{h})\rightarrow(\bm{u}_{h},p_{h}+P_{h}\psi),

where $P_{h}:W\rightarrow W_{h}$ is the $L^{2}$ projection (see 3.1).

Proof.

Since $\nabla\cdot V_{h}=W_{h}$ , we have

(2.14)

(p_{h}+P_{h}\psi,\nabla\cdot\bm{v}_{h})=(p_{h}+\psi,\nabla\cdot\bm{v}_{h})=(p_{h},\nabla\cdot\bm{v}_{h})-(\nabla\psi,\bm{v}_{h}).

Substituting 2.14 into 2.6 gives


(2.15a)	$\displaystyle a_{h}(\bm{u}_{h},\bm{v}_{h})-b(\bm{v}_{h},p_{h}+P_{h}\psi)$	$\displaystyle=(\bm{f}+\nabla\psi,\bm{v}_{h})\quad\forall~\bm{v}_{h}\in V_{h},$
(2.15b)	$\displaystyle b(\bm{u}_{h},q_{h})$	$\displaystyle=0\qquad\quad\quad\quad\,\quad\forall~q_{h}\in W_{h},$

which implies that $(\bm{u}_{h},p_{h}+P_{h}\psi)$ is a solution corresponding to $\bm{f}+\nabla\psi$ . Together with the unique solvability of our methods (see next section), the invariance property 2.13 follows. ∎

2.2. A pressure-robust and stabilized ${P_{1}^{c}}-P0$ scheme

A block form of the fully discrete problem 2.6 or 2.11 can be written as follows,

(2.16)

\mathcal{A}_{S}\begin{pmatrix}U_{L}\\ U_{R}\\ P\end{pmatrix}=\begin{pmatrix}F_{L}\\ F_{R}\\ \bm{0}\end{pmatrix},\quad\text{ with }\mathcal{A}_{S}=\begin{pmatrix}A_{LL}&0&-G_{L}^{T}\\ 0&A_{RR}&-G_{R}^{T}\\ G_{L}&G_{R}&0\end{pmatrix},

where ${U}_{L},{U}_{R},$ and ${P}$ are the unknown vectors corresponding to the continuous linear component of the velocity, the $RT0$ component of the velocity, and the pressure, respectively. The blocks $A_{LL}$ , $A_{RR}$ , $G_{L}$ , $G_{R}$ , $F_{L}$ and $F_{R}$ correspond to the bilinear forms $a(\bm{u}_{h}^{1},\bm{v}_{h}^{1})$ , $a^{R}(\bm{u}_{h}^{R},\bm{v}_{h}^{R})$ , $b(\bm{u}_{h}^{1},q_{h})$ , $b(\bm{u}_{h}^{R},q_{h})$ , $(\bm{f},\bm{v}_{h}^{1})$ and $(\bm{f},\bm{v}_{h}^{R})$ , respectively.

For a classical Bernardi and Raugel discretization, the cross region of the linear component and the bubble component of the velocity in the coefficient matrix is non-zero. In other words, the matrices from $a(\bm{u}_{h}^{1},\bm{v}_{h}^{b})$ and $a(\bm{u}_{h}^{b},\bm{v}_{h}^{1})$ are not zero matrices, where $\bm{u}_{h}^{b}$ and $\bm{v}_{h}^{b}$ denote the bubble parts. Thus, the number of the non-zero entries in the velocity-velocity region from our method is about half of the Bernardi and Raugel finite element method.

When $a^{D}$ or $a^{\operatorname{div}}$ is used, the block $A_{RR}$ is diagonal. In this situation, the component $\bm{u}_{h}^{R}$ can be locally expressed in terms of the discrete pressure $p_{h}$ . In other words, the inverse of $A_{RR}$ is easy to obtain and also diagonal (sparse). It is then easy to eliminate the unknowns corresponding to $RT0$ via static condensation, which leads to the following system:

\widehat{\mathcal{A}}_{S}^{D}\begin{pmatrix}U_{L}\\ P\end{pmatrix}=\begin{pmatrix}F_{L}\\ -G_{R}A_{RR}^{-1}F_{R}\end{pmatrix}\quad\text{ with }\widehat{\mathcal{A}}_{S}^{D}=\begin{pmatrix}A_{LL}&-G_{L}^{T}\\ G_{L}&G_{R}A_{RR}^{-1}G_{R}^{T}\end{pmatrix}.

Since $\bm{u}_{h}$ is pressure-robust, $\bm{u}_{h}^{1}$ is also pressure-robust (see the error estimate section below). In the next section we will prove that although $\bm{u}_{h}^{1}$ is not divergence-free, it has optimal convergence properties in $L^{2}$ and $H^{1}$ norms. Thus, the above system can be regarded as a stabilized and pressure-robust ${\bm{P}_{1}^{c}}-P0$ discretization of the Stokes equations.

Remark 3.

When $a^{R}=a^{\operatorname{div}}$ , all the terms in the left-hand side of 2.16 can be computed out exactly by the lowest order quadrature rule (barycentric quadrature rule). And the system can be reduced to a ${\bm{P}_{1}^{c}}-P0$ system. Thus, we prefer the bilinear form $a^{\operatorname{div}}$ in practice.

Remark 4.

To obtain a divergence-free solution from the above ${\bm{P}_{1}^{c}}-P0$ system, one can get $\bm{u}_{h}^{1}$ and $p_{h}$ first, then compute $\bm{u}_{h}^{R}$ by $p_{h}$ . Since $A_{RR}$ is diagonal, the computation of $\bm{u}_{h}^{R}$ is very low-cost. Compared to a direct ${\bm{P}_{1}^{c}}-P0$ discretization (unstable) of the Stokes equations, our method only changes the pressure-pressure region (lower right) of the coefficient matrix and the pressure region of the right-hand side.

3. Error estimates

3.1. Preliminaries

We define the standard Raviart-Thomas interpolation operator $\Pi_{h}^{R}:V\rightarrow RT0$ and the piecewise integral mean operator $P_{h}:W\rightarrow{W}_{h}$ by

\left(\left(\bm{v}-{\Pi}_{h}^{R}\bm{v}\right)\cdot\bm{n}_{e},1\right)_{e}=0\quad\forall~e\in\mathcal{E},

and

(3.1)

\left(r-P_{h}r,q\right)=0\quad\forall~q\in{W}_{h},

respectively.

Lemma 2.

The operators $\Pi_{h}^{R}$ and $P_{h}$ satisfy the following properties:

(3.2)

\|\bm{v}-\Pi_{h}^{R}\bm{v}\|_{T}\leq C_{R}h_{T}|\bm{v}|_{1,T}\quad\forall~T\in\mathcal{T}_{h},

(3.3)

\left\|r-P_{h}r\right\|_{T}\leq C_{P}h_{T}^{m}|r|_{m,T}\quad\forall~T\in\mathcal{T}_{h},\quad m=0,1,

and

(3.4)

\nabla\cdot\Pi_{h}^{R}\bm{v}=P_{h}\nabla\cdot\bm{v},

for any $\bm{v}\in V$ and $r\in W\cap H^{m}(\Omega)$ , where $C_{R}$ and $C_{P}$ depend on the shape regularity of $\mathcal{T}_{h}$ .

Proof.

See [5, §2.5]. ∎

Next, let $V(h)=V\oplus RT0.$ be a larger space containing $V$ and $V_{h}$ . Recall that $V=\left[H_{0}^{1}(\Omega)\right]^{d}$ . For any $\bm{v}\in V(h)$ , we have the unique decomposition $\bm{v}=\bm{v}^{1}+\bm{v}^{R}$ with $\bm{v}^{1}\in{V},\bm{v}^{R}\in RT0.$ An extension of the definition of $a_{h}(\cdot,\cdot)$ , 2.7, to the whole $V(h)\times V(h)$ is needed prior to analysis. We define that, for any $\bm{w},\bm{v}\in V(h)$ ,

(3.5)

a_{h}(\bm{w},\bm{v}):=a(\bm{w}^{1},\bm{v}^{1})+a^{R}(\bm{w}^{R},\bm{v}^{R}).

Restricted to the finite element space $V_{h}\times V_{h}$ , the definition 3.5 coincides with the definition 2.7. And on $V\times V$ we have $a_{h}(\cdot,\cdot)|_{V\times V}=a(\cdot,\cdot).$ Thus the definition of $a_{h}(\cdot,\cdot)$ is logical and natural.

To investigate the properties of the discretization formula 2.6, we introduce a norm $|||\cdot|||$ for the set $V(h)$ as follows,

(3.6)

|||\bm{v}|||^{2}=a_{h}(\bm{v},\bm{v})\quad\forall~\bm{v}\in V(h).

In this definition, for any $\bm{v}^{1}\in{V}$ , we have $|||\bm{v}^{1}|||=|\bm{v}^{1}|_{1}.$

For the case that $a^{R}=a^{0}$ , it is not difficult to verify that $\|\bm{v}\|\leq\|\bm{v}^{1}\|+\|\bm{v}^{R}\|\leq|\bm{v}^{1}|_{1}+\|\bm{v}^{R}\|\leq C|||\bm{v}|||$ for any $\bm{v}\in V(h)$ . Thus, 3.6 really defines a norm for $a^{R}=a^{0}$ . For the other two choices of $a^{R}$ , that is, $a^{D}$ and $a^{\operatorname{div}}$ , we have the following results.

Lemma 3.

For any $\bm{v}^{R}\in RT0$ , we have

(3.7)

C_{*}a^{0}(\bm{v}^{R},\bm{v}^{R})\leq a^{D}(\bm{v}^{R},\bm{v}^{R})\leq C^{*}a^{0}(\bm{v}^{R},\bm{v}^{R}),

and

(3.8)

\gamma_{*}a^{\operatorname{div}}(\bm{v}^{R},\bm{v}^{R})\leq a^{D}(\bm{v}^{R},\bm{v}^{R})\leq\gamma^{*}a^{\operatorname{div}}(\bm{v}^{R},\bm{v}^{R}),

where the constant $C_{*}$ only depends on the parameters $\alpha_{T},T\in\mathcal{T}_{h}$ , and $C^{*}$ , $\gamma_{*}$ , $\gamma^{*}$ only depend on the shape regularity of $\mathcal{T}_{h}$ and the parameters.

Proof.

For simplicity, we will only prove the lemma in the case that $\alpha_{T}\equiv\alpha$ for all $T\in\mathcal{T}_{h}$ , where $\alpha$ is an arbitrary but fixed positive constant. Let us prove the first inequality in 3.7. Applying the same technique in [45, Lemma 4.2], by the Cauchy-Schwarz inequality we can similarly obtain

(3.9)

\left(\bm{v}^{R},\bm{v}^{R}\right)_{T}\leq(d+1)\sum_{e\in\mathcal{E}^{0}\cap\partial T}v_{e}^{2}(\bm{\Phi}_{e},\bm{\Phi}_{e})_{T}\quad\forall~T\in\mathcal{T}_{h},

where $v_{e}$ is the coefficient of $\bm{\Phi}_{e}$ such that $\bm{v}^{R}=\sum_{e\in\mathcal{E}^{0}}v_{e}\bm{\Phi}_{e}$ . Summation of the above inequality over all the elements $T\in\mathcal{T}_{h}$ gives the first inequality in 3.7 with $C_{*}=1/(d+1)$ .

Let us prove the second inequality of 3.7. For any basis $\bm{\Phi}_{e}\in\mathcal{RT}$ , $e\subset\partial T$ , since $\bm{\Phi}_{e}\cdot\bm{n}_{e^{\prime}}$ vanishes on any face $e^{\prime}\neq e$ , we have $\bm{\Phi}_{e}|_{e}\cdot\bm{n}_{e}=(|T|/|e|)\nabla\cdot\bm{\Phi}_{e}|_{T},$ where $|\cdot|$ denotes the measure and we use the divergence theorem here. By a series of simple calculations we arrive at

(3.10)

\|\bm{\Phi}_{e}\cdot\bm{n}_{e}\|_{e}^{2}=(|T|/|e|)\|\nabla\cdot\bm{\Phi}_{e}\|_{T}^{2}.

From [5, Eqs. 2.1.75&2.1.78] and scaling arguments one can obtain

(3.11)

||\bm{\Phi}_{e}||_{T}\leq Ch_{T}||\nabla\cdot\bm{\Phi}_{e}||_{T}.

Since we assume the partition is shape-regular, the combination of 3.10 and 3.11 gives that $||\bm{\Phi}_{e}||_{T}\leq Ch_{e}^{1/2}||\bm{\Phi}_{e}\cdot\bm{n}_{e}||_{e}$ . Then for the elements $T_{1}$ and $T_{2}$ which share the face $e$ , we have

v_{e}^{2}||\bm{\Phi}_{e}||_{T_{i}}^{2}\leq Ch_{e}v_{e}^{2}||\bm{\Phi}_{e}\cdot\bm{n}_{e}||_{e}^{2}=Ch_{e}||v_{e}\bm{\Phi}_{e}\cdot\bm{n}_{e}||_{e}^{2}=Ch_{e}||\bm{v}^{R}\cdot\bm{n}_{e}||_{e}^{2}\leq C||\bm{v}^{R}||_{T_{i}},

for $i=1,2$ , where in the last inequality we use the trace inequality. Sum over all the elements and apply the fact that the mesh is shape-regular, then the second inequality in 3.7 can be proven.

Finally, the two inequalities in 3.8 follow immediately from an inverse inequality and the summation of 3.11 over all the elements, respectively. Thus we complete the proof. ∎

3.2. Approximation

In this subsection, we will construct an interpolation operator $\Pi_{h}:V\cap[C^{0}(\bar{\Omega})]^{d}\rightarrow V_{h}$ which satisfies the following properties

(3.12)

\nabla\cdot\Pi_{h}\bm{v}=P_{h}\nabla\cdot\bm{v}

(3.13)

||\bm{v}-\Pi_{h}\bm{v}||_{T}+h_{T}|||\bm{v}-\Pi_{h}\bm{v}|||_{T}\leq{C_{I}h_{T}^{2}|\bm{v}|_{2,T},}

for all $T\in\mathcal{T}_{h},\bm{v}\in[H^{2}({\Omega})]^{d}$ . Here $|||\cdot|||_{T}$ is defined as

|||\bm{v}|||_{T}^{2}=a_{h}^{T}(\bm{v},\bm{v}),

where $a_{h}^{T}(\cdot,\cdot)$ is defined in the natural sense such that $a_{h}(\cdot,\cdot)=\sum_{T\in\mathcal{T}_{h}}a_{h}^{T}(\cdot,\cdot)$ .

To this end, denote by $\Pi_{h}^{1}:V\cap[C^{0}(\bar{\Omega})]^{d}\rightarrow V_{h}^{1}$ the usual nodal interpolation and we define $\Pi_{h}$ as

(3.14)

\Pi_{h}\bm{v}=\Pi_{h}^{1}\bm{v}+\Pi_{h}^{R}(\bm{v}-\Pi_{h}^{1}\bm{v})\quad\forall~\bm{v}\in V\cap[C^{0}(\bar{\Omega})]^{d}.

For the properties of $\Pi_{h}^{R}$ , see Lemma 2. We have the following lemma:

Lemma 4.

The interpolation operator $\Pi_{h}:V\cap[C^{0}(\bar{\Omega})]^{d}\rightarrow V_{h}$ defined in 3.14 satisfies the properties 3.12 and 3.13.

Proof.

Since $\nabla\cdot V_{h}^{1}\subset W_{h},$ we have

(3.15)

P_{h}\nabla\cdot\bm{v}_{h}^{1}=\nabla\cdot\bm{v}_{h}^{1}\quad\forall~\bm{v}_{h}^{1}\in V_{h}^{1}.

Eqs. 3.4 and 3.15 imply

\nabla\cdot\Pi_{h}\bm{v}=\nabla\cdot\Pi_{h}^{1}\bm{v}+\nabla\cdot\Pi_{h}^{R}(\bm{v}-\Pi_{h}^{1}\bm{v})=\nabla\cdot\Pi_{h}^{1}\bm{v}+P_{h}\nabla\cdot(\bm{v}-\Pi_{h}^{1}\bm{v})=P_{h}\nabla\cdot\bm{v},

Thus 3.12 has been proved.

Let us prove 3.13. For simplicity, we will only prove it for $a^{0}(\cdot,\cdot)$ . From the approximation theory we know

(3.16)

\|\bm{v}-\Pi_{h}^{1}\bm{v}\|_{T}+h_{T}|\bm{v}-\Pi_{h}^{1}\bm{v}|_{1,T}\leq{C_{I_{1}}h_{T}^{2}|\bm{v}|_{2,T}}\quad\forall~T\in\mathcal{T}_{h},

which, together with 3.2, gives

(3.17)

\|\Pi_{h}^{R}(\bm{v}-\Pi_{h}^{1}\bm{v})\|_{T}\leq{C_{R}C_{I_{1}}h_{T}^{2}|\bm{v}|_{2,T}.}

Then 3.13 follows immediately from 3.16 and 3.17, together with the triangle inequality. ∎

3.3. Stability and boundedness

From 3.6, the coercivity and the boundedness of $a_{h}(\cdot,\cdot)$ are trivial. We just state the following lemmas without proof.

Lemma 5 (Coercivity).

For any $\bm{v}\in V(h)$ we have

(3.18)

a_{h}(\bm{v},\bm{v})=|||\bm{v}|||^{2}.

Lemma 6 (Boundedness).

For any $\bm{w},\bm{v}\in V(h)$ we have

(3.19)

a_{h}(\bm{v},\bm{w})\leq|||\bm{v}||||||\bm{w}|||.

Next, let us prove the discrete inf-sup condition here. Similarly to the Bernardi and Raugel element, here we can not use the interpolation operator $\Pi_{h}$ to prove the inf-sup conditions since it is not well-defined on $[H^{1}(\Omega)]^{d}$ . However, using the similar technique in [4] and Lemma 4, we can construct an interpolation $\mathcal{\pi}_{h}:V\rightarrow V_{h}$ which satisfies

(3.20)

\nabla\cdot\mathcal{\pi}_{h}\bm{v}=P_{h}\nabla\cdot\bm{v}\quad\forall~\bm{v}\in V,

and

(3.21)

|||\mathcal{\pi}_{h}\bm{v}|||\leq C_{F}||\bm{v}||_{1}\quad\forall~\bm{v}\in V.

For more details, we refer the readers to [4] and [20, pp. 136-138]. Based on this interpolation, we will prove the following inf-sup conditions.

Lemma 7 (LBB Stability).

There exists a positive constant $\beta_{1},$ independent of $h,$ such that

(3.22)

\sup_{\bm{v}_{h}\in V_{h}\setminus\{\bm{0}\}}\frac{b(\bm{v}_{h},q_{h})}{|||\bm{v}_{h}|||}\geq\beta_{1}\|q_{h}\|\quad\forall~q_{h}\in W_{h}.

Proof.

We first use the operator $\mathcal{\pi}_{h}$ to obtain

(3.23)

\sup_{\bm{v}_{h}\in V_{h}\setminus\{\bm{0}\}}\frac{b(\bm{v}_{h},q_{h})}{|||\bm{v}_{h}|||}\geq\sup_{\bm{v}\in V\setminus\{\bm{0}\}}\frac{b\left(\mathcal{\pi}_{h}\bm{v},q_{h}\right)}{|||\mathcal{\pi}_{h}\bm{v}|||}=\sup_{\bm{v}\in V\setminus\{\bm{0}\}}\frac{b(\bm{v},q_{h})}{|||\mathcal{\pi}_{h}\bm{v}|||}.

Thus, substituting 3.21 into the inequality 3.23 gives

\sup_{\bm{v}_{h}\in V_{h}\setminus\{\bm{0}\}}\frac{b(\bm{v}_{h},q_{h})}{|||\bm{v}_{h}|||}\geq C_{F}^{-1}\sup_{\bm{v}\in V\setminus\{\bm{0}\}}\frac{b(\bm{v},q_{h})}{\|\bm{v}\|_{1}}\geq\beta_{1}\|q_{h}\|,

where we have used the inf-sup condition for the continuous case 2.3. ∎

3.4. Consistency Error

The continuous solutions $(\bm{u},p)$ satisfy

(3.24)

\begin{split}\nu\left(-\Delta\bm{u},\bm{v}\right)-b\left(\bm{v},p\right)&=\left(\bm{f},\bm{v}\right)\quad\forall~\bm{v}\in V(h),\\ b\left(\bm{u},q\right)&=0\quad\forall~q\in W,\end{split}

which can be rewritten as

(3.25)

\begin{split}\nu a_{h}(\bm{u},\bm{v})-b(\bm{v},p)&=(\bm{f},\bm{v})-\nu\left((-\Delta\bm{u},\bm{v})-a_{h}(\bm{u},\bm{v})\right)\quad\forall~\bm{v}\in V(h),\\ b(\bm{u},q)&=0\quad\forall~q\in W.\end{split}

Then $\left(-\Delta\bm{u},\bm{v})-a_{h}(\bm{u},\bm{v}\right)$ is the consistency error of our method. From the arguments in [48, Sect. 3], one can see, in general, $\left|(-\Delta\bm{u},\bm{v})-a_{h}\left(\bm{u},\bm{v}\right)\right|\neq 0.$

Lemma 8 (Consistency).

The following inequality holds:

(3.26)

\left|(-\Delta\bm{u},\bm{v})-a_{h}\left(\bm{u},\bm{v}\right)\right|\leq(\sum_{T\in\mathcal{T}_{h}}h_{T}^{2}\|\Delta\bm{u}\|_{T}^{2})^{1/2}(\sum_{T\in\mathcal{T}_{h}}h_{T}^{-2}\left\|\bm{v}^{R}\right\|_{T}^{2})^{1/2}\quad\forall~\bm{v}\in{{V}(h)}.

Proof.

Note that $a_{h}(\bm{u},\bm{v}_{h})=a(\bm{u},\bm{v}^{1})=(-\Delta\bm{u},\bm{v}^{1})$ . Thus, we have


(3.27a)	$\displaystyle\left\|(-\Delta\bm{u},\bm{v})-a_{h}\left(\bm{u},\bm{v}\right)\right\|$	$\displaystyle=\left\|(-\Delta\bm{u},\bm{v}^{R})\right\|\leq\sum_{T\in\mathcal{T}_{h}}h_{T}\\|\Delta\bm{u}\\|_{T}h_{T}^{-1}\left\\|\bm{v}^{R}\right\\|_{T}$
(3.27b)		$\displaystyle\leq(\sum_{T\in\mathcal{T}_{h}}h_{T}^{2}\\|\Delta\bm{u}\\|_{T}^{2})^{1/2}(\sum_{T\in\mathcal{T}_{h}}h_{T}^{-2}\left\\|\bm{v}^{R}\right\\|_{T}^{2})^{1/2}.$

∎

3.5. Error estimates

Based on the above results, the error analysis can be done now. In addition, the following error equations are needed:


(3.28a)	$\displaystyle\nu a_{h}\left(\bm{u}-\bm{u}_{h},\bm{v}_{h}\right)-b\left(\bm{v}_{h},p-p_{h}\right)$	$\displaystyle=-\nu\delta_{h}(\bm{u},\bm{v}_{h})\quad\forall~\bm{v}_{h}\in V_{h},$
(3.28b)	$\displaystyle b\left(\bm{u}-\bm{u}_{h},q_{h}\right)$	$\displaystyle=0\quad\forall~q_{h}\in W_{h},$

where $\delta_{h}:V(h)\times V(h)\rightarrow\mathbb{R}$ is the consistency error defined as

\delta_{h}(\bm{u},\bm{v})=(-\Delta\bm{u},\bm{v})-a_{h}\left(\bm{u},\bm{v}\right).

These error equations can be obtained by subtracting 2.6 from 3.25. In what follows we shall use $C_{\alpha}$ to denote the constant in the upper bound of the consistency error such that

|\delta_{h}(\bm{u},\bm{v})|\leq C_{\alpha}h||\Delta\bm{u}||\left\{a^{R}(\bm{v}^{R},\bm{v}^{R})\right\}^{1/2}.

It follows from 3.26 that $C_{\alpha}$ only depends on the parameters $\alpha_{T}$ when $a^{R}=a^{0}$ . Now we are in a position to present an error estimate for the finite element scheme 2.6. Our attention is mainly paid to the velocity estimates, which are independent of the pressure and the viscosity.

Theorem 3.1.

Let $(\bm{u},p)$ be the solution of 2.1 and $\left(\bm{u}_{h},p_{h}\right)\in V_{h}\times W_{h}$ be the solution of 2.6. Assume that $(\bm{u},p)\in[H^{2}(\Omega)]^{d}\times H^{1}(\Omega)$ . Then there exists a constant $C$ independent of $h$ , $\nu$ and $\beta_{1}$ such that

(3.29)

|||\bm{u}-\bm{u}_{h}|||\leq Ch|\bm{u}|_{2},\quad\|p_{h}-P_{h}p\|\leq\frac{\nu}{\beta_{1}}Ch|\bm{u}|_{2},

and

(3.30)

\left\|p-p_{h}\right\|\leq Ch\left(\frac{\nu}{\beta_{1}}|\bm{u}|_{2}+\left|p\right|_{1}\right).

Here $\beta_{1}$ is the coefficient of the discrete inf-sup condition.

Proof.

Introduce the notation:

\xi_{h}=\bm{u}_{h}-\Pi_{h}\bm{u},\quad\eta_{h}=p_{h}-P_{h}p,\quad\xi=\bm{u}-\Pi_{h}\bm{u},\quad\eta=p-P_{h}p.

We note that 3.12 implies $\nabla\cdot\Pi_{h}\bm{u}=0$ , which, together with 2.12 gives that

(3.31)

\nabla\cdot\xi_{h}=0.

It follows from the error equations 3.28 that


(3.32a)	$\displaystyle\nu a_{h}\left(\xi_{h},\bm{v}_{h}\right)-b\left(\bm{v}_{h},\eta_{h}\right)$	$\displaystyle=\nu a_{h}(\xi,\bm{v}_{h})-b(\bm{v}_{h},\eta)+\nu\delta_{h}(\bm{u},\bm{v}_{h})\quad\forall~\bm{v}_{h}\in V_{h},$
(3.32b)	$\displaystyle b\left(\xi_{h},q_{h}\right)$	$\displaystyle=b(\xi,q_{h})=0\quad\forall~q_{h}\in W_{h}.$

By taking $\bm{v}_{h}=\xi_{h}$ in 3.32a and $q_{h}=\eta_{h}$ in 3.32b, the sum of 3.32a and 3.32b together with 3.31 gives

(3.33)

a_{h}\left(\xi_{h},\xi_{h}\right)=a_{h}\left(\xi,\xi_{h}\right)+\delta_{h}(\bm{u},\xi_{h}),

where $\nu$ has been eliminated. Thus, it follows from the coercivity 3.18, the boundedness 3.19 and the consistency error 3.26 that

|||\xi_{h}|||^{2}\leq|||\xi|||\times|||\xi_{h}|||+C_{\alpha}h|\bm{u}|_{2}|||\xi_{h}|||,

which implies

|||\xi_{h}|||\leq|||\xi|||+C_{\alpha}h|\bm{u}|_{2}.

The above estimate is equivalent to

|||\bm{u}_{h}-\Pi_{h}\bm{u}|||\leq|||\bm{u}-\Pi_{h}\bm{u}|||+C_{\alpha}h|\bm{u}|_{2}.

By taking $s=2$ in 3.13, a combination of the triangle inequality and the above error estimate gives

(3.34)

|||\bm{u}-\bm{u}_{h}|||\leq 2|||\bm{u}-\Pi_{h}\bm{u}|||+C_{\alpha}h|\bm{u}|_{2}\leq\left(2C_{I}+C_{\alpha}\right)h|\bm{u}|_{2}.

Let us estimate the pressure error. Since $\nabla\cdot V_{h}=W_{h}$ , it follows from 3.1 that $b(\bm{v}_{h},p-P_{h}p)=0$ for all $\bm{v}_{h}\in V_{h},$ which, together with the discrete inf-sup condition 3.22 gives

	$\displaystyle\left\\|p_{h}-P_{h}p\right\\|$	$\displaystyle\leq\frac{1}{\beta_{1}}\sup_{\bm{v}_{h}\in V_{h}\setminus\{\bm{0}\}}\frac{b\left(\bm{v}_{h},p_{h}-P_{h}p\right)}{\|\|\|\bm{v}_{h}\|\|\|}=\frac{1}{\beta_{1}}\sup_{\bm{v}_{h}\in V_{h}\setminus\{\bm{0}\}}\frac{b\left(\bm{v}_{h},p_{h}-p\right)}{\|\|\|\bm{v}_{h}\|\|\|}$
		$\displaystyle=\frac{1}{\beta_{1}}\sup_{\bm{v}_{h}\in V_{h}\setminus\{\bm{0}\}}\frac{-\nu a_{h}\left(\bm{u}-\bm{u}_{h},\bm{v}_{h}\right)-\nu\delta_{h}(\bm{u},\bm{v}_{h})}{\|\|\|\bm{v}_{h}\|\|\|}$
		$\displaystyle\leq\frac{\nu}{\beta_{1}}\sup_{\bm{v}_{h}\in V_{h}\setminus\{\bm{0}\}}\frac{1}{\|\|\|\bm{v}_{h}\|\|\|}\|\|\|\bm{v}_{h}\|\|\|\left(\|\|\|\bm{u}-\bm{u}_{h}\|\|\|+C_{\alpha}h\|\bm{u}\|_{2}\right)\leq\frac{\nu}{\beta_{1}}\left(\|\|\|\bm{u}-\bm{u}_{h}\|\|\|+C_{\alpha}h\|\bm{u}\|_{2}\right).$

Now using the above estimate and the estimate 3.34 we get

\left\|p_{h}-P_{h}p\right\|\leq 2\frac{\nu}{\beta_{1}}\left(C_{I}+C_{\alpha}\right)h|\bm{u}|_{2}.

Thus we complete the proof of 3.29. Then the error estimate 3.30 follows from the standard triangle inequality, the above estimate and taking $m=1$ in 3.3. ∎

Remark 5.

Classical mixed methods which do not satisfy divergence-free (or pressure-robust) property can only obtain an abstract estimate as

\left\|\nabla\left(\bm{u}-\bm{u}_{h}^{cla}\right)\right\|_{L^{2}(\Omega)}\leq C\left(\inf_{{\bm{v}}_{h}\in V_{h}^{cla}}\left\|\nabla\left(\bm{u}-{\bm{v}}_{h}\right)\right\|_{L^{2}(\Omega)}+{\nu}^{-1}\inf_{q_{h}\in W_{h}^{cla}}\left\|p-q_{h}\right\|_{L^{2}(\Omega)}\right).

See [25, eq. 3.5]. In classical mixed methods, when $\nu$ is small, or when the best approximation error of the pressure is large, the error bound of the velocity may be large.

To derive an optimal order $L^{2}$ -error estimate for the velocity approximation, we seek $(\bm{z},\lambda)$ satisfying

-\Delta\bm{z}+\nabla\lambda=\bm{u}-\bm{u}_{h},\ \nabla\cdot\bm{z}=0\text{ in }\Omega,\quad\bm{z}=\bm{0}\text{ on }\partial\Omega,

which is also a Stokes system with the viscosity equal to 1.

Here we assume that $\Omega$ is convex, and since $\bm{u}-\bm{u}_{h}\in\left[L^{2}\left(\Omega\right)\right]^{d}$ , we have $(\bm{z},\lambda)\in[H^{2}(\Omega)]^{d}\times H^{1}(\Omega)$ and the following regularity estimate holds true (see [14, 19]):

(3.35)

\|\bm{z}\|_{2}+\|\lambda\|_{1}\leq C\left\|\bm{u}-\bm{u}_{h}\right\|.

Define $\delta_{h}^{*}:V(h)\times V(h)\rightarrow\mathbb{R}$ as $\delta_{h}^{*}(\bm{v},\bm{z})=\delta_{h}(\bm{z},\bm{v}),$ which is the adjoint consistency error [3]. In fact, since $-\Delta$ is a self-adjoint operator and $a_{h}(\cdot,\cdot)$ is symmetric, the properties of consistency and adjoint consistency are equivalent.

Similarly to 3.25, note that $a_{h}(\cdot,\cdot)$ is symmetric, then we have

(3.36)

a_{h}(\bm{v},\bm{z})-b(\bm{v},\lambda)=\left(\bm{u}-\bm{u}_{h},\bm{v}\right)-\delta_{h}^{*}(\bm{v},\bm{z})\quad\forall~\bm{v}\in V(h).

Theorem 3.2.

(3.37)

\left\|\bm{u}-\bm{u}_{h}\right\|\leq Ch^{2}|\bm{u}|_{2}.

Proof.

By taking $\bm{v}=\bm{u}-\bm{u}_{h}$ in 3.36 we get

(3.38)

a_{h}\left(\bm{u}-\bm{u}_{h},\bm{z}\right)-b\left(\bm{u}-\bm{u}_{h},\lambda\right)=\left\|\bm{u}-\bm{u}_{h}\right\|^{2}-\delta_{h}^{*}(\bm{u}-\bm{u}_{h},\bm{z}).

Eq. 3.28b implies that

(3.39)

b\left(\bm{u}-\bm{u}_{h},\lambda\right)=b\left(\bm{u}-\bm{u}_{h},\lambda-P_{h}\lambda\right).

Meanwhile, Eq. 3.28a and the fact $\nabla\cdot\Pi_{h}\bm{z}=0$ give

(3.40)	$\displaystyle a_{h}\left(\bm{u}-\bm{u}_{h},\bm{z}\right)$	$\displaystyle=a_{h}\left(\bm{u}-\bm{u}_{h},\bm{z}-\Pi_{h}\bm{z}\right)+a_{h}\left(\bm{u}-\bm{u}_{h},\Pi_{h}\bm{z}\right)$
		$\displaystyle=a_{h}\left(\bm{u}-\bm{u}_{h},\bm{z}-\Pi_{h}\bm{z}\right)+\frac{1}{\nu}b\left(\Pi_{h}\bm{z},p-p_{h}\right)-\delta_{h}(\bm{u},{\Pi_{h}\bm{z}})$
		$\displaystyle=a_{h}\left(\bm{u}-\bm{u}_{h},\bm{z}-\Pi_{h}\bm{z}\right)-\delta_{h}(\bm{u},{\Pi_{h}\bm{z}}).$

Substituting 3.39 and 3.40 into 3.38 one can obtain

\left\|\bm{u}-\bm{u}_{h}\right\|^{2}=a_{h}\left(\bm{u}-\bm{u}_{h},\bm{z}-\Pi_{h}\bm{z}\right)-b\left(\bm{u}-\bm{u}_{h},\lambda-P_{h}\lambda\right)+\delta_{h}^{*}(\bm{u}-\bm{u}_{h},\bm{z})-\delta_{h}(\bm{u},{\Pi_{h}\bm{z}}).

From 3.14 and 3.26 we have

|\delta_{h}(\bm{u},\Pi_{h}\bm{z})|\leq(\sum_{T\in\mathcal{T}_{h}}h_{T}^{2}\|\Delta\bm{u}\|_{T}^{2})^{1/2}(\sum_{T\in\mathcal{T}_{h}}h_{T}^{-2}\left\|\Pi_{h}^{R}(\bm{z}-\Pi_{h}^{1}\bm{z})\right\|_{T}^{2})^{1/2}\leq Ch^{2}|\bm{u}|_{2}||\bm{z}||_{2}\quad({\rm see~\lx@cref{creftype~refnum}{rt0part}}),

and $|\delta_{h}^{*}(\bm{u}-\bm{u}_{h},\bm{z})|\leq C_{\alpha}h||\bm{z}||_{2}|||\bm{u}-\bm{u}_{h}|||,$ which imply

\left\|\bm{u}-\bm{u}_{h}\right\|^{2}\leq C|||\bm{u}-\bm{u}_{h}|||\left(|||\bm{z}-\Pi_{h}\bm{z}|||+\left\|\lambda-P_{h}\lambda\right\|\right)+Ch(h|\bm{u}|_{2}+|||\bm{u}-\bm{u}_{h}|||)||\bm{z}||_{2}.

It follows from the above estimate, the interpolation error estimates 3.3 and 3.13, and the regularity estimate 3.35 that

\left\|\bm{u}-\bm{u}_{h}\right\|^{2}\leq Ch\left(|||\bm{u}-\bm{u}_{h}|||+h|\bm{u}|_{2}\right)\left\|\bm{u}-\bm{u}_{h}\right\|,

which implies that

\left\|\bm{u}-\bm{u}_{h}\right\|\leq Ch\left(|||\bm{u}-\bm{u}_{h}|||+h|\bm{u}|_{2}\right).

Then the estimate 3.37 follows immediately from the above inequality and the estimate 3.29. Thus we complete the proof. ∎

Remark 6.

From Theorem 3.1 and 3.6 we have $\|\nabla(\bm{u}-\bm{u}_{h}^{1})\|\leq Ch|u|_{2}$ . Thus, $\bm{u}_{h}^{1}$ is a conforming approximation of the continuous solution $\bm{u}$ with optimal order in $H^{1}$ norm. In addition, Theorem 3.1 also imply that $||\bm{u}_{h}^{R}||\leq Ch^{2}|\bm{u}|_{2}$ , which, together with Theorem 3.2 and a triangle inequality, demonstrates that $\bm{u}_{h}^{1}$ is also convergent to $\bm{u}$ with the optimal order in $L^{2}$ norm.

We summarize the above remark in the following corollary.

Corollary 1.

\left\|\bm{u}-\bm{u}_{h}^{1}\right\|+h\left\|\nabla(\bm{u}-\bm{u}_{h}^{1})\right\|\leq Ch^{2}|\bm{u}|_{2}.

Remark 7.

Although our method is stable and optimally convergent as long as the parameters $\alpha_{T},T\in\mathcal{T}_{h}$ are positive, here we can not give a set of optimal choices in the sense of the accuracy. This is still an open and interesting topic. Consider the simplest case. Assuming that $\alpha_{T}\equiv\alpha$ and taking $a^{R}=a^{0}$ may give some rough answers for this topic. Note that 3.33, 3.26 and the definition of $|||\cdot|||$ imply

|||\xi_{h}|||\leq\left\{||\nabla(\bm{u}-\Pi_{h}^{1}\bm{u})||^{2}+\alpha\sum_{T\in\mathcal{T}_{h}}h_{T}^{-2}||\Pi_{h}^{R}(\bm{u}-\Pi_{h}^{1}\bm{u})||_{T}^{2}\right\}^{1/2}+\left\{\alpha^{-1}\sum_{T\in\mathcal{T}_{h}}h_{T}^{2}||\Delta\bm{u}||_{T}^{2}\right\}^{1/2}.

Since $||\Pi_{h}^{R}(\bm{u}-\Pi_{h}^{1}\bm{u})||_{T}\leq C_{R}C_{I_{1}}h_{T}^{2}|\bm{u}|_{2,T}$ (see 3.17) and $||\Delta\bm{u}||_{T}^{2}\leq d|\bm{u}|_{2,T}^{2}$ , by a series of simple calculations we arrive at

|||\xi_{h}|||^{2}\leq 2||\nabla(\bm{u}-\Pi_{h}^{1}\bm{u})||^{2}+2(\alpha C_{R}^{2}C_{I_{1}}^{2}+\alpha^{-1}d)\sum_{T\in\mathcal{T}_{h}}h_{T}^{2}|\bm{u}|_{2,T}^{2}.

The above inequality implies that taking $\alpha=\sqrt{d}/(C_{R}C_{I_{1}})$ might be a good choice. This requires a severe estimate of $C_{R}$ and $C_{I_{1}}$ , which depends on the shape regularity of the mesh. In the numerical experiments we will study the performance of our method with different parameters on various meshes.

Remark 8.

In Remark 3 we point out that in practice we prefer the bilinear form $a^{\operatorname{div}}(\cdot,\cdot)$ . However, in the numerical experiments we find that this bilinear form is a little sensitive to the change of the parameters. To reduce the sensitivity, some strategies could be applied. For example, we can modify the bilinear form $a_{h}(\cdot,\cdot)$ to

(3.41)

a_{h}^{m}(\bm{u},\bm{v})=a_{h}(\bm{u},\bm{v})-\sum_{T\in\mathcal{T}_{h}}\gamma_{T}(\nabla\cdot\bm{u}^{1},\nabla\cdot\bm{v}^{1})_{T},

where $0\leq\gamma_{T}<\alpha_{T}/(d+1),T\in\mathcal{T}_{h}$ are the parameters. The analysis of $a_{h}^{m}(\cdot,\cdot)$ is very similar to $a_{h}(\cdot,\cdot)$ . The main difference is that we can not prove the coercivity of $a_{h}^{m}(\cdot,\cdot)$ on the whole $V(h)$ . However, we can prove that $a_{h}^{m}(\cdot,\cdot)$ is coercive on the subspace of $V(h)$

Z(h)=\{\bm{v}\in V(h):b(\bm{v},q)=0\quad\forall\ q\in W\}=\{\bm{v}\in V(h):\nabla\cdot\bm{v}=0\}.

Indeed, for any $\bm{v}\in Z(h)$ , we have $\nabla\cdot\bm{v}^{1}=-\nabla\cdot\bm{v}^{R}$ , which implies that


(3.42a)	$\displaystyle\gamma_{T}(\nabla\cdot\bm{v}^{1},\nabla\cdot\bm{v}^{1})_{T}$	$\displaystyle=\gamma_{T}(\nabla\cdot\bm{v}^{R},\nabla\cdot\bm{v}^{R})_{T}\leq(d+1)\gamma_{T}\sum_{e\in\partial T\cap\mathcal{E}^{0}}v_{e}^{2}(\nabla\cdot\bm{\Phi}_{e},\nabla\cdot\bm{\Phi}_{e})_{T}$
(3.42b)		$\displaystyle<\alpha_{T}\sum_{e\in\partial T\cap\mathcal{E}^{0}}v_{e}^{2}(\nabla\cdot\bm{\Phi}_{e},\nabla\cdot\bm{\Phi}_{e})_{T}.$

Then the coercivity of $a_{h}^{m}(\cdot,\cdot)$ on $Z(h)$ follows immediately from the above inequality. Note that $Z(h)$ include the finite element space $Z_{h}=\{\bm{v}_{h}\in V_{h}:b(\bm{v}_{h},q_{h})=0\quad\forall\ q_{h}\in W_{h}\}=\{\bm{v}_{h}\in V_{h}:\nabla\cdot\bm{v}_{h}=0\}.$ In addition, the boundedness of the bilinear forms and the inf-sup conditions could be similarly obtained. Then from the classical theory on saddle point systems [6], the existence and uniqueness of the solutions are guaranteed. The error estimate also follows similarly. The numerical experiments below (see section 4.2) demonstrate that this modified scheme does reduce the sensitivity.

4. Numerical experiments

We perform some numerical experiments on four examples. The examples in section 4.1 and section 4.3 are from [25]. In some examples, we also give some results computed by the classical low order Bernardi and Raugel pair [4] to make a comparison.

4.1. Convergence test

The exact solutions are prescribed as follows,

\bm{u}=200\left(x^{2}(1-x)^{2}y(1-y)(1-2y),-x(1-x)(1-2x)y^{2}(1-y)^{2}\right)^{\top},

and

p=10\left(\left(x-1/2\right)^{3}y^{2}+(1-x)^{3}\left(y-1/2\right)^{3}\right).

The velocity field has the form of a large vortex [25]. We set $\Omega=(0,1)\times(0,1)$ and consider a small viscosity case, $\nu=10^{-6}$ , to show the robustness of our methods. The corresponding results are shown in Table 1, where we also give some results from the Bernardi and Raugel element. For the choice of parameters, we take $\alpha_{T}\equiv 20,T\in\mathcal{T}_{h}$ for $a^{R}=a^{0}$ and $a^{R}=a^{D}$ . In the case that $a^{R}=a^{\operatorname{div}}$ , we take $\alpha_{T}\equiv 1.5$ . One could see, the velocity errors of our methods (“compact H(div) method”) are roughly $10^{4}$ times smaller than the Bernardi and Raugel method. Moreover, comparing with the bottom line of Table 1, one can find that the discrete pressures obtained by our methods are almost equal to the best approximation of the continuous pressure in $L^{2}(\Omega)$ .

Table 1. section 4.1. The

L^{2}

and

H^{1}

(seminorm) velocity errors and

L^{2}

pressure errors. CHdiva0: the method 2.6 with

a^{0}

2.8; CHdivaD: the method 2.6 with

a^{D}

2.9; CHdivadiv: the method 2.6 with

a^{\operatorname{div}}

2.10; B-R: the Bernardi and Raugel finite element method [4]. The

L^{2}

projection errors of the pressure are listed at the bottom.

		$h=0.1$	$0.05$	$0.025$	$0.0125$	$0.00625$
$\left\\|\bm{u}-\bm{u}_{h}\right\\|$	${\rm CHdiva0}$	1.51E-2	3.87E-3	9.73E-4	2.43E-4	6.09E-5
	${\rm CHdivaD}$	1.60E-2	4.11E-3	1.03E-3	2.60E-4	6.52E-5
	${\rm CHdivadiv}$	1.97E-2	5.05E-3	1.27E-3	3.19E-4	8.16E-5
	${\rm B-R}$	130.45	33.80	8.55	2.14	5.37E-1
$\left\\|\nabla(\bm{u}-\bm{u}_{h}^{1})\right\\|$	${\rm CHdiva0}$	1.12	5.67E-1	2.84E-1	1.42E-1	7.11E-2
	${\rm CHdivaD}$	1.13	5.68E-1	2.84E-1	1.42E-1	7.11E-2
	${\rm CHdivadiv}$	1.13	5.69E-1	2.84E-1	1.42E-1	7.11E-2
	${\rm B-R}$	7257.64	3895.73	2010.36	1020.13	513.69
$\left\\|p-p_{h}\right\\|$	${\rm CHdiva0}$	2.52E-2	1.26E-2	6.32E-3	3.16E-3	1.58E-3
	${\rm CHdivaD}$	2.52E-2	1.26E-2	6.32E-3	3.16E-3	1.58E-3
	${\rm CHdivadiv}$	2.52E-2	1.26E-2	6.32E-3	3.16E-3	1.58E-3
	${\rm B-R}$	2.61E-2	1.31E-2	6.55E-3	3.27E-3	1.63E-3
$\left\\|p-P_{h}p\right\\|$		2.52E-2	1.26E-2	6.32E-3	3.16E-3	1.58E-3

4.2. Effects of different values of the parameters $\alpha_{T}$

In this example we focus on the case that $a^{R}=a^{\operatorname{div}}$ and investigate the effects of different choices of the parameters on four meshes (see Fig. 3). The exact solutions are the same as section 4.1. We set $\alpha_{T}\equiv\alpha$ with $\alpha=0.1:0.1:10$ , that is, we test 100 different values of the parameters which are from 0.1 to 10. The results from the original scheme 2.6 and the modified scheme 3.41 are shown in Fig. 4 and Fig. 5, respectively. It is demonstrated that the most stable quantity is $||\bm{u}-\bm{u}_{h}^{1}||$ ( $L^{2}u^{1}$ ), which is very close to the nodal interpolation error ( $L^{2}\Pi^{1}u$ ) in the whole region. The quantities $||\bm{u}-\bm{u}_{h}||$ ( $L^{2}u$ ) and $||\nabla(\bm{u}-\bm{u}_{h}^{1})||$ ( $H^{1}u^{1}$ ) reduce violently before $\alpha<1$ . All the minimum points of $\|\bm{u}-\bm{u}_{h}\|$ are in $(1,5)$ . The comparison between the results in Fig. 4 and Fig. 5 shows that the modified scheme 3.41 does reduce the sensitivity to the change of $\alpha$ when $\alpha>1$ .

Based on the above observations, we recommend that one choose $\alpha>1$ in practice.

4.3. Stokes flow with Coriolis force

In this example, we consider the flows with strong Coriolis forces which appear in several real-world applications such as meteorology. The simplest model is described as follows (see [34, 25]),

(4.1)

-\nu\Delta\bm{u}+\nabla p+2\bm{\omega}\times\bm{u}=\bm{f},\quad\nabla\cdot\bm{u}=0\quad\text{in}\ \Omega,

where $\bm{\omega}$ is a constant angular velocity vector.

Here $\bm{u}=(u_{1},u_{2})^{T}$ . We embed it into three-dimensional space as $\bm{u}=(u_{1},u_{2},0)^{T}$ , and similarly for other variables. We take $\bm{\omega}=(0,0,\omega)^{T}$ then we get a two-dimensional model since $\bm{\omega}\times\bm{u}=\omega(-u_{2},u_{1},0)^{T}.$ The computational domain $\Omega$ is taken as a L-shaped domain: $(0,4)\times(0,2)\setminus[2,4)\times(0,1]$ . We choose $\nu=1$ . And the value of $\omega$ changes from $1$ , $100$ to $10000$ . We note that changing the magnitude $\omega$ of the Coriolis force will change only the pressure, and not the velocity (see [25, Example 6.5]). The inlet is set at $x=0$ and outlet is set at $x=4$ . Pure Dirichlet boundary condition is considered, where the volume preserving parabolic inflow and outflow profiles are given by

\bm{u}_{in}=\left(y(2-y)/2,0\right)^{\top},\quad\bm{u}_{out}=\left(4(2-y)(y-1),0\right)^{\top},

and homogeneous conditions are used in other parts of the boundary. Here the Coriolis term is discretized by $(2\bm{\omega}\times\bm{u}_{h},\bm{v}_{h})$ , and the nonhomogeneous boundary condition is imposed by enforcing $\bm{u}_{h}|_{\Gamma}=\Pi_{h}\bm{u}|_{\Gamma}$ . Of course, we do not know the exact solution $\bm{u}$ but the boundary value of $\Pi_{h}\bm{u}$ can be obtained by the known boundary value of $\bm{u}$ . The computational results are shown in Fig. 6.

One could see, our method is robust with respect to the magnitude $\omega$ of the Coriolis force, while undesired instabilities occur for the Bernardi and Raugel element when a strong Coriolis force exists ( $\omega=100,10000$ ).

4.4. The robustness of the pressure-robust ${P_{1}^{c}-P0}$ scheme

In this example we want to show that the construction method of the ${\bm{P}_{1}^{c}-P0}$ discretization on the Stokes equations could be extended to general incompressible fluid flows. We will use the Brinkman equation and the Coriolis force problem (like section 4.3 but the exact velocity is known) as two examples. The (scaled) Brinkman equation is described as

(4.2)

-\nu\Delta\bm{u}+\bm{u}+\nabla p=\bm{f},\quad\nabla\cdot\bm{u}=0\quad\text{in}~\Omega,

and the Coriolis force problem is described as 4.1. In two problems the exact velocity are the same as section 4.1 and we apply the homogenous velocity boundary conditions. In 4.2 we choose $\nu=10^{-6}$ and the pressure is also the same as section 4.1. In the Coriolis problem we choose $\nu=1$ and $\omega=10^{4}$ with $\bm{f}\equiv 0$ (resulting in a very large pressure gradient that balances the Coriolis force). The corresponding discretizations respectively read: (in $a_{h}(\cdot,\cdot)$ we choose $a^{R}=a^{\operatorname{div}}$ with $\alpha_{T}\equiv 1.5$ )


	$\displaystyle{\rm Find}~(\bm{u}_{h},p_{h})\in V_{h}\times W_{h}~{\rm such~that}~~~~~~~~$
(4.3a)		$\displaystyle\nu a_{h}(\bm{u}_{h},\bm{v}_{h})+a^{B}(\bm{u}_{h},\bm{v}_{h})-b(\bm{v}_{h},p_{h})$	$\displaystyle+b(\bm{u}_{h},q_{h})=(\bm{f},\bm{v}_{h})\quad\forall~\bm{v}_{h}\in V_{h},~q_{h}\in W_{h},$

and


	$\displaystyle{\rm Find}~(\bm{u}_{h},p_{h})\in V_{h}\times W_{h}~{\rm such~that}~~~~~~~~$
(4.4a)		$\displaystyle\nu a_{h}(\bm{u}_{h},\bm{v}_{h})+a^{C}(\bm{u}_{h},\bm{v}_{h})-b(\bm{v}_{h},p_{h})$	$\displaystyle+b(\bm{u}_{h},q_{h})=(\bm{f},\bm{v}_{h})\quad\forall~\bm{v}_{h}\in V_{h},~q_{h}\in W_{h},$

where

a^{B}(\bm{u}_{h},\bm{v}_{h})=(\bm{u}_{h}^{1},\bm{v}_{h}^{1})+(\bm{u}_{h}^{R},\bm{v}_{h}^{1})+(\bm{u}_{h}^{1},\bm{v}_{h}^{R})+(d+1)\sum_{T\in\mathcal{T}_{h}}\sum_{e\in\partial T\cap\mathcal{E}^{0}}u_{e}v_{e}(\bm{\Phi}_{e},\bm{\Phi}_{e})_{T},

and

a^{C}(\bm{u}_{h},\bm{v}_{h})=(2\bm{\omega}\times\bm{u}_{h}^{1},\bm{v}_{h}^{1})+(2\bm{\omega}\times\bm{u}_{h}^{R},\bm{v}_{h}^{1})+(2\bm{\omega}\times\bm{u}_{h}^{1},\bm{v}_{h}^{R}).

The final term in $a^{B}$ is to preserve the stability since the second term and the third term in $a^{B}$ are two indefinite terms. A similar term does not exist in $a^{C}$ since one could verify that $a^{C}$ itself has a skew-symmetric structure, that is, $a^{C}(\bm{v}_{h},\bm{v}_{h})=0$ for all $\bm{v}_{h}\in V_{h}$ . This time the $RT0$ - $RT0$ regions of the coefficient matrices are both diagonal for the Brinkman and Coriolis force equations. Thus, we can easily obtain the corresponding ${\bm{P}_{1}^{c}-P0}$ schemes via static condensation. Unlike the scheme for the Stokes equations, in these two cases, the $\bm{u}_{h}^{R}$ parts are expressed in terms of $\bm{u}_{h}^{1}$ and $p_{h}$ together, not only $p_{h}$ . Some convergence results are shown in Table 2 and Table 3.

Table 2. section 4.4. Errors for the Brinkman equation with

\nu=10^{-6}

h	$\|\|\bm{u}-\bm{u}_{h}\|\|$	order	$\|\|\nabla(\bm{u}-\bm{u}_{h}^{1})\|\|$	order	$\|\|p-p_{h}\|\|$	order
1E-1	1.23E-2		1.15		2.51E-2
5E-2	2.93E-3	2.07	5.73E-1	1.01	1.26E-2	0.99
2.5E-2	7.10E-4	2.04	2.85E-1	1.00	6.32E-3	0.99
1.25E-2	1.74E-4	2.02	1.42E-1	1.00	3.16E-3	0.99
6.25E-3	4.32E-5	2.01	7.11E-2	1.00	1.58E-3	0.99

Table 3. section 4.4. Errors for the Stokes equation with strong Coriolis force (

\omega=10^{4}

h	$\|\|\bm{u}-\bm{u}_{h}\|\|$	order	$\|\|\nabla(\bm{u}-\bm{u}_{h}^{1})\|\|$	order
1E-1	4.27E-2		1.18
5E-2	8.41E-3	2.34	5.71E-1	1.04
2.5E-2	1.37E-3	2.61	2.84E-1	1.00
1.25E-2	3.19E-4	2.10	1.42E-1	1.00
6.25E-3	8.11E-5	1.97	7.11E-2	1.00

5. Conclusions

We have developed a class of divergence-free mixed methods for the Stokes problem based on ${\bm{P}_{1}^{c}}\oplus RT0-P0$ pair. By putting this pair into a modified conforming-like formulation, our methods combine some advantages of conforming methods and $H(\operatorname{div})$ -conforming methods. In addition, a stabilized and pressure-robust ${\bm{P}_{1}^{c}}-P0$ discretization is proposed in this paper which is computationally cheap. Future work includes:

•

Extending our methods to the mixed velocity and stress boundary value problems which are also an important class of problems in real-world applications. We note that, in this extension, a minor modification to our methods is required to preserve the convergence speed.
•

Extending our methods to the Navier-Stokes equations (NSEs). The key issue is the discretization of the nonlinear term in NSE. In this regard, we will propose a modified convective formulation which conserves the energy, linear momentum and angular momentum in some cases. This is based on the construction of our velocity space.
•

We are planning to propose a general framework on constructing stable and pressure-robust ${\bm{P}_{1}^{c}}-P0$ discretizations for incompressible fluid flows, including the Brinkman equations, the Oseen equations, the Navier-Stokes equations etc. Some related numerical results have been shown in this paper (see section 4.4).

Funding

This work was supported by the National Natural Science Foundation of China grant 12131014.

References

[1] Naveed Ahmed, Alexander Linke, and Christian Merdon, On really locking-free mixed finite element methods for the transient incompressible Stokes equations, SIAM J. Numer. Anal. 56 (2018), no. 1, 185–209.
[2] D. N. Arnold and J. Qin, Quadratic velocity/linear pressure Stokes elements, Advances in Computer Methods for Partial Differential Equations-VII, R. Vichnevetsky, D. Knight & G. Richter, eds., IMACS, New Brunswick, NJ (1992), 28–34.
[3] Douglas N. Arnold, Franco Brezzi, Bernardo Cockburn, and L. Donatella Marini, Unified analysis of discontinuous Galerkin methods for elliptic problems, SIAM J. Numer. Anal. 39 (2002), no. 5, 1749–1779 (en).
[4] Christine Bernardi and Genevieve Raugel, Analysis of some finite elements for the Stokes problem, Math. Comp. 44 (1985), no. 169, 71–79.
[5] Daniele Boffi, Franco Brezzi, and Michel Fortin, Mixed finite element methods and applications, Springer Series in Computational Mathematics, vol. 44, Springer Berlin Heidelberg, 2013.
[6] F. Brezzi, On the existence, uniqueness and approximation of saddle-point problems arising from lagrangian multipliers, ESAIM: Math. Model. Numer. Anal. 8 (1974), no. R2, 129–151.
[7] Franco Brezzi, Jim Douglas, and L. D. Marini, Two families of mixed finite elements for second order elliptic problems, Numer. Math. 47 (1985), no. 2, 217–235 (en).
[8] Erik Burman, Snorre H. Christiansen, and Peter Hansbo, Application of a minimal compatible element to incompressible and nearly incompressible continuum mechanics, Comput. Methods Appl. Mech. Engrg. 369 (2020), 113224.
[9] Snorre H. Christiansen and Kaibo Hu, Generalized finite element systems for smooth differential forms and Stokes’ problem, Numer. Math. 140 (2018), no. 2, 327–371.
[10] Philippe G. Ciarlet, The finite element method for elliptic problems, North-Holland, New York, 1978.
[11] Bernardo Cockburn, Jayadeep Gopalakrishnan, and Raytcho Lazarov, Unified hybridization of discontinuous galerkin, mixed, and continuous galerkin methods for second order elliptic problems, SIAM J. Numer. Anal. 47 (2009), no. 2, 1319–1365.
[12] Bernardo Cockburn, Guido Kanschat, and Dominik Sch ${\rm\ddot{o}}$ tzau, A note on discontinuous Galerkin divergence-free solutions of the Navier-Stokes equations, J. Sci. Comput. 31 (2007), no. 1-2, 61–73.
[13] Crouzeix, M. and Raviart, P.-A., Conforming and nonconforming finite element methods for solving the stationary Stokes equations I, R.A.I.R.O. 7 (1973), 33–75.
[14] Monique Dauge, Stationary Stokes and Navier-Stokes systems on two- or three-dimensional domains with corners. Part I. linearized equations, SIAM J. Math. Anal. 20 (1989), no. 1, 74–97.
[15] John A. Evans and Thomas J. R. Hughes, Isogeometric divergence-conforming B-splines for the steady Navier-Stokes equations, Math. Models Methods Appl. Sci. 23 (2013), no. 08, 1421–1478.
[16] by same author, Isogeometric divergence-conforming B-splines for the unsteady Navier-Stokes equations, J. Comput. Phys. 241 (2013), 141–167.
[17] Keith J. Galvin, Alexander Linke, Leo G. Rebholz, and Nicholas E. Wilson, Stabilizing poor mass conservation in incompressible flow problems with large irrotational forcing and application to thermal convection, Comput. Methods Appl. Mech. Engrg. 237-240 (2012), 166–176.
[18] Nicolas R. Gauger, Alexander Linke, and Phillip W. Schroeder, On high-order pressure-robust space discretisations, their advantages for incompressible high Reynolds number generalised Beltrami flows and beyond, SMAI J. Comput. Math. 5 (2019), 89–129.
[19] Vivette Girault and Jacques-Louis Lions, Two-grid finite-element schemes for the transient Navier-Stokes problem, ESAIM Math. Model. Numer. Anal. 35 (2001), no. 5, 945–980.
[20] Vivette Girault and Pierre-Arnaud Raviart, Finite element methods for Navier-Stokes equations, Springer Series in Computational Mathematics, vol. 5, Springer Berlin Heidelberg, Berlin, Heidelberg, 1986.
[21] J. Guzmán and M. Neilan, Conforming and divergence-free Stokes elements in three dimensions, IMA J. Numer. Anal. 34 (2014), 1489–1508.
[22] Johnny Guzmán and Michael Neilan, A family of nonconforming elements for the Brinkman problem, IMA J. Numer. Anal. 32 (2012), no. 4, 1484–1508.
[23] by same author, Inf-sup stable finite elements on barycentric refinements producing divergence–free approximations in arbitrary dimensions, SIAM J. Numer. Anal. 56 (2018), no. 5, 2826–2844.
[24] Volker John, Petr Knobloch, and Julia Novo, Finite elements for scalar convection-dominated equations and incompressible flow problems: a never ending story?, Comput. Visual Sci. 19 (2018), no. 5-6, 47–63.
[25] Volker John, Alexander Linke, Christian Merdon, Michael Neilan, and Leo G. Rebholz, On the divergence constraint in mixed finite element methods for incompressible flows, SIAM Rev. 59 (2017), no. 3, 492–544.
[26] Guido Kanschat and Natasha Sharma, Divergence-conforming discontinuous Galerkin methods and $c^{0}$ interior penalty methods, SIAM J. Numer. Anal. 52 (2014), no. 4, 1822–1842.
[27] Juho Könnö and Rolf Stenberg, $H(\operatorname{div})$ -conforming finite elements for the Binkman problem, Math. Models Methods Appl. Sci. 21 (2011), no. 11, 2227–2248.
[28] Philip L. Lederer, Christoph Lehrenfeld, and Joachim Schöberl, Hybrid discontinuous Galerkin methods with relaxed H(div)-conformity for incompressible flows. Part I, SIAM J. Numer. Anal. 56 (2018), no. 4, 2070–2094.
[29] by same author, Hybrid discontinuous Galerkin methods with relaxed H(div)-conformity for incompressible flows. Part II, ESAIM: Math. Model. Numer. Anal. 53 (2019), no. 2, 503–522.
[30] Philip L. Lederer, Alexander Linke, Christian Merdon, and Joachim Schöberl, Divergence-free reconstruction operators for pressure-robust Stokes discretizations with continuous pressure finite elements, SIAM J. Numer. Anal. 55 (2017), no. 3, 1291–1314.
[31] Philip L. Lederer and Sander Rhebergen, A pressure-robust embedded discontinuous Galerkin method for the Stokes problem by reconstruction operators, SIAM J. Numer. Anal. 58 (2020), no. 5, 2915–2933.
[32] Christoph Lehrenfeld, Hybrid discontinuous galerkin methods for solving incompressible flow problems, Master’s thesis, RWTH Aachen University, 2010.
[33] Christoph Lehrenfeld and Joachim Schöberl, High order exactly divergence-free hybrid discontinuous Galerkin methods for unsteady incompressible flows, Comput. Methods Appl. Mech. Engrg. 307 (2016), 339–361.
[34] A. Linke and C. Merdon, On velocity errors due to irrotational forces in the Navier-Stokes momentum balance, J. Comput. Phys. 313 (2016), 654–661.
[35] by same author, Pressure-robustness and discrete Helmholtz projectors in mixed finite element methods for the incompressible Navier-Stokes equations, Comput. Methods Appl. Mech. Engrg. 311 (2016), 304–326.
[36] Alexander Linke, Collision in a cross-shaped domain-A steady 2d Navier-Stokes example demonstrating the importance of mass conservation in CFD, Comput. Methods Appl. Mech. Engrg. 198 (2009), no. 41, 3278–3286.
[37] Alexander Linke, On the role of the Helmholtz decomposition in mixed methods for incompressible flows and a new variational crime, Comput. Methods Appl. Mech. Engrg. 268 (2014), 782–800.
[38] Alexander Linke, Gunar Matthies, and Lutz Tobiska, Robust arbitrary order mixed finite element methods for the incompressible Stokes equations with pressure independent velocity errors, ESAIM: M2AN 50 (2016), no. 1, 289–309.
[39] Kent Andre Mardal, Xue-Cheng Tai, and Ragnar Winther, A robust finite element method for Darcy–Stokes flow, SIAM J. Numer. Anal. 40 (2002), no. 5, 1605–1631.
[40] J. C. Nédélec, Mixed finite elements in R3, Numer. Math. 35 (1980), 315–341.
[41] Phillip Joseph Phillips and Mary F. Wheeler, A coupling of mixed and continuous Galerkin finite element methods for poroelasticity I: the continuous in time case, Comput. Geosci. 11 (2007), no. 2, 131–144.
[42] P. A. Raviart and J. M. Thomas, A mixed finite element method for 2nd order elliptic problems, in Mathematical Aspects of Finite Element Methods, I. Galligani, E. Magenes, eds., Lectures Notes in Math. 606 (1977), 292–315.
[43] Sander Rhebergen and Garth N. Wells, Analysis of a hybridized/interface stabilized finite element method for the Stokes equations, SIAM J. Numer. Anal. 55 (2017), no. 4, 1982–2003.
[44] Sander Rhebergen and Garth N. Wells, An embedded-hybridized discontinuous Galerkin finite element method for the Stokes equations, Comput. Methods Appl. Mech. Engrg. 358 (2020), 112619.
[45] C. Rodrigo, X. Hu, P. Ohm, J.H. Adler, F.J. Gaspar, and L.T. Zikatanov, New stabilized discretizations for poroelasticity and the Stokes’ equations, Comput. Methods Appl. Mech. Engrg. 341 (2018), 467–484 (en).
[46] L. R. Scott and M. Vogelius, Norm estimates for a maximal right inverse of the divergence operator in spaces of piecewise polynomials, RAIRO Modél. Math. Anal. Numér. 19 (1985), 111–143.
[47] X.-C. Tai and R. Winther, A discrete de Rham complex with enhanced smoothness, Calcolo 43 (2006), no. 4, 287–306.
[48] Junping Wang and Xiu Ye, New finite element methods in computational fluid dynamics by ${H}(\operatorname{div})$ elements, SIAM J. Numer. Anal. 45 (2007), no. 3, 1269–1286 (en).
[49] Xiaoping Xie, Jinchao Xu, and Guangri Xue, Uniformly-stable finite element methods for Darcy-Stokes-Brinkman models, J. Comput. Math. 26 (2008), no. 003, 437–455.
[50] Shangyou Zhang, A new family of stable mixed finite elements for the 3D Stokes equations, Math. Comp. 74 (2005), no. 250, 543–554.


(3.27a)	$\displaystyle\left\|(-\Delta\bm{u},\bm{v})-a_{h}\left(\bm{u},\bm{v}\right)\right\|$	$\displaystyle=\left\|(-\Delta\bm{u},\bm{v}^{R})\right\|\leq\sum_{T\in\mathcal{T}_{h}}h_{T}\\|\Delta\bm{u}\\|_{T}h_{T}^{-1}\left\\|\bm{v}^{R}\right\\|_{T}$
(3.27b)		$\displaystyle\leq(\sum_{T\in\mathcal{T}_{h}}h_{T}^{2}\\|\Delta\bm{u}\\|_{T}^{2})^{1/2}(\sum_{T\in\mathcal{T}_{h}}h_{T}^{-2}\left\\|\bm{v}^{R}\right\\|_{T}^{2})^{1/2}.$

	$\displaystyle\left\\|p_{h}-P_{h}p\right\\|$	$\displaystyle\leq\frac{1}{\beta_{1}}\sup_{\bm{v}_{h}\in V_{h}\setminus\{\bm{0}\}}\frac{b\left(\bm{v}_{h},p_{h}-P_{h}p\right)}{\|\|\|\bm{v}_{h}\|\|\|}=\frac{1}{\beta_{1}}\sup_{\bm{v}_{h}\in V_{h}\setminus\{\bm{0}\}}\frac{b\left(\bm{v}_{h},p_{h}-p\right)}{\|\|\|\bm{v}_{h}\|\|\|}$
		$\displaystyle=\frac{1}{\beta_{1}}\sup_{\bm{v}_{h}\in V_{h}\setminus\{\bm{0}\}}\frac{-\nu a_{h}\left(\bm{u}-\bm{u}_{h},\bm{v}_{h}\right)-\nu\delta_{h}(\bm{u},\bm{v}_{h})}{\|\|\|\bm{v}_{h}\|\|\|}$
		$\displaystyle\leq\frac{\nu}{\beta_{1}}\sup_{\bm{v}_{h}\in V_{h}\setminus\{\bm{0}\}}\frac{1}{\|\|\|\bm{v}_{h}\|\|\|}\|\|\|\bm{v}_{h}\|\|\|\left(\|\|\|\bm{u}-\bm{u}_{h}\|\|\|+C_{\alpha}h\|\bm{u}\|_{2}\right)\leq\frac{\nu}{\beta_{1}}\left(\|\|\|\bm{u}-\bm{u}_{h}\|\|\|+C_{\alpha}h\|\bm{u}\|_{2}\right).$

A low order divergence-free H(div)-conforming finite element method for Stokes flows

Abstract.

Key words and phrases:

1. Introduction

2. Model, notation and method

2.1. Discretizations

Lemma 1.

Proof.

Remark 1 (The motivation of aRa^{R}).

Remark 2.

Theorem 2.1 (Mass Conservation).

Proof.

Theorem 2.2.

Proof.

2.2. A pressure-robust and stabilized P1c−P​0{P_{1}^{c}}-P0 scheme

Remark 3.

Remark 4.

3. Error estimates

3.1. Preliminaries

Lemma 2.

Proof.

Lemma 3.

Proof.

3.2. Approximation

Lemma 4.

Proof.

3.3. Stability and boundedness

Lemma 5 (Coercivity).

Lemma 6 (Boundedness).

Lemma 7 (LBB Stability).

Proof.

3.4. Consistency Error

Lemma 8 (Consistency).

Proof.

3.5. Error estimates

Theorem 3.1.

Proof.

Remark 5.

Theorem 3.2.

Proof.

Remark 6.

Corollary 1.

Remark 7.

Remark 8.

4. Numerical experiments

4.1. Convergence test

4.2. Effects of different values of the parameters αT\alpha_{T}

4.3. Stokes flow with Coriolis force

4.4. The robustness of the pressure-robust P1c−P​0{P_{1}^{c}-P0} scheme

5. Conclusions

Funding

References

Remark 1 (The motivation of $a^{R}$ ).

2.2. A pressure-robust and stabilized ${P_{1}^{c}}-P0$ scheme

4.2. Effects of different values of the parameters $\alpha_{T}$

4.4. The robustness of the pressure-robust ${P_{1}^{c}-P0}$ scheme