Pressure-robustness in quasi-optimal a priori estimates for the Stokes problem

TL;DR

This paper establishes pressure-robust a priori error estimates for the Stokes problem that are valid under minimal regularity assumptions, improving the understanding of discretization accuracy in incompressible flow simulations.

Contribution

It provides the first general a priori estimate for pressure-robust discretizations of the Stokes problem applicable to data with minimal regularity.

Findings

Pressure-robust estimates depend only on the Helmholtz--Hodge projector of data.

The new estimates are valid without assuming high regularity of the velocity.

Numerical examples confirm the theoretical results.

Abstract

Recent analysis of the divergence constraint in the incompressible Stokes/Navier--Stokes problem has stressed the importance of equivalence classes of forces and how it plays a fundamental role for an accurate space discretization. Two forces in the momentum balance are velocity--equivalent if they lead to the same velocity solution, i.e., if and only if the forces differ by only a gradient field. Pressure-robust space discretizations are designed to respect these equivalence classes. One way to achieve pressure-robust schemes is to introduce a non-standard discretization of the right-side forcing term for any inf-sup stable mixed finite element method. This modification leads to pressure-robust and optimal-order discretizations, but a proof was only available for smooth situations and remained open in the case of minimal regularity, where it cannot be assumed that the vector Laplacian of the velocity is at least square-integrable. This contribution closes this gap by delivering a general estimate for the consistency error that depends only on the regularity of the data term. Pressure-robustness of the estimate is achieved by the fact that the new estimate only depends on the $L^2$ norm of the Helmholtz--Hodge projector of the data term and not on the $L^2$ norm of the entire data term. Numerical examples illustrate the theory.

Tables

Table 1. Table 1: Errors for the classical Bernardi–Raugel finite element method for 𝒇 = ∇ ( sin ⁡ ( x ​ y ​ π ) ) 𝒇 ∇ 𝑥 𝑦 𝜋 \bm{f}=\nabla(\sin(xy\pi)) and ν = 1 𝜈 1 \nu=1 .

ndof	${\Vert \nabla \left(𝒗-{𝒗}_h\right)\Vert }_{L^2\left(\Omega \right)}$	order	${\Vert \nabla \left({𝒗}_h-{𝑺}_h𝒗\right)\Vert }_{L^2\left(\Omega \right)}$	order
379	1.4151e+00	-	5.0351e-02	-
1414	9.7300e-01	0.542	3.0576e-02	0.722
5458	6.7235e-01	0.535	1.6366e-02	0.905
21442	4.6297e-01	0.540	8.4114e-03	0.964
84994	3.1819e-01	0.543	4.2567e-03	0.986
338434	2.1844e-01	0.545	2.1402e-03	0.996
1350658	1.4988e-01	0.546	1.0729e-03	1.000

Table 2. Table 2: Errors for the modified Bernardi–Raugel finite element method for 𝒇 = ∇ ( sin ⁡ ( x ​ y ​ π ) ) 𝒇 ∇ 𝑥 𝑦 𝜋 \bm{f}=\nabla(\sin(xy\pi)) and ν = 1 𝜈 1 \nu=1 .

ndof	${\Vert \nabla \left(𝒗-{𝒗}_h\right)\Vert }_{L^2\left(\Omega \right)}$	order	${\Vert \nabla \left({𝒗}_h-{𝑺}_h𝒗\right)\Vert }_{L^2\left(\Omega \right)}$
379	1.4142e+00	-	8.5800e-11
1414	9.7261e-01	0.542	1.2467e-13
5458	6.7218e-01	0.535	1.9887e-14
21442	4.6290e-01	0.540	4.3878e-14
84994	3.1816e-01	0.543	9.8787e-14
338434	2.1844e-01	0.545	2.2136e-13
1350658	1.4988e-01	0.546	4.4909e-13

Table 3. Table 3: Errors for the classical Bernardi–Raugel finite element method for 𝒇 = ∇ ( sin ⁡ ( x ​ y ​ π ) ) 𝒇 ∇ 𝑥 𝑦 𝜋 \bm{f}=\nabla(\sin(xy\pi)) and ν = 10 − 2 𝜈 superscript 10 2 \nu=10^{-2} .

ndof	${\Vert \nabla \left(𝒗-{𝒗}_h\right)\Vert }_{L^2\left(\Omega \right)}$	order	${\Vert \nabla \left({𝒗}_h-{𝑺}_h𝒗\right)\Vert }_{L^2\left(\Omega \right)}$	order
379	4.5794e+00	-	5.0351e+00	-
1414	2.8168e+00	0.704	3.0576e+00	0.722
5458	1.5663e+00	0.850	1.6366e+00	0.905
21442	8.6350e-01	0.862	8.4114e-01	0.964
84994	4.8756e-01	0.828	4.2567e-01	0.986
338434	2.8682e-01	0.768	2.1402e-01	0.996
1350658	1.7650e-01	0.703	1.0729e-01	1.000

Table 4. Table 4: Errors for the modified Bernardi–Raugel finite element method for 𝒇 = ∇ ( sin ⁡ ( x ​ y ​ π ) ) 𝒇 ∇ 𝑥 𝑦 𝜋 \bm{f}=\nabla(\sin(xy\pi)) and ν = 10 − 2 𝜈 superscript 10 2 \nu=10^{-2} .

ndof	${\Vert \nabla \left(𝒗-{𝒗}_h\right)\Vert }_{L^2\left(\Omega \right)}$	order	${\Vert \nabla \left({𝒗}_h-{𝑺}_h𝒗\right)\Vert }_{L^2\left(\Omega \right)}$
379	1.4142e+00	-	8.5800e-09
1414	9.7261e-01	0.542	1.2516e-11
5458	6.7218e-01	0.535	6.5365e-13
21442	4.6290e-01	0.540	1.3425e-12
84994	3.1816e-01	0.543	2.7291e-12
338434	2.1844e-01	0.545	5.5018e-12
1350658	1.4988e-01	0.546	1.1034e-11

Table 5. Table 5: Errors for the classical Bernardi–Raugel finite element method for 𝒇 = ∇ ( sin ⁡ ( x ​ y ​ π ) ) 𝒇 ∇ 𝑥 𝑦 𝜋 \bm{f}=\nabla(\sin(xy\pi)) and ν = 10 − 4 𝜈 superscript 10 4 \nu=10^{-4} .

ndof	${\Vert \nabla \left(𝒗-{𝒗}_h\right)\Vert }_{L^2\left(\Omega \right)}$	order	${\Vert \nabla \left({𝒗}_h-{𝑺}_h𝒗\right)\Vert }_{L^2\left(\Omega \right)}$	order
379	4.3469e+02	-	5.0351e+02	-
1414	2.6424e+02	0.721	3.0576e+02	0.722
5458	1.4134e+02	0.906	1.6366e+02	0.905
21442	7.2822e+01	0.960	8.4114e+01	0.964
84994	3.6908e+01	0.984	4.2567e+01	0.986
338434	1.8571e+01	0.995	2.1402e+01	0.996
1350658	9.3137e+00	0.999	1.0729e+01	1.000

Table 6. Table 6: Errors for the modified Bernardi–Raugel finite element method for 𝒇 = ∇ ( sin ⁡ ( x ​ y ​ π ) ) 𝒇 ∇ 𝑥 𝑦 𝜋 \bm{f}=\nabla(\sin(xy\pi)) and ν = 10 − 4 𝜈 superscript 10 4 \nu=10^{-4} .

ndof	${\Vert \nabla \left(𝒗-{𝒗}_h\right)\Vert }_{L^2\left(\Omega \right)}$	order	${\Vert \nabla \left({𝒗}_h-{𝑺}_h𝒗\right)\Vert }_{L^2\left(\Omega \right)}$
379	1.4142e+00	-	7.9830e-07
1414	9.7261e-01	0.542	1.1388e-09
5458	6.7218e-01	0.535	6.5440e-11
21442	4.6290e-01	0.540	1.3407e-10
84994	3.1816e-01	0.543	2.7257e-10
338434	2.1844e-01	0.545	5.4872e-10
1350658	1.4988e-01	0.546	1.0965e-09

Equations

$$\begin{array}{rl}-\nu \Delta \mathbf{v}+\nabla p& =\mathbf{f},\\ -\mathrm{div}\mathbf{v}& =g\end{array}$$

$${\mathbf{f}}_1\simeq {\mathbf{f}}_2,$$

$${\mathbf{f}}_1\simeq {\mathbf{f}}_2\iff \exists \varphi \in H^1(\Omega )/\mathbb{R}:{\mathbf{f}}_2={\mathbf{f}}_1+\nabla \varphi .$$

$$\mathbb{P}(\mathbf{f}):=\underset{\varphi \in H^1(\Omega )}{\mathrm{argmin}}\Vert \mathbf{f}-\nabla \varphi {\Vert }_{{\mathbf{L}}^2(\Omega )},$$

$$\Vert \nabla ({\mathbf{v}}_h-{\mathbf{S}}_h(\mathbf{v})){\Vert }_{L^2(\Omega )}\le C_A{h}^s\Vert \mathbb{P}(-\Delta \mathbf{v}){\Vert }_{L^2(\Omega )}+\frac{C_{B}h}{\nu }\Vert \mathbf{f}-\mathbb{P}(\mathbf{f}){\Vert }_{L^2(\Omega )},$$

$$\Vert \nabla ({\mathbf{v}}_h-{\mathbf{S}}_h(\mathbf{v})){\Vert }_{L^2(\Omega )}\le (C_{1}h+C_A{h}^s)\Vert \mathbb{P}(-\Delta \mathbf{v}){\Vert }_{L^2(\Omega )},$$

$$\mathbf{H}(\mathrm{div};\Omega )$$

$${\mathbf{H}}_0(\mathrm{div};\Omega )$$

$$-\nu \Delta \mathbf{v}+\nabla p=\mathbf{f},\text{and}\mathrm{div}\mathbf{v}=0\text{in }\Omega .$$

$$\begin{array}{rl}\nu (\nabla \mathbf{v},\nabla \mathbf{w})-(p,\mathrm{div}\mathbf{w})& =(\mathbf{f},\mathbf{w}),\\ (\mathrm{div}\mathbf{v},q)& =0\end{array}$$

$${\mathbf{V}}^0:=\{\mathbf{w}\in {\mathbf{H}}_0^1(\Omega ):\nabla \cdot \mathbf{w}=0\}.$$

$$\mathbf{f}=\nabla \alpha +\mathbb{P}(\mathbf{f}),$$

$$\mathbb{P}(\mathbf{f})\in {\mathbf{L}}_{\sigma }^2(\Omega ):=\{\mathbf{w}\in {\mathbf{L}}^2(\Omega ):(\nabla q,\mathbf{w})=0\text{ for all }q\in H^1(\Omega )\}$$

$$(\mathbb{P}(\mathbf{f}),\mathbf{w})=(\mathbf{f},\mathbf{w})\text{for all }\mathbf{w}\in {\mathbf{L}}_{\sigma }^2(\Omega ).$$

$$\nu (\nabla \mathbf{v},\nabla \mathbf{w})=(\mathbf{f},\mathbf{w})=(\mathbb{P}(\mathbf{f}),\mathbf{w})\text{for all }\mathbf{w}\in {\mathbf{V}}^0.$$

$$<\mathbb{P}(\mathbf{f}),\mathbf{w}>=<\mathbf{f},\mathbf{w}>\text{for all }\mathbf{w}\in {\mathbf{V}}^0.$$

$$<\mathbb{P}(\hat{\mathbf{f}}),\mathbf{w}>=<\hat{\mathbf{f}},\mathbf{w}>=(\mathbf{f},\mathbf{w})=(\mathbb{P}(\mathbf{f}),\mathbf{w}).$$

$$<-\Delta \mathbf{w},\mathbf{\psi }>:=(\nabla \mathbf{w},\nabla \mathbf{\psi })\text{for all }\mathbf{\psi }\in {\mathbf{H}}_0^1(\Omega ).$$

$$\mathbb{P}(-\Delta \mathbf{v})=\frac{1}{\nu }\mathbb{P}(\mathbf{f}).$$

$$\Vert \nabla \mathbf{v}{\Vert }_{L^2(\Omega )}\le \frac{C_{PF}}{\nu }\Vert \mathbb{P}(\mathbf{f}){\Vert }_{L^2(\Omega )}=C_{PF}\Vert \mathbb{P}(\Delta \mathbf{v}){\Vert }_{L^2(\Omega )}$$

$$({\pi }_{h}r,q_h)=(r,q_h)\text{for all qh∈Qh.}$$

$$\Vert r-{\pi }_{h}r{\Vert }_{L^2(\Omega )}=\underset{q_h\in Q_h}{\inf }\Vert r-q_h{\Vert }_{L^2(\Omega )}\le C_{{\pi }_h,s}{h}^s\Vert r{\Vert }_{H^s(\Omega )}$$

$${\mathbf{V}}_h^0:=\{{\mathbf{v}}_h\in {\mathbf{X}}_h:{\mathrm{div}}_h{\mathbf{v}}_h=0\}.$$

$$(i)$$

$$(ii)$$

$$\nu (\nabla {\mathbf{v}}_h,\nabla {\mathbf{w}}_h)-({\mathrm{div}}_h{\mathbf{w}}_h,p_h)$$

$$({\mathrm{div}}_h{\mathbf{v}}_h,q_h)$$

$$\nu (\nabla {\mathbf{v}}_h,\nabla {\mathbf{w}}_h)=(\mathbf{f},{\mathbf{I}}_h\mathbf{w})=(\mathbb{P}(\mathbf{f}),{\mathbf{I}}_h{\mathbf{w}}_h)\text{for all }{\mathbf{w}}_h\in {\mathbf{V}}_h^0,$$

$$\Vert \nabla {\mathbf{v}}_h{\Vert }_{L^2(\Omega )}\le (C_{PF}+C_{1}h)\Vert {\nu }^{-1}\mathbb{P}(\mathbf{f}){\Vert }_{L^2(\Omega )}=(C_{PF}+C_{1}h)\Vert \mathbb{P}(\Delta \mathbf{v}){\Vert }_{L^2(\Omega )}.$$

$$\Vert \nabla {\mathbf{v}}_h{\Vert }_{L^2(\Omega )}$$

$$\le C_{PF}\Vert \mathbb{P}(\Delta \mathbf{v}){\Vert }_{L^2(\Omega )}+\frac{Ch}{\nu }\Vert \mathbf{f}-\mathbb{P}(\mathbf{f}){\Vert }_{L^2(\Omega )}.$$

$$\nu \Vert \nabla {\mathbf{v}}_h{\Vert }_{L^2(\Omega )}^2=(\mathbb{P}(\mathbf{f}),{\mathbf{I}}_h{\mathbf{v}}_h)$$

Full text

1 Introduction

Classical mixed finite element theory for the steady Stokes problem

[TABLE]

with inhomogeneous Dirichlet boundary data, $\bm{f}\in\bm{L}^{2}(\Omega)$ and $g\in L^{2}_{0}(\Omega)$ emphasizes that the divergence constraint $-{\rm div}\,{\bm{v}}=g$ requires an appropriate discrete mimicking of the surjectivity of the divergence operator $\mathrm{div}\!\!:{\bm{H}}^{1}_{0}(\Omega)\to L^{2}_{0}(\Omega)$ in order to guarantee optimal convergence properties, see e.g. [3, 10]. Recently it has been stressed that the divergence constraint in the Stokes problem naturally induces a semi-norm and corresponding equivalence classes of forces, which require a second challenge for an accurate space discretization: two forces $\bm{f}_{1}\in\bm{L}^{2}(\Omega)$ and $\bm{f}_{2}\in\bm{L}^{2}(\Omega)$ are velocity-equivalent [7]

[TABLE]

if they lead to the same velocity solution ${\bm{v}}$ in the Stokes problem (1) — and this happens if and only if both forces differ by a gradient field [10, 1], i.e.,

[TABLE]

The argument is straightforward: denote by $({\bm{v}}_{1},p_{1})$ and $({\bm{v}}_{2},p_{2})$ the pairs of velocity and pressure solutions corresponding to the forces $\bm{f}_{1}$ and $\bm{f}_{2}=\bm{f}_{1}+\nabla\phi$. Then, the difference of the solutions $(\delta{\bm{v}},\delta p):=({\bm{v}}_{2}-{\bm{v}}_{1},p_{2}-p_{1})\in\bm{H}^{1}_{0}(\Omega)\times L^{2}_{0}(\Omega)$ fulfills the incompressible Stokes equations $-\nu\Delta(\delta{\bm{v}})+\nabla(\delta p)=\nabla\phi$, $\mathrm{div}(\delta{\bm{v}})=0$ with homogeneous Dirichlet boundary data. This problem has the unique solution $(\delta{\bm{v}},\delta p)=(\bm{0},\phi)$, and thus $\bm{f}_{1}$ and $\bm{f}_{2}=\bm{f}_{1}+\nabla\phi$ are velocity-equivalent due to $\delta{\bm{v}}=\bm{0}$.

In conclusion one observes that the velocity solution ${\bm{v}}$ of (1) is determined by the following data:

Dirichlet boundary data, 2. 2.

the data $g$, 3. 3.

and the Helmholtz–Hodge projector of the data $\bm{f}$, which is defined by

[TABLE]

while the data term $\bm{f}-\mathbb{P}(\bm{f})$ only influences the pressure.

The recently introduced notion pressure-robustness [14] allows to discriminate between space discretizations for (1), whose discrete velocity solutions ${\bm{v}}_{h}$ depend on $\mathbb{P}(\bm{f})$ and not on the entire data $\bm{f}$. Such schemes lead to a priori error estimates for the discrete velocity that depend only on ${\bm{v}}$ and not on $({\bm{v}},\frac{1}{\nu}p)$ — as in nearly all classical mixed finite element methods [10].

This contribution focuses now on applying the improved understanding of relevant data in the Stokes problem, in order to derive a priori error estimates for various discretely inf–sup stable mixed methods in cases of minimal regularity. A special focus is set on a recent modified pressure-robust mixed method [13, 11], where the modification introduces a consistency error that can be optimally estimated in a straightforward manner by $Ch^{k}|\Delta{\bm{v}}|_{H^{k-1}(\Omega)}$ provided that ${\bm{v}}\in{\bm{H}}^{k-1}(\Omega)$. For the lowest–order methods $(k=1)$ this requires $\Delta{\bm{v}}\in{\bm{L}}^{2}(\Omega)$. In situations of minimal regularity, i.e., ${\bm{v}}\in{\bm{H}}^{1+s}(\Omega)$ with $0<s<1$, we provide an estimation of the consistency error by a more sophisticated argument involving the Helmholtz–Hodge projector of the data $\nu^{-1}\mathbb{P}(\bm{f})$. This term is obviously in ${\bm{L}}^{2}(\Omega)$, whenever it holds $\bm{f}\in{\bm{L}}^{2}(\Omega)$ and it is shown to be equal to $\mathbb{P}(-\Delta{\bm{v}})$. Thus, although it holds in general that $\Delta{\bm{v}}\not\in{\bm{L}}^{2}(\Omega)$ one can exploit in the numerical analysis that at least the divergence–free part of $\Delta{\bm{v}}$ is in ${\bm{L}}^{2}(\Omega)$.

This observation also leads to a seemingly new estimate for classical mixed methods, which can be sharper than classical a priori estimates, see Theorem 6.1. Eventually, all classical conforming finite element methods yield an estimate of the form

[TABLE]

while their pressure-robust siblings allow for estimates of the form

[TABLE]

with $C_{1}>0$, $C_{A}>0$ and $C_{B}>0$ are constants that do not depend on $h$. Note that for divergence-free conforming methods, see e.g. [16, 9], it holds $C_{1}=C_{A}=C_{B}=0$, but for them the only nontrivial part of the numerical analysis is the proof of the discrete inf-sup stability. Further, structurally identical results are obtained for the classical and a modified pressure-robust nonconforming Crouzeix–Raviart finite element method.

The rest of this paper is structured as follows. Section 2 introduces the Stokes problem as well as the framework for the modified finite element method and the assumptions that are crucial for the theoretical results. Section 3 focusses on the Helmholtz–Hodge projector and its application in stability estimates. Section 5 introduces the continuous and discrete Stokes projectors and their properties. Section 6 applies the tools of the previous sections to obtain quasi-optimal estimates for classical finite element methods that only depend on the data. Section 7 does the same for the modified pressure-robust finite element methods where now the error is additionally independent of the pressure and the inverse of the viscosity $\nu$. Section 8 revisits quasi-optimal and pressure-robust error estimates for the nonconforming Crouzeix–Raviart finite element method. Finally we perform some numerical experiments in Section 9 and compare these empirical results with the theory.

2 Preliminary results

This section introduces some notation, recalls some preliminaries and formulates an assumption that is fundamental for the presented theory. We adopt standard space notation and denote vector–valued functions and vector–valued function spaces in boldface. We use $(\cdot,\cdot)$ to denote the $L^{2}$-inner product over $\Omega\subset\mathbb{R}^{n}$, and by $\langle\cdot,\cdot\rangle$ the duality pairing between some Hilbert space and its dual. We denote by $L^{2}_{0}(\Omega)$ the Hilbert-space of square-integrable scalar functions with zero average, and

[TABLE]

where $\bm{n}$ denotes the outward unit normal of ${\partial}\Omega$.

2.1 Stokes problem and weak elliptic regularity assumption

In the following, we study finite element methods for the model problem: for $\bm{f}\in\bm{L}^{2}(\Omega)$ seek $({\bm{v}},p)\in{\bm{H}}^{1}_{0}(\Omega)\times L^{2}_{0}(\Omega)$ such that it holds

[TABLE]

The extension to the more general divergence constraint ${\rm div}\,{\bm{v}}=g$ with $g\in L^{2}_{0}(\Omega)$ is straightforward, and we refer to [10] for details.

A weak formulation of the problem is given by: search for $({\bm{v}},p)\in{\bm{H}}^{1}_{0}(\Omega)\times L^{2}_{0}(\Omega)$ such that it holds

[TABLE]

for all $(\bm{w},q)\in{\bm{H}}^{1}_{0}(\Omega)\times L^{2}_{0}(\Omega)$.

The space of divergence-free ${\bm{H}}^{1}_{0}(\Omega)$ vector fields is denoted as

[TABLE]

Assumption 2.1

Throughout the paper, we assume that the Stokes problem inherits ${\bm{H}}^{1+s}(\Omega)\times H^{s}(\Omega)$ elliptic regularity for some $s\in(0,1]$ and that $\nu\|{\bm{v}}\|_{H^{1+s}(\Omega)}+\|p\|_{H^{s}(\Omega)}\leq C_{\mathrm{ell},s}\|\bm{f}\|_{L^{2}(\Omega)}$.

3 Helmholtz–Hodge projector

According to the $\bm{L}^{2}$-orthogonal Helmholtz–Hodge decomposition (see e.g. [8]) any vector field $\bm{f}\in\bm{L}^{2}(\Omega)$ can be uniquely decomposed into

[TABLE]

where $\alpha\in H^{1}(\Omega)/\mathbb{R}$, and

[TABLE]

is the Helmholtz–Hodge projector of $\bm{f}$. Note, that the Helmholtz–Hodge projector of ${\bm{f}}$ is divergence–free and is the orthogonal ${\bm{L}}^{2}$ projection of $\bm{f}$ onto $\bm{L}^{2}_{\sigma}(\Omega)$, i.e.,

[TABLE]

Moreover, for the Stokes velocity solution ${\bm{v}}$ it holds

[TABLE]

The domain of the Helmholtz–Hodge projector can be extended to ${\bm{H}}^{-1}(\Omega)$ with range in $({\bm{V}}^{0})^{*}$, the space of bounded linear functionals on ${\bm{V}}^{0}$. Indeed, for every functional $\bm{f}\in{\bm{H}}^{-1}(\Omega)$ the Helmholtz–Hodge projector can be defined as the restriction to ${\bm{V}}^{0}$, i.e., it holds

[TABLE]

Condition (9) defines an extension of the Helmholtz–Hodge projector from $\bm{L}^{2}(\Omega)$ to ${\bm{H}}^{-1}(\Omega)$. Assume that the functional $\hat{\bm{f}}\in{\bm{H}}^{-1}(\Omega)$ has a representation $\bm{f}\in{\bm{L}}^{2}(\Omega)$ with $\bm{f}=\nabla\alpha+\mathbb{P}\bm{f}$. Then it holds for all $\bm{w}\in{\bm{V}}^{0}$

[TABLE]

Lemma 3.1**.**

Denote by $-\Delta:{\bm{H}}^{1}_{0}(\Omega)\to{\bm{H}}^{-1}(\Omega)$ via

[TABLE]

Then the weak velocity solution ${\bm{v}}$ of (5) satisfies

[TABLE]

Proof 3.2**.**

This follows directly from a combination of (8) and (9).

Thus, although the regularity of the functional $-\Delta{\bm{v}}$ is not better in general than $-\Delta{\bm{v}}\in{\bm{H}}^{-1}(\Omega)$, its divergence–free part $\mathbb{P}(-\Delta{\bm{v}})$ has the better regularity ${\bm{L}}^{2}(\Omega)$.

Remark 3.3**.**

We emphasize that Lemma 3.1 is of central importance for the derivation of pressure-robust a priori error estimates in case of minimal regularity. We also stress that the quantity $\nu^{-1}\mathbb{P}(\bm{f})$, which appears naturally in the analysis of pressure-robust methods, does in fact not scale with the inverse of $\nu$, since it only depends on ${\bm{v}}$.

An immediate consequence from Lemma 3.1 is the following result that bounds the norm of the velocity field by the norm of the Helmholtz–Hodge projector of the data $\bm{f}$.

Lemma 3.4** (Continuous stability estimate).**

The exact solution of problem  (4) satisfies

[TABLE]

where $C_{PF}$ is the constant from the Poincaré–Friedrichs inequality.

Proof 3.5**.**

The result follows directly from testing (8) with $\bm{w}={\bm{v}}$ and using the Poincaré–Friedrichs inequality.

Remark 3.6**.**

Here, we emphasize that the right hand side of the stability estimate is given by a semi-norm of the data $\bm{f}$. This is a crucial point, which arguably has not been fully exploited in classical mixed theory [3, 8].

4 Notation and setting for conforming finite element methods

In the following, we introduce some notation for the finite element methods used in this contribution. We denote by ${\bm{X}}_{h}\times Q_{h}$, a discretely inf-sup stable finite element pair [3] for the Stokes problem with homogeneous Dirichlet boundary conditions with respect to a conforming, shape–regular and simplicial triangulation $\mathcal{T}_{h}$ with $h=\max_{T\in\mathcal{T}_{h}}{\rm diam}(T)$. The $L^{2}$ best approximation onto the discrete pressure space $L^{2}_{0}(\Omega)$ is denoted by $\pi_{h}:L^{2}_{0}(\Omega)\to Q_{h}$, i.e., for all $r\in L^{2}_{0}(\Omega)$ it holds

[TABLE]

We assume that $Q_{h}$ has the approximation property

[TABLE]

for all $r\in H^{s}(\Omega)\cap L^{2}_{0}(\Omega)$ and $s\in(0,1]$.

Let $\mathrm{div}_{h}:{\bm{X}}_{h}\to Q_{h}$ with $\mathrm{div}_{h}=\pi_{h}{\rm div}\,$ denote the discrete divergence operator. Due to the assumed discrete inf–sup stability of the pair ${\bm{X}}_{h}\times Q_{h}$, $\mathrm{div}_{h}$ is surjective with bounded right–inverse [3]. We define the space of discretely divergence–free functions as

[TABLE]

4.1 Some modified finite element methods

As shown in [12, 13, 11], a certain modification of the discrete right–hand side of the incompressible Stokes problem renders inf-sup stable mixed methods pressure-robust. These pressure-robust finite element methods employ a reconstruction operator with the properties stated in the following assumption.

Assumption 4.1

We assume that there exists an auxiliary finite element space $\bm{Y}_{h}\subset{\bm{H}}_{0}({\rm div};\Omega)$ and a reconstruction operator ${\bm{I}}_{h}:{\bm{H}}^{1}_{0}(\Omega)\to\bm{Y}_{h}$ such that

[TABLE]

where $C_{1}$ depends only on the shape regularity of the mesh.

The modified finite element method for the Stokes problem applies the reconstruction operator in the right-hand side. The resulting scheme seeks $({\bm{v}}_{h},p_{h})\in{\bm{X}}_{h}\times Q_{h}$ such that

[TABLE]

Testing (16) with discretely divergence-free velocity test functions yields

[TABLE]

since for $\bm{w}_{h}\in{\bm{V}}^{0}_{h}$ it holds ${\bm{I}}_{h}\bm{w}_{h}\in\bm{L}^{2}_{\sigma}(\Omega)$. This last identity is characteristic for pressure-robustness and in general not true for non-divergence-free classical finite element methods. It tells us that the discrete velocity solution ${\bm{v}}_{h}$ of (16) depends on the appropriate continuous data $\nu^{-1}\mathbb{P}(\bm{f})$ of the problem.

In the case of discontinuous pressure spaces $Q_{h}$, the standard interpolation operators of the Raviart-Thomas or Brezzi–Douglas–Marini finite element spaces can be employed as a reconstruction operator ${\bm{I}}_{h}$, see [10, 14, 13] for details. For instance, in the case of the Bernardi–Raugel finite element method [2], the standard interpolator into the BDM space of order one can be used. For continuous pressure spaces, the design of the reconstruction operator is more involved; see [11] for details in case of the Taylor–Hood or MINI finite element family.

Remark 4.1**.**

Note, that for ${\bm{I}}_{h}=\bm{1}$ (the identity operator) in (16) the classical finite element method is obtained. However, only divergence-free $H^{1}$-conforming classical finite element methods, see e.g. [16, 9], satisfy Assumption 4.1 with $C_{1}=0$. In the results below it will be specified which results rely on this assumption.

Lemma 4.2** (Discrete stability estimates).**

Let $({\bm{v}}_{h},p_{h})\in{\bm{X}}_{h}\times Q_{h}$ satisfy (16) and write $\bm{f}=\nabla\alpha+\mathbb{P}\bm{f}$. Then if the discrete scheme satisfies Assumption 4.1, it holds the estimate

[TABLE]

If the discrete scheme with ${\bm{I}}_{h}=\bm{1}$ does not satisfy Assumption 4.1, it only holds

[TABLE]

Proof 4.3**.**

Testing (17) with $\bm{w}_{h}={\bm{v}}_{h}$, a discrete Poincaré–Friedrichs inequality and (15) yield

[TABLE]

If ${\bm{I}}_{h}=\bm{1}$ and Assumption 4.1 is not satisfied then inserting the Helmholtz–Hodge decomposition of $\bm{f}$ and an integration by parts give

[TABLE]

Property (13) shows $\|\alpha-\pi_{h}\alpha\|_{L^{2}(\Omega)}\leq Ch\|\nabla\alpha\|_{L^{2}(\Omega)}=Ch\|\bm{f}-\mathbb{P}(\bm{f})\|_{L^{2}(\Omega)}$. This concludes the proof.

5 Continuous and discrete Stokes projectors

In preparation for the a priori error estimates, this section studies the continuous and the discrete Stokes projectors. They are defined as the ${\bm{H}}^{1}$-seminorm best-approximations into the (discretely) divergence-free functions, i.e. ${\bm{S}}_{h}:{\bm{H}}^{1}_{0}(\Omega)\to{\bm{V}}^{0}_{h}$ and ${\bm{S}}:{\bm{X}}_{h}\to{\bm{V}}^{0}$ are defined by

[TABLE]

The rest of this section collects useful properties of these projectors.

Lemma 5.1** (Stokes projector identity).**

For any ${\bm{v}}\in{\bm{H}}^{1}_{0}(\Omega)$ and ${\bm{v}}_{h}\in{\bm{X}}_{h}$, it holds the identity

[TABLE]

Proof 5.2**.**

This follows directly from the combination of the definitions of ${\bm{S}}_{h}$ and ${\bm{S}}$.

Lemma 5.3**.**

Suppose that the Stokes problem satisfies Assumption 2.1. Then there holds

[TABLE]

Proof 5.4**.**

Let $(\bm{\psi},r)\in{\bm{H}}^{1}_{0}(\Omega)\times L^{2}_{0}(\Omega)$ solve the Stokes problem with source $\bm{w}_{h}-{\bm{S}}(\bm{w}_{h})$ and unit viscosity:

[TABLE]

Testing the first equation with $\bm{z}=\bm{w}_{h}-{\bm{S}}(\bm{w}_{h})$ and employing (19) leads to

[TABLE]

Recall that $\pi_{h}r$ is the $L^{2}$–projection of $r$ defined by (12), and note that it holds $(\nabla\cdot\bm{w}_{h},\pi_{h}{r}_{h})=0$ since $\bm{w}_{h}\in{\bm{V}}^{0}_{h}$. Consequently, by (13), we have

[TABLE]

Finally, the elliptic regularity Assumption 2.1 implies $\|r\|_{H^{s}(\Omega)}\leq C_{\mathrm{ell},s}\|\bm{w}_{h}-{\bm{S}}(\bm{w}_{h})\|_{L^{2}(\Omega)}$, and so

[TABLE]

Dividing the last inequality by $\|\bm{w}_{h}-{\bm{S}}(\bm{w}_{h})\|_{L^{2}(\Omega)}$ gets the desired result.

6 Quasi-optimal a priori error estimates for classical

finite element methods

This section derives a priori error estimates for classical finite element methods that are not pressure-robust, i.e. do not satisfy Assumption 4.1 with ${\bm{I}}_{h}=\bm{1}$ like the Bernardi–Raugel, MINI or Taylor–Hood finite element methods. The proof of the estimate bounds the error of the best-approximation by the right-hand side data.

Theorem 6.1**.**

Suppose that the Stokes problem satisfies Assumption 2.1, the reconstruction operator is taken to be the identity ${\bm{I}}_{h}=\bm{1}$, and that ${\bm{I}}_{h}$ does not satisfy Assumption 4.1. Then there holds

[TABLE]

with $C_{2}>0$ given by (20).

Proof 6.2**.**

Write ${\bm{e}}_{h}:={\bm{v}}_{h}-{\bm{S}}_{h}({\bm{v}})$ and note that ${\bm{e}}_{h}\in{\bm{V}}^{0}_{h}$. Hence, it follows from Lemmas 5.1 and 5.3 that

[TABLE]

where $\alpha$ stems from the Helmholtz–Hodge decomposition (7) of $\bm{f}$. The best approximation property of $\pi_{h}\alpha$ shows $\|\alpha-\pi_{h}\alpha\|_{L^{2}(\Omega)}\leq C_{\pi_{h},1}h\|\nabla\alpha\|_{L^{2}(\Omega)}=C_{\pi_{h},1}h\|\bm{f}-\mathbb{P}(\bm{f})\|_{L^{2}(\Omega)}$. This concludes the proof.

Remark 6.3**.**

Classical results for conforming mixed methods [8] show the a priori estimate

[TABLE]

*which scales like $\nu^{-1}h^{s}$ under the given regularity assumptions. Such an estimate is sometimes sharper than Theorem 6.1, but can also be less sharp.

i) If it holds, e.g., $p\in Q_{h}$, then the error on the right hand side of the classical estimate is zero. This is also preserved in the computations for the new estimate until (21), since in the special case $\bm{f}=-\nu\Delta{\bm{v}}$ (21) can be shown to vanish identically.

ii) If it holds $p\not\in Q_{h}$ and if the solution $({\bm{v}},p)\in{\bm{H}}^{1+s}(\Omega)\times H^{s}(\Omega)$ has a low regularity with $s<1$, then the new estimate can be sharper e.g. for $\nu\ll 1$, since it predicts an a priori error $\mathcal{O}(h^{s}+\nu^{-1}h)$, while the classical estimate predicts an error decay like $\mathcal{O}(\nu^{-1}h^{s})$. We remark that the pressure-dependent consistency error is influenced by two different contributions, one determined by $-\Delta{\bm{v}}$ and another one determined by $\frac{1}{\nu}({\bm{f}}-\mathbb{P}({\bm{f}}))$.*

Theorem 6.4** (A priori error estimate).**

Under the assumptions of Theorem 6.1, it holds

[TABLE]

Proof 6.5**.**

The proof starts with the Pythagoras theorem (using (18))

[TABLE]

The second term can be estimated by Theorem 6.1 and the first term can be bounded by the best-approximation error in ${\bm{X}}_{h}$ by the standard argument

[TABLE]

where $C_{F}\geq 1$ denotes the stability constant of the Fortin operator of the mixed method, see e.g. [10, 8].

7 Quasi-optimal pressure-robust a priori error estimates

This section concerns novel quasi-optimal a priori error estimates for conforming divergence-free and pressure-robustly modified finite element methods. Here, the distance between the discrete solution and the discrete Stokes projector can be bounded by $\|\mathbb{P}(-\Delta{\bm{v}})\|_{L^{2}(\Omega)}$ which is in general much smaller than the bound in Theorem 6.1, especially for small $\nu$.

Theorem 7.1**.**

Suppose that the Stokes problem satisfies Assumption 2.1 and that the reconstruction operator ${\bm{I}}_{h}$ satisfies Assumption 4.1. Then there holds

[TABLE]

with $C_{1}>0$ and $C_{2}>0$ given by (15) and (20), respectively. Note, that there is no dependency on $\nu^{-1}$.

Proof 7.2**.**

Write ${\bm{e}}_{h}:={\bm{v}}_{h}-{\bm{S}}_{h}({\bm{v}})$ and note that ${\bm{e}}_{h}\in{\bm{V}}^{0}_{h}$. Hence,

[TABLE]

The latter term is split up into (using also Lemma 5.1)

[TABLE]

It then follows from Lemma 5.3 and (15) that

[TABLE]

This concludes the proof.

Theorem 7.3** (A priori error estimate).**

Under the assumptions of Theorem 7.1, it holds

[TABLE]

Proof 7.4**.**

The proof starts with the Pythagoras theorem (using (18))

[TABLE]

The second term can be estimated by Theorem 7.1 and the first term can be bounded by the best-approximation error in ${\bm{X}}_{h}$ by the standard argument

[TABLE]

where $C_{F}\geq 1$ denotes the stability constant of the Fortin operator of the mixed method, see e.g. [10, 8].

8 Estimates for the nonconforming Crouzeix–Raviart finite element method

In this section we consider the space ${\bm{X}}_{h}\not\subset{\bm{H}}^{1}_{0}(\Omega)$ of nonconforming Crouzeix-Raviart functions, i.e., piecewise affine vector fields that are weakly continuous across edges (2D) or faces (3D) in the triangulation, see e.g. [6, 5]. To describe this space in detail we require some notation. Recall that $\mathcal{T}_{h}$ is a conforming, shape–regular, and simplicial triangulation of $\Omega$ parameterized by $h=\max_{T\in\mathcal{T}_{h}}{\rm diam}(T)$. We denote by $\mathcal{E}_{h}$ the set of $(n-1)$–dimensional simplices in $\mathcal{T}_{h}$, i.e., $\mathcal{E}_{h}$ is either the set of edges (2D) or faces (3D) in $\mathcal{T}_{h}$. Let $P_{m}(T)$ denote the space of polynomials of degree $\leq m$ on $T$, and let ${\bm{P}}_{m}(T)=(P_{m}(T))^{n}$. Then the Crouzeix-Raviart space ${\bm{X}}_{h}$ consists of all functions $\bm{w}_{h}\in{\bm{L}}^{2}(\Omega)$ with the properties $\bm{w}_{h}|_{T}\in{\bm{P}}_{1}(T)$, $\int_{E}\bm{w}_{h}$ is single–valued for all $E\in\mathcal{E}_{h}$, and $\int_{E}\bm{w}_{h}=0$ for all boundary $E\in\mathcal{E}_{h}$. The discrete pressure space $Q_{h}$ is the space of piecewise constants with vanishing mean. It is well–known that the pair ${\bm{X}}_{h}\times Q_{h}$ is inf–sup stable.

Note that Crouzeix-Raviart functions $\bm{w}_{h}\in{\bm{V}}^{0}_{h}$ are not divergence-free in a ${\bm{H}}(\mathrm{div})$-sense (as their normal traces are not continuous), but their piecewise divergence vanishes. Possible ${\bm{H}}(\mathrm{div})$-conforming reconstruction operators ${\bm{I}}_{h}$ for this method are the lowest-order Raviart–Thomas or BDM interpolation operators, see [4] for details.

In order to show the same quasi-optimal a priori error estimates for the Crouzeix–Raviart method some arguments have to be slightly modified. First, the Stokes projectors ${\bm{S}}_{h}:{\bm{H}}^{1}_{0}(\Omega)\to{\bm{V}}^{0}_{h}$ and ${\bm{S}}:{\bm{X}}_{h}\to{\bm{V}}^{0}$ are now defined by using the piecewise gradients $\nabla_{h}$, i.e.,

[TABLE]

Recall the Crouzeix–Raviart Fortin interpolation

[TABLE]

which satisfies the approximation property

[TABLE]

for all $s\in[0,1]$. This definition of the interpolant yields the well–known property [5]

[TABLE]

and in particular $\int_{T}\mathrm{div}({\bm{v}}-{\bm{I}}_{\text{CR}}{\bm{v}})=0$ for any $T\in\mathcal{T}_{h}$. Since $\nabla_{h}\bm{w}_{h}$ is piecewise constant this also reveals that we have ${\bm{S}}_{h}={\bm{I}}_{\text{CR}}$, i.e., the Crouzeix–Raviart interpolator is the discrete Stokes projector. Also note that the Stokes projector identity holds in the form

[TABLE]

However, in general ${\bm{I}}_{\text{CR}}{\bm{v}}\in{\bm{X}}_{h}$ does not imply ${\bm{I}}_{\text{CR}}{\bm{v}}\in{\bm{H}}(\mathrm{div},\Omega)$ and therefore Lemma 5.3 has to be modified as well.

The analysis also needs another mapping that projects a discretely divergence–free Crouzeix–Raviart function to some $H^{1}$-conforming divergence-free function. Such an operator was introduced in [15] and is based on rational bubble functions.

Lemma 8.1**.**

Suppose that the Stokes problem satisfies Assumption 2.1. Then there holds

[TABLE]

Proof 8.2**.**

Consider the $H^{1}_{0}$-conforming and $H^{1}$-stable operator $\bm{E}_{h}$ from [15] with the properties

[TABLE]

The second property follows from [15, $I_{2}=0$ in proof of Theorem 5.1].As in Lemma 5.3 we look at the solution $(\bm{\psi},r)\in{\bm{H}}^{1}_{0}(\Omega)\times L^{2}_{0}(\Omega)$ of the Stokes problem with modified source $\bm{E}_{h}\bm{w}_{h}-{\bm{S}}(\bm{w}_{h})$ and unit viscosity:

[TABLE]

Testing the first equation with $\bm{z}=\bm{E}_{h}\bm{w}_{h}-{\bm{S}}(\bm{w}_{h})\in{\bm{H}}^{1}_{0}(\Omega)$ and using (23), (26) and (24) leads to

[TABLE]

The elliptic regularity assumption implies $\|\bm{\psi}\|_{H^{1+s}(\Omega)}\leq C_{\mathrm{ell},s}\|\bm{E}_{h}\bm{w}_{h}-{\bm{S}}(\bm{w}_{h})\|_{L^{2}(\Omega)}$ and yields

[TABLE]

Finally, a triangle inequality gives

[TABLE]

This concludes the proof.

The previous result and similar arguments as in the conforming case enable us to prove the following theorem.

Theorem 8.3**.**

Suppose that the Stokes problem satisfies Assumption 2.1 and that the reconstruction operator ${\bm{I}}_{h}$ satisfies Assumption 4.1. Then there holds

[TABLE]

with $C_{1}>0$ and $C_{3}>0$ given by (15) and (25), respectively. Without Assumption 4.1, a result similar to Theorem 6.1 is valid, i.e.,

[TABLE]

Proof 8.4**.**

The proof of the first result is nearly identical to the proof of Theorem 7.1 with slight changes concerning the application of $\nabla_{h}$ and the replacement of Lemma 5.3 by Lemma 8.1. Likewise, the proof of the second result is almost identical to the proof of Theorem 6.1. However, one term has to be estimated differently, as follows. With $\mathrm{div}{\bm{S}}({\bm{e}}_{h})=0$ and (26), it holds

[TABLE]

9 Numerical Example

This sections gives a short numerical example to illustrate the theory. We consider an L-shaped domain $\Omega:=(-1,1)^{2}\setminus((0,1)\times(-1,0))$ and the manufactured solution

[TABLE]

where

[TABLE]

and $\gamma=856399/1572864\approx 0.54$, $\omega=3\pi/2$ taken from [17]. Note, that this yields $-\nu\Delta{\bm{v}}+\nabla p_{0}=0$. To have a nonzero right-hand side we add $p_{+}:=\sin(xy\pi)$ to the pressure, i.e. $p:=p_{0}+p_{+}$ and $\bm{f}:=\nabla(p_{+})$. Note that the exact solutions satisfy ${\bm{v}}\in{\bm{H}}^{1+s}(\Omega)$ and $p\in H^{s}(\Omega)$ for any $s<\gamma$. Moreover, we set the viscosity parameter to either $\nu=1$, $\nu=10^{-2}$ or $\nu=10^{-4}$.

Tables 1-6 compare the ${\bm{H}}^{1}$ errors of the classical Bernardi–Raugel finite element method and its pressure-robust sibling on a series of unstructured uniformly red-refined meshes for $\nu=1$ (Tables 1 and 2), $\nu=10^{-2}$ (Tables 3 and 4) and $\nu=10^{-4}$ (Tables 5 and 6). For the classical method the distance between the discrete Stokes projector and the discrete solution is non-zero and really scales with $\nu^{-1}$, but asymptotically converges with $h$ instead of $h^{s}$. At first glance this seems better than expected in Theorem 6.1, but the first term vanishes due to $\mathbb{P}(-\Delta{\bm{v}})=\nu^{-1}\mathbb{P}(\nabla p_{0})=0$ in this example. This also pre-asymptotically leads to a slightly higher convergence order of the full error than in case of $\nu=1$ at least for $\nu=10^{-2}$ and $\nu=10^{-4}$ where the $\mathcal{O}(h)$ error dominates at first. The numbers of the modified pressure-robust variant convey that the discrete solution of the modified method and the discrete Stokes projector are identical as predicted by Lemma 7.1 (again due to $\mathbb{P}(-\Delta{\bm{v}})=\nu^{-1}\mathbb{P}(\nabla p_{0})=0$). The numerical results confirm that for pressure-robust methods, the discrete velocity is independent of $\nu$. However, this $\nu$-independence only holds up to a quadrature error in the right-hand side, which scales with $\nu^{-1}$, and up to round-off errors.