Keywords: Generalized Newtonian fluid; finite element method; a priori error estimates; divergence reconstruction operator; inhomogeneous Dirichlet boundary condition.

AMS MSC (2020): 35J60; 35Q35; 65N12; 65N15; 65N30; 76A05.

1 Introduction↩︎

The steady motion of a homogeneous, incompressible generalized Newtonian fluid can be modeled by the generalized Navier–Stokes equations, i.e., \[\begin{align} \label{eq:fem:main95problem1} - \operatorname{div} \mathbf{S}(\mathbf{Dv}) + \operatorname{div}(\mathbf{v}\otimes\mathbf{v}) + \nabla q = \mathbf{f} \quad \textrm{ in } \Omega\,. \end{align}\tag{1}\] Here, \(\Omega\subseteq \mathbb{R}^d\), \(d\in \{2,3\}\), denotes either a bounded polygonal (if \(d=2\)) or polyhedral (if \({d=3}\)) Lipschitz domain and \(\mathbf{f}\colon \Omega\to \mathbb{R}^d\) denotes an external force. The velocity vector field \(\mathbf{v}\colon \overline{\Omega}\to \mathbb{R}^d\) and the scalar kinematic pressure \(q\colon \Omega\to \mathbb{R}\) are the unknowns in 1 . \(\mathbf{D}\) denotes the symmetric gradient operator. Moreover, we assume that the viscosity of the fluid can be expressed as a function of the shear rate. Therefore, we consider a non-degenerate extra-stress tensor \(\mathbf{S}\colon \mathbb{R}^{d\times d}\to \mathbb{R}_{\mathrm{sym}}^{d\times d} \, \mathrel{\vcenter{:}}= {\{\mathbf{A}\in \mathbb{R}^{d\times d}\mid \mathbf{A}=\mathbf{A}^\top\}}\) that has \((p,\delta)\)-structure (cf.Assumption 1), where \(\delta > 0\) and \(p\in (1,+\infty)\) is the shear stress exponent. A prototypical example for \(\mathbf{S}\colon \mathbb{R}^{d\times d}\to \mathbb{R}_{\mathrm{sym}}^{d\times d}\) is \[\begin{align} \mathbf{S}(\mathbf{A}) \mathrel{\vcenter{:}}= \nu_0\,(\delta+\vert \mathbf{A}^{\mathrm{sym}}\vert)^{p-2}\mathbf{A}^{\mathrm{sym}}\,, \end{align}\] for \(\mathbf{A}\in \mathbb{R}^{d\times d}\), where \(\nu_0>0\), \(\delta > 0\), \(p\in (1,+\infty)\) and \(\mathbf{A}^{\mathrm{sym}}\mathrel{\vcenter{:}}= \frac{1}{2}(\mathbf{A}+\mathbf{A}^\top)\in \mathbb{R}^{d\times d}_{\mathrm{sym}}\).

While the particular case \(p=2\) results in the well-known Navier–Stokes equations for Newtonian fluids, the cases \(p<2\) and \(p>2\) correspond to shear-thinning and shear-thickening fluids, respectively. Note that shear-thinning is the most common type of generalized Newtonian behavior of fluids and is seen in many industrial and everyday applications (cf.[1]). This motivates the study of finite element approximations of generalized Navier–Stokes equations 1 for small shear stress exponents \(p<2\).

The system 1 is completed by the inhomogeneous Dirichlet boundary condition \[\begin{align} {2} \mathbf{v} &= \mathbf{g}_2 &&\quad \textrm{ on } \partial\Omega \,,\tag{2} \\ \intertext{and, for the sake of generality, the \textit{inhomogeneous divergence constraint}} \operatorname{div}\mathbf{v} &= g_1&& \quad \textrm{ in } \Omega\,.\tag{3} \end{align}\]

This work is dedicated to finite element (FE) approximations of 1 –3 , with the focus on convergence rates for velocity and kinematic pressure.

1.1 Related contributions↩︎

A FE discretization of 1 and a proof of weak convergence to a solution of 1 is given in [2]. More precisely, in [2], homogeneous Dirichlet boundary conditions and the maximal range of the shear stress exponent (i.e., \(p\in (\frac{2d}{d+2},+\infty)\)) are considered: in the case \(p > \tfrac{2d}{d+1}\), the standard Temam modification (cf.[3]) of the weak convective term is admissible, while in the case \(p > \tfrac{2d}{d+2}\), exactly divergence-free FE pairs can be employed. A fully-discrete FE discretization of the unsteady Boussinesq system and a proof of weak convergence is given in [4]. There, in the case \({p \in (\tfrac{2d}{d+2}, \tfrac{2d}{d+1}]}\), in order to recover the crucial cancellation property of the weak convective term, the usage of a divergence reconstruction operator is proposed. On the theoretical side, this covers and generalizes the case of divergence-free FE pairs. On the practical side, the large computational costs needed by common divergence-free FE pairs are reduced in this way. To the best of the authors’ knowledge, the only contribution on a priori error estimates for a finite discretization of 1 is [5]. There, an inhomogeneous Dirichlet boundary condition and the restricted range \(p \geq \tfrac{2d}{d+1}\) (for which Temam’s modification (cf.[3]) is still admissible) is considered. The inhomogeneous Dirichlet boundary condition is imposed using the Fortin interpolation operator of the velocity FE space, which is (for many FE spaces; we refer the reader to [5] for a short list) given as a Scott–Zhang interpolation operator plus a divergence correction operator. The divergence correction operator does not need to be computed explicitly as it does not contribute to the Dirichlet boundary traces. Nevertheless, we still consider the computation of the Scott–Zhang interpolation operator as computationally expensive.

1.2 New contributions↩︎

The present paper is intended to improve the FE discretization in [5] in such a way that the above mentioned gaps are closed. Thus, the new contribution of the present paper is two-fold:

• First, we vary the discrete Dirichlet boundary condition in such a way that any approximation of the (continuous) Dirichlet boundary condition can be employed. In particular, if the Dirichlet boundary data are sufficiently regular, then this gives a theoretical justification for the usage of nodal interpolation, which is new in this context.

• Second, we use the method of divergence-free reconstruction operators in order to treat the full range of the shear stress exponent (i.e., \(p\in (\frac{2d}{d+2},+\infty)\)). To the best of the authors’ knowledge, this is the first a priori error analysis for a discretization of the generalized Navier–Stokes equations 1 in the case \(p \in (\tfrac{2d}{d+2}, \tfrac{2d}{d+1})\). More precisely, we derive a priori error estimates for the velocity vector field with error decay rates that are optimal for any \(p \in (\tfrac{2d}{d+2}, \tfrac{2d}{d+1})\) and a priori error estimates for the kinematic pressure with error decay rates that are optimal for any \(p \in (\tfrac{2d}{d+2},2]\). The derivation of a priori error estimates for the pressure turns out to be more challenging than in the case of a laminar flow (cf. [6]), i.e., 1 without the convective term. More precisely, from the a priori error estimates for the velocity vector field and a bootstrap argument, we deduce improved stability estimates for the discrete velocity vector field, which we then employ to prove the a priori error estimates for the kinematic pressure.

Since these two developments are independent of each other, we decided to present our results in a rather general way. More precisely, we also consider the Temam case, since this –in connection with nodal interpolation of Dirichlet data– is still new and probably the most relevant case for application. Therefore, certain steps of the proofs work out analogously to [5] and we resort to the respective assertions.

The paper is structured as follows: in Section 2, the necessary notation and a precise formulation of the continuous problem are given. In Section 3, assumptions on the finite element spaces are stated and the discretization is developed. Section 4 is dedicated to the proofs of a priori error estimates for the velocity vector field and the kinematic pressure. In Section 5, numerical experiments are presented and the theoretical findings are compared to the experimental results.

2 Preliminaries↩︎

2.1 Notation↩︎

Throughout the entire article, we denote by \(\Omega \subseteq \mathbb{R}^d\), \(d \in \{2,3\}\), either a bounded polygonal (if \(d=2\)) or a bounded polyhedral (if \(d=3\)) Lipschitz domain. We employ \(c>0\) to denote a generic constant, that may change from line to line, but does not depend on the crucial quantities. Moreover, we write \(u\sim v\) if and only if there exist a constant \(c>0\) such that \(c^{-1}\, u \leq v\leq c\, u\).

For \(k\in \mathbb{N}\) and \(p\in [1,\infty]\), we employ the customary Lebesgue spaces \({(L^p(\Omega), \|\cdot\|_{p,\Omega})}\) and Sobolev spaces \((W^{k,p}(\Omega),\|\cdot\|_{k,p,\Omega})\). There exists a linear, continuous, and surjective trace operator \(\tr_{\partial\Omega}\colon W^{1,p}(\Omega) \to W^{1-\smash{\frac{1}{p}},p}(\partial\Omega)\). If it gets clear from the context, we do not denote it explicitly for the sake of readability. The space \(\smash{W^{1,p}_{0}(\Omega)}\) is defined as those functions from \(W^{1,p}(\Omega)\) whose traces vanish on \(\partial\Omega\). The conjugate exponent is denoted by \(p' = \tfrac{p}{p-1} \in [1, \infty]\), whereas \(p^*\) denotes the critical Sobolev exponent with respect to the embedding \(W^{1,p}(\Omega) \hookrightarrow L^{q}(\Omega)\), i.e., if \({p > d}\), then \({p^* = \infty}\), while \(p^*\) is used as a placeholder for any finite \({q \in [1, \infty)}\) if \(p=d\). If \(p=d\), we formally understand \((p^*)'\) as a placeholder for \(1+\varepsilon\), where \(\varepsilon > 0\) is fixed, but arbitrary.

We always denote vector-valued functions by lowercase boldface letters and tensor-valued functions by capital boldface letters. The Euclidean inner product between two vectors \(\mathbf{a} =(a_1,\dots,a_d),\mathbf{b} =(b_1,\dots,b_d)\in \mathbb{R}^d\) is denoted by \({\mathbf{a} \cdot\mathbf{b}\mathrel{\vcenter{:}}= \sum_{i=1}^d{a_ib_i}}\), while the Frobenius inner product between two tensors \(\mathbf{A}=(A_{ij})_{i,j\in \{1,\dots,d\}},\mathbf{B}=(B_{ij})_{i,j\in \{1,\dots,d\}}\in \mathbb{R}^{d\times d }\) is denoted by \(\mathbf{A}: \mathbf{B}\mathrel{\vcenter{:}}= \sum_{i,j=1}^d{A_{ij}B_{ij}}\). For a tensor \(\mathbf{A}\in \mathbb{R}^{d\times d}\), we denote its symmetric part by \(\mathbf{A}^{\mathrm{sym}}\mathrel{\vcenter{:}}= \frac{1}{2}(\mathbf{A}+\mathbf{A}^\top)\in \mathbb{R}^{d\times d}_{\mathrm{sym}}\mathrel{\vcenter{:}}= {\{\mathbf{A}\in \mathbb{R}^{d\times d}\mid \mathbf{A}=\mathbf{A}^\top\}}\). Moreover, for a (Lebesgue) measurable set \(\omega\subseteq \mathbb{R}^d\) and (Lebesgue) measurable functions \(u,v\colon \omega\to \mathbb{R}\) (written \(u,v\in L^0(\omega)\)), we employ the product \[\begin{align} (u,v)_\omega \mathrel{\vcenter{:}}= \int_\omega u\, v\,\mathrm{d}x\,, \end{align}\] whenever the right-hand side is well-defined. The integral mean of an integrable function \(u\in L^1(\omega)\) over a (Lebesgue) measurable set \(\omega\subseteq \mathbb{R}^d\) with \(\vert \omega\vert>0\) is denoted by \[\begin{align} \langle u\rangle_\omega \mathrel{\vcenter{:}}= \frac{1}{|\omega|}\int_\omega u \,\mathrm{d}x\,. \end{align}\]

For \(r \in (1, \infty)\), we abbreviate the function spaces \[\begin{align} \begin{aligned} X^r &\mathrel{\vcenter{:}}= (W^{1,r}(\Omega))^d\,, &&\quad V^r \mathrel{\vcenter{:}}= (W_{0}^{1,r}(\Omega))^d\,, \\ Y^r &\mathrel{\vcenter{:}}= \smash{ L^{r'}(\Omega)}\,, &&\quad Q^r \mathrel{\vcenter{:}}= \smash{ L_0^{r'}(\Omega)} \mathrel{\vcenter{:}}= \big \{ z \in \smash{L^{r'}(\Omega) }\mid \langle z \rangle_{\Omega} = 0 \big\} \,, \end{aligned} \end{align}\] as well as \(X\mathrel{\vcenter{:}}= X^p\), \(V\mathrel{\vcenter{:}}= V^p\), \(Y\mathrel{\vcenter{:}}= Y^p\), and \(Q\mathrel{\vcenter{:}}= Q^p\).

2.2 N-functions and Orlicz spaces↩︎

A convex function \(\psi \colon \mathbb{R}_{\geq 0} \to \mathbb{R}_{\geq 0}\) is called N-function if \(\psi(0)=0\), \(\psi>0\) in \(\mathbb{R}_{> 0}\), \(\lim_{t\rightarrow0}{\frac{\psi(t)}{t}}=0\) and \({\lim_{t\rightarrow\infty}{\frac{\psi(t)}{t}}=\infty}\). The conjugate N-function \(\psi^*\colon \mathbb{R}_{\geq 0} \to \mathbb{R}_{\geq 0}\) is defined by \(\psi^*(s)\mathrel{\vcenter{:}}=\sup_{t \geq 0}\{s\,t - \psi(t)\}\) for all \(s \geq 0\). An N-function \({\psi\colon \mathbb{R}_{\geq 0} \to \mathbb{R}_{\geq 0}}\) satisfies the \(\Delta_2\)-condition, if there exists \(K> 2\) such that \({\psi(2\,t) \leq K\,\psi(t)}\) for all \(t \geq 0\). The smallest such constant is denoted by \(\Delta_2(\psi) > 0\). We need the following refined version of the \(\varepsilon\)-Young inequality: for every \(\varepsilon>0\), there exists a constant \(c_\varepsilon>0\), depending only on \(\Delta_2(\psi),\Delta_2( \psi ^*)<\infty\), such that for every \({s,t\geq 0}\), there holds \[\begin{align} \label{eq:fem:orliczyoung} s\,t&\leq c_\varepsilon \,\psi^*(s)+ \varepsilon \, \psi(t)\,. \end{align}\tag{4}\]

For \(p \in (1,\infty)\) and \(\delta\geq 0\), we define the special N-function \(\phi \mathrel{\vcenter{:}}= \phi_{p,\delta}\colon\smash{\mathbb{R}_{\geq 0}\to \mathbb{R}_{\geq 0}}\) by \[\begin{align} \label{eq:fem:nfunction} \phi(t)\mathrel{\vcenter{:}}= \int _0^t \phi'(s)\, \mathrm{d}s\,, \quad\text{where}\quad \phi'(t) \mathrel{\vcenter{:}}= (\delta +t)^{p-2} t\,,\quad\mathrm{ for all }t\geq 0\,. \end{align}\tag{5}\]

For an N-function \(\psi\colon\mathbb{R}_{\geq 0}\to \mathbb{R}_{\geq 0}\), we define shifted N-functions \({\psi_a\colon\mathbb{R}_{\geq 0}\to \mathbb{R}_{\geq 0}}\) via \[\begin{align} \label{eq:fem:phi95shifted} \psi_a(t)\mathrel{\vcenter{:}}= \int _0^t \psi_a'(s)\, \mathrm{d}s\,,\quad\text{where }\quad \psi'_a(t)\mathrel{\vcenter{:}}= \psi'(a+t)\frac{t}{a+t}\,,\quad\mathrm{ for all }a, t\geq 0\,. \end{align}\tag{6}\] The shifted special N-functions \(\phi_a\colon \mathbb{R}_{\geq 0}\to \mathbb{R}_{\geq 0}\), \(a\ge 0\), and their conjugate N-functions \((\phi_a)^*\colon \mathbb{R}_{\geq 0}\to \mathbb{R}_{\geq 0}\), \(a\ge 0\), satisfy the \(\Delta_2\)-condition with \[\begin{align} \sup_{a\ge 0}{\Delta_2(\phi_a)}&\leq c\,2^{\max \{2,p\}}\,,\\ \sup_{a\ge 0}{\Delta_2((\phi_a)^*)}&\leq c\,2^{\max \{2,p'\}}\,. \end{align}\] In addition, uniformly in \(t,\delta\geq 0\), we have that \[\begin{align} \phi_a(t) &\sim (\delta + a + t)^{p-2} t^2\,,\tag{7}\\ (\phi_a)^*(t)& \sim ((\delta+a)^{p-1} + t)^{p'-2} t^2\tag{8}\,. \end{align}\]

Let \(\omega\subseteq \mathbb{R}^d\), \(d \in \mathbb{N}\), be a (Lebesgue) measurable set. A Carathéodory function \(\psi \colon \omega \times \mathbb{R}_{\geq 0} \to \mathbb{R}_{\geq 0}\), such that \(\psi(x,\cdot)\) for a.e. \(x \in \omega\) is an N-function, is called generalized N-function. For a given generalized N-function \(\psi \colon \omega \times \mathbb{R}_{\geq 0} \to \mathbb{R}_{\geq 0}\), the modular (with respect to \(\psi\) on \(\omega\)) of a (Lebesgue) measurable function \(u\in L^0(\omega)\) is defined by \[\begin{align} \rho_{\psi,\omega}(u)\mathrel{\vcenter{:}}= \int_\omega {\psi(\cdot,\vert u\vert)\,\mathrm{d}x}\,. \end{align}\]

2.3 The extra-stress tensor and related functions↩︎

In this subsection, we specify the assumptions on the extra-stress tensor and important consequences thereof.

Since we use estimates that are only available for \(\delta>0\) in Section 4, we exclude the degenerate case \(\delta=0\) in Assumption 1. Closely related to the extra-stress tensor \(\mathbf{S}\colon \mathbb{R}^{d \times d} \to \mathbb{R}^{d \times d}_{\mathrm{sym}}\) with \((p,\delta)\)-structure (cf.Assumption 1) is the non-linear mapping \(\mathbf{F} \colon\mathbb{R}^{d\times d}\to \mathbb{R}^{d\times d}_{\mathrm{sym}}\), for every \(\mathbf{A}\in \mathbb{R}^{d\times d}\) defined by \[\begin{align} \begin{aligned} \mathbf{F}(\mathbf{A})&\mathrel{\vcenter{:}}= (\delta+\vert \mathbf{A}^{\mathrm{sym}}\vert)^{\smash{\frac{p-2}{2}}}\mathbf{A}^{\mathrm{sym}}\,. \end{aligned} \label{eq:def95F} \end{align}\tag{9}\] The non-linear mappings \(\mathbf{S},\mathbf{F}\colon \mathbb{R}^{d \times d} \to \mathbb{R}^{d\times d}_{\mathrm{sym}}\) and \(\phi_a,(\phi_a)^*\colon \mathbb{R}^{\ge 0}\to \mathbb{R}^{\ge 0}\), \({a\ge 0}\), are closely related.

2.4 Continuous weak formulation↩︎

In this subsection, we give a precise weak formulation of the strong formulation 1 –3 .

The convective term is consistently reformulated as \[\begin{align} \label{eq:fem:conv95reform} (\operatorname{div}(\mathbf{v} \otimes \mathbf{v}), \mathbf{z})_{\Omega} = - (\mathbf{v} \otimes \mathbf{v}, \nabla \mathbf{z})_\Omega \,. \end{align}\tag{10}\] Admissibility of the right hand side of 10 will be ensured by the restriction \(p > \tfrac{2d}{d+2}\) and \(\mathbf{z} \in V^s\).

Then, a first weak formulation of the strong formulation 1 –3 is given via:

Problem (Q). Under the Assumption 2, find \((\mathbf{v}, q) \in X \times Q^{s}\) with \(\mathbf{v} = \mathbf{g}_2\) a.e.on \(\partial\Omega\) such that for every \((\mathbf{z},z)\in V^{s} \times Y^{s}\), there holds \[\begin{align} (\mathbf{S}(\mathbf{Dv}), \mathbf{Dz})_\Omega - (\mathbf{v} \otimes \mathbf{v}, \nabla \mathbf{z})_\Omega - (q, \operatorname{div}\mathbf{z})_\Omega &= (\mathbf{f}, \mathbf{z})_\Omega \,, \tag{11} \\ (\operatorname{div}\mathbf{v},z)_\Omega &= (g_1,z)_\Omega \,. \tag{12} \end{align}\]

We temporarily assume that there exists a divergence and trace lift \(\mathbf{g} \in X^s\), which solves the system \[\begin{align} \begin{aligned} \label{eq:fem:div95eq} \operatorname{div}\mathbf{g} &= g_1 &&\quad\text{ in }L^s(\Omega)\,, \\[-0.5mm] \mathbf{g} &= \mathbf{g}_2 &&\quad\text{ in }(W^{\smash{1-\frac{1}{s}},s}(\partial\Omega))^d\,. \end{aligned} \end{align}\tag{13}\] An explicit construction of such \(\mathbf{g}\) is given later in Lemma 9. Then, setting \[\mathbf{u} \mathrel{\vcenter{:}}= \mathbf{v} - \mathbf{g} \in V_0\,,\] where \(V_0 \mathrel{\vcenter{:}}= V_0^p\) and for \(r\in [1,\infty)\) \[\begin{align} V_0^r \mathrel{\vcenter{:}}= \big \{ \mathbf{z} \in V^r \mid\operatorname{div}\mathbf{z} = 0 \text{ a.e.\;in }\Omega\big\}\,, \end{align}\] we obtain the following weak formulation:

Problem (P). Under the Assumption 2, find \(\mathbf{u} \in V_0\) such that for every \(\mathbf{z}\in V_0^{s}\), there holds \[\begin{align} (\mathbf{S}(\mathbf{Du}+\mathbf{Dg}), \mathbf{Dz})_\Omega - ((\mathbf{u}+\mathbf{g}) \otimes (\mathbf{u}+\mathbf{g}), \nabla \mathbf{z})_\Omega = (\mathbf{f}, \mathbf{z})_\Omega\,. \end{align}\]

First, the existence of the velocity vector field in Problem (P) can be guaranteed via the celebrated Lipschitz truncation method (cf.[8]). Second, the existence of the kinematic pressure in Problem (Q) can be inferred on the basis of the following inf-sup stability result. This, in turn, implies equivalence of Problem (Q) and Problem (P).

3 Discretization↩︎

3.1 Triangulations and finite element spaces↩︎

The following assumption is needed in the construction of a discrete (divergence-corrected) Lipschitz truncation (cf.[2] or [11]).

Assumption 6 implies the following discrete inf-sup stability result.

3.2 Discrete Dirichlet boundary condition↩︎

We assume that among the degrees of freedom parametrizing \(X_h\), those which are relevant for the Dirichlet boundary condition 2 can be uniquely determined (cf. [10]). We interpret the Dirichlet boundary condition 2 as Dirichlet boundary condition at the Dirichlet degrees of freedom in \(X_h\). As a consequence, an interpolation of boundary data \(\mathbf{g}_b \in X^s\) needs to be computed. The Fortin interpolation operator has ideal analytical properties for that purpose (cf. [5] and Remark 3), but this comes with significant computational costs. Therefore, we develop a setting which allows for interpolation operators which are computationally inexpensive and straightforward to implement. Since regularity of boundary data is rarely critical in applications, the nodal interpolation operator is probably the most interesting example. In this subsection, we formulate requirements on continuous and discrete boundary data. Then, we discuss the applicability of some interpolation operators.

The following assumptions will be sufficient for existence and convergence of discrete solutions:

The regularity \(\mathbf{g}_b \in X^s\) complies with the best known assumptions for existence of weak solutions to Problem (Q) (at least for small values of \(p\)), cf.[14].

Higher regularity of \(\mathbf{g}_b\in X^s\) and approximability by \(\mathbf{g}_b^h\in X_h\) will be employed in the proof of a priori error estimates in Section 4:

Note that Assumption 8 implies Assumption 7 (cf. reasoning in [15]). While \(k=2\), \(\beta=p\) in Assumption 8 is certainly a natural choice in order to establish linear convergence rates, different regularities may be imposed due to the admissibility of interpolation operators and technical restrictions in the analysis.

The next Lemma shows that despite the nodal interpolation operator \(\Pi_{h, N}^X\) is not \(W^{1,1}(\Omega)\)-stable, it still satisfies Assumption 8 (under appropriate regularity assumptions).

We note that the difference \(\Pi_h^X \mathbf{g}_b - \mathbf{g}_b^h\) approaches zero with an at least linear rate in \(h\): due to Hölder’s inequality, Sobolev approximability (cf.[21]), and ?? , for every \(t \in [1, \infty)\) and \(\smash{\widetilde{\beta}} \mathrel{\vcenter{:}}= \max \{t, \beta\}\), there holds \[\begin{align} \begin{aligned} \label{eq:fem:ph-gb} \|\Pi_h^X \mathbf{g}_b - \mathbf{g}_b^h\|_{1,t,\Omega} &\leq c \, \|\Pi_h^X \mathbf{g}_b - \mathbf{g}_b\|_{1,\smash{\widetilde{\beta}},\Omega} + c \, \|\mathbf{g}_b - \mathbf{g}_b^h\|_{1,\smash{\widetilde{\beta}},\Omega} \\&\leq c\, h\, \|\nabla^k \mathbf{g}_b\|_{\smash{\widetilde{\beta}},\Omega} \,. \end{aligned} \end{align}\tag{14}\]

3.3 Discrete convective term↩︎

Classically, skew-symmetry plays a crucial role in the analysis of the convective term. While this property seems out of reach in the case of a non-zero divergence constraint on the velocity, we may at least preserve certain cancellation properties. Therefore, in the case \(p \geq \smash{\tfrac{2d}{d+1}}\), we employ discretizations \({b,\widetilde{b}\colon (X)^2\times X^s \to\mathbb{R}}\), commonly referred to as Temam’s modifications, for every \({(\mathbf{u},\mathbf{v},\mathbf{w})\in (X)^2\times X^s}\), (cf.[3]) defined by \[\begin{align} \label{eq:fem:defi95conv95temam} \begin{aligned} b(\mathbf{u},\mathbf{v},\mathbf{w}) &\mathrel{\vcenter{:}}= \tilde{b}(\mathbf{u},\mathbf{v},\mathbf{w}) + \tfrac{1}{2} (g_1 \mathbf{u}, \mathbf{w})_\Omega \\ &\mathrel{\vcenter{:}}= \tfrac{1}{2} (\mathbf{w}\otimes\mathbf{u}, \nabla \mathbf{v} + g_1 \mathbf{I})_\Omega - \tfrac{1}{2} (\mathbf{v}\otimes\mathbf{u}, \nabla \mathbf{w})_\Omega\,. \end{aligned} \end{align}\tag{15}\]

The lower bound on the range of \(p\) is clearly suboptimal, as the (continuous) existence theory for weak solutions allows for values \(p > \tfrac{2d}{d+2}\) (cf.[8]). The underlying issue can be circumvented by the usage of a (divergence) reconstruction operator, which maps discretely to exactly divergence-free functions.

Next, we present short lists of common mixed finite element spaces \(\{X_h\}_{h>0}\) and \(\{Y_h\}_{h>0}\) with projectors \(\{\Pi_h^X\}_{h>0}\) and \(\{\Pi_h^Y\}_{h>0}\) on regular triangulations \(\{\mathcal{T}_h\}_{h>0}\) satisfying both Assumption 5 and Assumption 6, for which reconstruction operators \(\{\Sigma_h\}_{h>0}\) in the sense of Definition 2 are available.

In the second case in Assumption 9, the discrete convective term, for every \(\smash{(\mathbf{u},\mathbf{v},\mathbf{w})}\in (X)^2\times X\), is defined by \[\begin{align} \label{eq:fem:defi95conv95recon} b(\mathbf{u},\mathbf{v},\mathbf{w}) \mathrel{\vcenter{:}}= - (\mathbf{v}\otimes \Sigma_h \mathbf{u}, \nabla \mathbf{w})_{\Omega}\,. \end{align}\tag{19}\]

3.4 Discrete problem formulation↩︎

Since the boundary conditions fix volume flux everywhere on the topological boundary \(\partial\Omega\), the combination of discretized boundary and divergence data have to fulfill a compatibility condition. Therefore, we define \[\begin{align} \label{def:g95195tilde} g_1^h \mathrel{\vcenter{:}}= g_1 + \langle\mathrm{div}\,\mathbf{g}_b^h - g_1\rangle_{\Omega} \in \smash{Y^{s'}} \,. \end{align}\tag{20}\]

We consider the following discrete counterpart of Problem (Q):

Problem (Q\(_h\)). Find \((\mathbf{v}_h, q_h) \in X_h \times Q_h\) with \(\mathbf{v}_h=\mathbf{g}_b^h\) in \(\tr X_h\) such that for every \((\mathbf{z}_h,y_h) \in V_h\times Y_h\), there holds \[\begin{align} (\mathbf{S}(\mathbf{Dv}_h), \mathbf{Dz}_h)_{\Omega} + b(\mathbf{v}_h, \mathbf{v}_h,\mathbf{z}_h ) - (q_h, \mathrm{div}\,\mathbf{z}_h )_{\Omega} &= (\mathbf{f}, \mathbf{z}_h )_{\Omega}\,, \tag{21} \\ (\mathrm{div}\,\mathbf{v}_h,y_h)_{\Omega} &= (g_1^h, y_h)_{\Omega}\,. \tag{22} \end{align}\]

In order to establish the existence of discrete solutions to Problem (Q\(_h\)), a discrete counterpart of Problem (P) is useful. As a preparation, we construct a discrete extension of the boundary data which has the correct divergence.

In order to establish the existence of a discrete velocity vector field in Problem (Q\(_h\)), we prefer a problem formulation “hiding” the discrete pressure. Taking the discrete lift \(\mathbf{g}_h\in X_h\) from Lemma 9 and making the ansatz \[\begin{align} \label{eq:discrete95ansatz} \mathbf{u}_h \mathrel{\vcenter{:}}= \mathbf{v}_h - \mathbf{g}_h \in V_{h,0}\,,\\[-6.5mm]\notag \end{align}\tag{23}\] leads to the following discrete counterpart of Problem (P):

Problem (P\(_h\)). Given a solution \(\mathbf{g}_h \in X_h\) of ?? , ?? , find \(\mathbf{u}_h \in V_{h,0}\) such that for every \(\mathbf{z}_h \in V_{h,0}\), there holds \[\begin{align} (\mathbf{S}(\mathbf{Du}_h+\mathbf{Dg}_h), \mathbf{Dz}_h)_\Omega + b(\mathbf{u}_h+\mathbf{g}_h, \mathbf{u}_h+\mathbf{g}_h,\mathbf{z}_h) = (\mathbf{f}, \mathbf{z}_h)_\Omega\,.\\[-6.5mm]\notag \end{align}\]

By Lemma 9 and Lemma 2, Problem (Q\(_h\)) and Problem (P\(_h\)) are equivalent.

3.5 Well-posedness and convergence analysis↩︎

Analogously to [5], we establish the well-posedness (i.e., existence of discrete solutions), stability (i.e., a priori bounds), and (weak) convergence of Problem (Q_\(h\)) (and Problem (P_\(h\)), respectively). While in [5], only the case \(p \geq \tfrac{2d}{d+1}\) is treated, one easily checks that the argumentation works out indeed completely analogously, as it is the case in the continuous existence theory. Due to the length of the proofs, we refrain from repeating them here.

4 A priori error estimates↩︎

In this section, we derive a priori error estimates for the approximation of a “regular” solution \((\mathbf{v}, q)\in X\times Q^s\) of Problem (Q) by a solution \((\mathbf{v}_h, q_h)\in X_h\times Q_h\) of Problem (Q\(_h\)).

Therefore, we define the auxiliary coefficients \[\begin{align} r & \mathrel{\vcenter{:}}= \min \{2, p\}\,,\\ \ell & \mathrel{\vcenter{:}}= \max \{2, p, s \} = \left\{ \begin{aligned} &s = \smash{(\tfrac{p^*}{2})'} &&\textrm{if } p < \tfrac{4d}{d+4}\,, \\ &2 &&\textrm{if } p \in [\tfrac{4d}{d+4}, 2]\,, \\ &p &&\textrm{else\,.} \end{aligned} \right. \end{align}\]

The next step is to bound the pressure error up to the velocity error.

The proof of Proposition 12 follows the same reasoning as in [5] using the following error estimate for the convective term.

If \(\{\mathcal{T}_h\}_{h>0}\) is quasi-uniform, a stronger stability statement follows from inverse estimates and the regularity of the velocity vector field (cf.Assumption 10).

5 Numerical experiments↩︎

In this section, we review the theoretical findings of Section 4 via numerical experiments. Since in [5], the case \(p\ge \frac{2d}{d+1}\) (i.e., \(p\ge \frac{4}{3}\) if \(d=2\)) is already studied, the numerical experiments focus on the case \({p<2}\).

5.1 Implementation details↩︎

All experiments were conducted employing the finite element software firedrake (version 0.13.0, cf.[36]). In the numerical experiments, we deploy the finite element spaces in Remark 5 with a discontinuous pressure space. In the case \(p<\frac{2d}{d+1}\), the discrete convective term is given via 19 and involves a (divergence) reconstruction operator \(\Sigma_h\colon X\to X_h+Z_h\) in the sense of Definition 2. We incorporated the latter not via constructing it explicitly, but weakly instead, i.e., via extending Problem (Q\(_h\)) with an unknown \(\mathbf{z}_h\in Z_h\) (serving as a placeholder for \(\Sigma_h \mathbf{v}_h\in Z_h\)) in the following ways:

\(\bullet\) P2P0 element/First-order Bernardi–Raugel element. According to Remark 5(i), a divergence reconstruction operator is given by the Fortin interpolation operator \({\Sigma_h\mathrel{\vcenter{:}}= \Pi_h^{rt,0}\colon}\) \(X\to Z_h\) (cf.[10]) of the lowest order Raviart–Thomas element \(Z_h\mathrel{\vcenter{:}}= \mathcal{R}T^0(\mathcal{T}_h)\). We incorporated this interpolation operator into Problem (Q\(_h\)) via extending this non-linear saddle point problem on the basis of the definition ?? . More precisely, we solve the following equivalent augmented problem:

Augmented Problem (AQ\(_h^0\)). Find \((\mathbf{v}_h,\mathbf{z}_h,q_h)\in X_h\times Z_h\times Q_h\) such that \(\mathbf{v}_h = \mathbf{g}_b^h\) in \(\tr X_h\) and for every \((\mathbf{w}_h,\mathbf{y}_h,y_h)\in X_h\times Z_h\times Q_h\), there holds \[\begin{align} (\mathbf{S}(\mathbf{Dv}_h), \mathbf{Dw}_h)_{\Omega} - (\mathbf{v}_h\otimes \mathbf{z}_h, \nabla \mathbf{w}_h)_{\Omega} - (q_h, \operatorname{div} \mathbf{w}_h )_{\Omega} &= (\mathbf{f}, \mathbf{w}_h )_{\Omega} \,, \\ \langle \mathbf{z}_h\cdot \mathbf{n}, \mathbf{y}_h\cdot \mathbf{n}\rangle_{\mathcal{F}_h}&=\langle \mathbf{v}_h\cdot \mathbf{n}, \mathbf{y}_h\cdot \mathbf{n}\rangle_{\mathcal{F}_h}\,,\\ (\operatorname{div} \mathbf{v}_h, y_h)_{\Omega} &= (g_1^h, y_h)_{\Omega}\,, \end{align}\] where we exploit that \((\mathbf{y}_h\mapsto \mathbf{y}_h\cdot \mathbf{n}|_F)\colon Z_h\to \mathbb{P}^0(F)\) for all \(F\in \mathcal{F}_h\) is surjective and we used the notation \(\langle \mathbf{z}_h\cdot \mathbf{n}, \mathbf{y}_h\cdot \mathbf{n}\rangle_{\mathcal{F}_h}\mathrel{\vcenter{:}}=\sum_{F\in \mathcal{F}_h}{\int_F{(\mathbf{z}_h\cdot \mathbf{n})\,(\mathbf{y}_h\cdot \mathbf{n})\,\mathrm{d}s}}\).

\(\bullet\) Conforming Crouzeix–Raviart element/Second-order Bernardi–Raugel element. According to Remark 5(ii), a divergence reconstruction operator is given via the Fortin interpolation operator \(\Sigma_h\mathrel{\vcenter{:}}= \Pi_h^{rt,1} \colon X\to Z_h\) (cf.[10]) of the first degree Raviart–Thomas element \(Z_h\mathrel{\vcenter{:}}= \mathcal{R}T^1(\mathcal{T}_h)\). We incorporated this operator into Problem (Q\(_h\)) via extending this non-linear saddle point problem based on the definitions ?? , ?? . More precisely, we solve the following equivalent augmented problem:

Augmented Problem (AQ\(_h^1\)). Find \((\mathbf{v}_h,\mathbf{z}_h,q_h) \in X_h\times Z_h\times Q_h\) such that \({\mathbf{v}_h = \mathbf{g}_b^h}\) in \(\tr X_h\) and for every \((\mathbf{w}_h,\mathbf{y}_h^1,\mathbf{y}_h^0,y_h) \in V_h\times Z_h|_{\cup\mathcal{F}_h}\times (\mathbb{P}^0(\mathcal{T}_h))^d \times Y_h\), there holds \[\begin{align} (\mathbf{S}(\mathbf{Dv}_h),\mathbf{Dw}_h)_{\Omega} - (\mathbf{v}_h\otimes \mathbf{z}_h, \nabla \mathbf{w}_h)_{\Omega} - (q_h, \operatorname{div} \mathbf{w}_h )_{\Omega} &= (\mathbf{f}, \mathbf{w}_h )_{\Omega} \,, \\ \langle \mathbf{z}_h\cdot \mathbf{n}, \mathbf{y}_h^1\cdot \mathbf{n}\rangle_{\mathcal{F}_h}&=\langle \mathbf{v}_h\cdot \mathbf{n}, \mathbf{y}_h^1\cdot \mathbf{n}\rangle_{\mathcal{F}_h}\,,\\ (\mathbf{z}_h, \mathbf{y}_h^0)_{\Omega}&=(\mathbf{v}_h, \mathbf{y}_h^0)_{\Omega}\,,\\ (\operatorname{div} \mathbf{v}_h, y_h)_{\Omega} &= (g_1^h, y_h)_{\Omega}\,, \end{align}\] where we exploit that \((\mathbf{y}_h\mapsto \mathbf{y}_h\cdot \mathbf{n}|_F)\colon Z_h|_F\to \mathbb{P}^1(F)\) for all \(F\in \mathcal{F}_h\) is surjective.

We emphasize the Augmented Problems (AQ\(_h^0\)) and (AQ\(_h^1\)) are designed in such a way that both the formulas ?? , ?? (where ?? is only relevant for Augmented Problem (AQ\(_h^1\))) for the reconstruction operator and the structural connection between the discrete pressure space \(Y_h\) and the discrete velocity space \(X_h\) are guaranteed, where \(\Sigma_h\colon X_h \to X_h + Z_h\) with \(\operatorname{div}(Z_h) = Y_h\) (cf. Remark 5).

As discretization of the boundary data, we employ the nodal interpolation of the trace lift \({\mathbf{g}_b\in X^s}\), i.e., we set \(\mathbf{g}_b^h\mathrel{\vcenter{:}}= \Pi_{h,N}^X\mathbf{g}_b\in X_h\), which, according to Lemma 4, given sufficient regularity of the trace lift \(\mathbf{g}_b\in X^s\), is justified. In addition, since we will restrict to the case \({g_1=0}\), we set \({g_1^h\mathrel{\vcenter{:}}= \langle \operatorname{div} \mathbf{g}_b^h\rangle_{\Omega}\in Y^{s'}}\).

We approximate the discrete solutions \((\mathbf{v}_h,\mathbf{z}_h,q_h)\in X_h\times Z_h\times Q_h\) of the Augmented Problems (AQ\(_h^0\)) and (AQ\(_h^1\)), respectively, using a standard Newton iteration which is deemed to have converged when the Euclidean norm of the residual falls below the tolerance \(\tau_{abs} \mathrel{\vcenter{:}}= 1\textrm{e}{-}8\). The linear system emerging in each Newton iteration is solved using a sparse direct solver from MUMPS (version 5.5.0, cf. [37]).

5.2 Experimental setup↩︎

We use the Augmented Problems (AQ\(_h^0\)) and (AQ\(_h^1\)), respectively, (or equivalently the Problems (Q\(_h\)),(P\(_h\))) to approximate the non-linear system 1 –3 on \(\Omega=(0,1)^2\) with \(\mathbf{S}\colon \mathbb{R}^{2\times 2}\to \mathbb{R}^{2\times 2}_{\mathrm{sym}}\), for every \(\mathbf{A}\in\mathbb{R}^{2\times 2}\) defined by \[\begin{align} \mathbf{S}(\mathbf{A}) \mathrel{\vcenter{:}}= \nu_0\,(\delta+\vert \mathbf{A}^{\mathrm{sym}}\vert)^{p-2}\mathbf{A}^{\mathrm{sym}}\,, \end{align}\] where \(\nu_0 =100\), \(\delta\mathrel{\vcenter{:}}=1\textrm{e}{-}5\), and \(p\in (1, 1.5]\).

As manufactured solutions serve the vector field \(\mathbf{v}\in X\) and the function \(q \in Q^s\), for every \(x=(x_1,x_2)\in \Omega\), defined by \[\begin{align} \mathbf{v}(x)\mathrel{\vcenter{:}}= \vert x\vert^\beta(-x_2,x_1)\,, \qquad q(x)\vcentcolon = \vert x\vert^{\gamma}-\langle\,\vert \cdot\vert^{\gamma}\,\rangle_\Omega\,, \end{align}\] i.e., we choose the right-hand side \(\mathbf{f}\in \smash{(L^{p'}(\Omega))^2}\), the divergence \(g_1=0\in L^s(\Omega)\), and boundary data \(\mathbf{g}_2\in \textcolor{black}{(W^{\smash{1-\frac{1}{s},s}}(\partial\Omega))^2}\) accordingly.

For the regularity of the velocity vector field and kinematic pressure, we choose \(\beta = 0.01\) and \(\gamma= 1-\smash{\frac{2}{p'}}+\beta\), which just yields that Assumption 10 is satisfied by \(\mathbf{v}\) and \(q\). The boundary data \(\mathbf{v}|_{ \partial\Omega}\) may, in general, not be regular enough to admit a lift \(\mathbf{g}_b \in (W^{2,\ell}(\Omega))^2\), nevertheless, the EOC presented in Tables [tab:1]–[tab:3] reaches the theoretically proven rates. This gives another indication that the extra regularity of the boundary data (cf.Assumption 10) is probably not necessary.

An initial triangulation \(\mathcal{T}_{h_0}\), where \(h_0=1\), is constructed via subdividing \(\overline{\Omega}\) into four triangles along its diagonals. Then, refined triangulations \(\mathcal{T}_{h_i}\), \(i=1,\ldots,7\), where \({h_{i+1}=\frac{h_i}{2}}\) for all \({i=0,\ldots,6}\), are obtained by applying uniform refinement; more precisely, the red-refinement rule (cf. [35]).

As estimation of the convergence rates, the experimental order of convergence (EOC) \[\begin{align} \texttt{EOC}_i(e_i)\mathrel{\vcenter{:}}=\frac{\log(e_{i+1})-\log(e_i)}{\log(h_{i+1})-\log(h_i)}\,, \quad i=1,\dots,6\,,\label{eoc} \end{align}\tag{50}\] is used, where, for every \(i=1,\ldots,7\), we denote by \(e_i\) a general error quantity.

5.3 Quasi-optimality of the error decay rates derived in Theorem 13↩︎

For \(p \in \{1.1,1.2,1.3, 4/3, 1.4, 1.5\}\), the first-order Bernardi–Raugel element (cf.Remark 5(i.a)) and the conforming Crouzeix–Raviart element (cf.Remark 5(ii.a)), for \(i=1,\ldots,7\), we compute the error quantities \[\begin{align} \smash{e_{\mathbf{v},i}^{\mathbf{F}}}\mathrel{\vcenter{:}}=\|\mathbf{F}(\mathbf{D}\mathbf{v}_{h_i})-\mathbf{F}(\mathbf{D}\mathbf{v})\|_{2,\Omega}\,, \end{align}\] and the corresponding EOCs, which are presented in Table [tab:1]:

For both the first-order Bernardi–Raugel element and the conforming Crouzeix–Raviart element, we report the expected convergence rate of \(\texttt{EOC}_i(\smash{e_{\mathbf{v},i}^{\mathbf{F}}}) \approx 1\), \(i=1,\dots,6\), which confirms the quasi-optimality of the a priori error estimates for the velocity vector field derived in Theorem 13. For the P2P0 element (cf.Remark 5(i.b)), we observed the same error decay rates as for the first-order Bernardi–Raugel element.

[H]
    \setlength\tabcolsep{4.5pt}
    \centering
    \begin{tabular}{c |c|c|c|c|c|c|c|c|c|c|c|c|c|c|} 
        \cmidrule[\heavyrulewidth](){2-13}
        & \multicolumn{6}{c||}{\cellcolor{lightgray}first-order Bernardi--Raugel}   & \multicolumn{6}{c|}{\cellcolor{lightgray}conforming Crouzeix--Raviart}\\ 
        \cmidrule(){2-13}
        & \multicolumn{3}{c|}{\cellcolor{lightgray}using \eqref{eq:fem:defi95conv95recon}} & \multicolumn{3}{c||}{\cellcolor{lightgray}using \eqref{eq:fem:defi95conv95temam}}
        & \multicolumn{3}{c|}{\cellcolor{lightgray}using \eqref{eq:fem:defi95conv95recon}} & \multicolumn{3}{c|}{\cellcolor{lightgray}using \eqref{eq:fem:defi95conv95temam}}\\
        \hline 
        \multicolumn{1}{|c||}{\cellcolor{lightgray}\diagbox[height=1.1\line,width=0.11\dimexpr\textwidth]{\vspace{-0.6mm}$i$}{\\[-5mm] $p$}}
        & \cellcolor{lightgray}1.1 & \cellcolor{lightgray}1.2 & \cellcolor{lightgray}1.3 & \cellcolor{lightgray}4/3  & \cellcolor{lightgray}1.4  &  \multicolumn{1}{c||}{\cellcolor{lightgray}1.5} &
    \multicolumn{1}{c|}{\cellcolor{lightgray}1.1}  & \cellcolor{lightgray}1.2 &\cellcolor{lightgray}1.3 & \cellcolor{lightgray}4/3 &\cellcolor{lightgray}1.4 & \cellcolor{lightgray}1.5   \\ \toprule\toprule
    \multicolumn{1}{|c||}{\cellcolor{lightgray}$1$}            & 1.012 & 1.011 & 1.009 & 1.009 & 1.008 & \multicolumn{1}{c||}{1.008} & \multicolumn{1}{c|}{1.001} & 1.002 & 1.002 & 1.002 & 1.002 & 1.002 \\ \hline
    \multicolumn{1}{|c||}{\cellcolor{lightgray}$2$}            & 1.010 & 1.010 & 1.010 & 1.010 & 1.010 & \multicolumn{1}{c||}{1.010} & \multicolumn{1}{c|}{1.009} & 1.010 & 1.010 & 1.010 & 1.010 & 1.010 \\ \hline
    \multicolumn{1}{|c||}{\cellcolor{lightgray}$3$}            & 1.008 & 1.008 & 1.009 & 1.009 & 1.009 & \multicolumn{1}{c||}{1.009} & \multicolumn{1}{c|}{1.007} & 1.007 & 1.008 & 1.008 & 1.008 & 1.008 \\ \hline
    \multicolumn{1}{|c||}{\cellcolor{lightgray}$4$}            & 1.007 & 1.007 & 1.008 & 1.008 & 1.008 & \multicolumn{1}{c||}{1.008} & \multicolumn{1}{c|}{1.006} & 1.006 & 1.007 & 1.007 & 1.007 & 1.008 \\ \hline
    \multicolumn{1}{|c||}{\cellcolor{lightgray}$5$}            & 1.006 & 1.007 & 1.007 & 1.007 & 1.007 & \multicolumn{1}{c||}{1.008} & \multicolumn{1}{c|}{1.006} & 1.006 & 1.007 & 1.007 & 1.007 & 1.008 \\ \hline
    \multicolumn{1}{|c||}{\cellcolor{lightgray}$6$}            & 1.006 & 1.006 & 1.007 & 1.007 & 1.007 & \multicolumn{1}{c||}{1.008} & \multicolumn{1}{c|}{1.006} & 1.006 & 1.007 & 1.007 & 1.007 & 1.008 \\ \toprule\toprule
    \multicolumn{1}{|c||}{\cellcolor{lightgray}\small theory}  & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 & \multicolumn{1}{c||}{1.000} & \multicolumn{1}{c|}{1.000} & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 \\ \toprule
    \end{tabular}\vspace{-2mm}
    \caption{Experimental order of convergence: $\texttt{EOC}_i(\smash{e_{\mathbf{v},i}^{\mathbf{F}}})$, ${i=1,\dots,6}$.}
    \label{tab:1}

5.4 Quasi-optimality of the error decay rates proved in Proposition 14↩︎

For \(p \in \{1.1,1.2,1.3, 4/3, 1.4, 1.5\}\), the first-order Bernardi–Raugel element (cf.Remark 5(i.a)) and the conforming Crouzeix–Raviart element (cf.Remark 5(ii.a)), we compute the error quantities \[\begin{align} \left.\begin{aligned} \smash{e_{q,i}^{L^{p'}}}&\mathrel{\vcenter{:}}= \|q_{h_i}-q\|_{p',\Omega}\,,\\ \smash{e_{q,i}^{L^2}}&\mathrel{\vcenter{:}}= \|q_{h_i}-q\|_{2,\Omega}\,, \end{aligned}\quad\right\}\quad i=0,\dots,6\,, \end{align}\] and the corresponding EOCs, which are presented in Table [tab:2] and Table [tab:3], respectively:

For the first-order Bernardi–Raugel element and the conforming Crouzeix–Raviart element, we report the expected convergence rates of \(\texttt{EOC}_i(\smash{e_{q,i}^{L^{p'}}}) \approx \frac{2}{p'}\), \(i=1,\dots,6\), and \(\texttt{EOC}_i(\smash{e_{q,i}^{L^2}}) \approx 1\), \({i=1,\dots,6}\), which confirms the quasi-optimality of the a priori error estimates for the pressure in Proposition 14. For the P2P0 element (cf. Remark 5(i.b)), we observed the same error decay rates as for the first-order Bernardi–Raugel element.

[H]
    \setlength\tabcolsep{4.5pt}
    \centering
    \begin{tabular}{c |c|c|c|c|c|c|c|c|c|c|c|c|c|c|} 
        \cmidrule[\heavyrulewidth](){2-13}
        & \multicolumn{6}{c||}{\cellcolor{lightgray}first-order Bernardi--Raugel}   & \multicolumn{6}{c|}{\cellcolor{lightgray}conforming Crouzeix--Raviart}\\ 
        \cmidrule(){2-13}
        & \multicolumn{3}{c|}{\cellcolor{lightgray}using \eqref{eq:fem:defi95conv95recon}} & \multicolumn{3}{c||}{\cellcolor{lightgray}using \eqref{eq:fem:defi95conv95temam}}
        & \multicolumn{3}{c|}{\cellcolor{lightgray}using \eqref{eq:fem:defi95conv95recon}} & \multicolumn{3}{c|}{\cellcolor{lightgray}using \eqref{eq:fem:defi95conv95temam}}\\
        \hline 
        \multicolumn{1}{|c||}{\cellcolor{lightgray}\diagbox[height=1.1\line,width=0.11\dimexpr\textwidth]{\vspace{-0.6mm}$i$}{\\[-5mm] $p$}}
        & \cellcolor{lightgray}1.1 & \cellcolor{lightgray}1.2 & \cellcolor{lightgray}1.3 & \cellcolor{lightgray}4/3  & \cellcolor{lightgray}1.4  &  \multicolumn{1}{c||}{\cellcolor{lightgray}1.5} &
        \multicolumn{1}{c|}{\cellcolor{lightgray}1.1}  & \cellcolor{lightgray}1.2 &\cellcolor{lightgray}1.3 & \cellcolor{lightgray}4/3 &\cellcolor{lightgray}1.4 & \cellcolor{lightgray}1.5   \\ \toprule\toprule
        \multicolumn{1}{|c||}{\cellcolor{lightgray}$1$}            & 0.351 & 0.057 & 0.010 & -0.02 & 0.061 & \multicolumn{1}{c||}{0.260} & \multicolumn{1}{c|}{0.052} & 0.227 & 0.387 & 0.436 & 0.527 & 0.644 \\ \hline
        \multicolumn{1}{|c||}{\cellcolor{lightgray}$2$}            & 0.602 & 0.558 & 0.403 & 0.432 & 0.524 & \multicolumn{1}{c||}{0.668} & \multicolumn{1}{c|}{0.179} & 0.328 & 0.457 & 0.496 & 0.570 & 0.671 \\ \hline
        \multicolumn{1}{|c||}{\cellcolor{lightgray}$3$}            & 0.321 & 0.282 & 0.422 & 0.470 & 0.558 & \multicolumn{1}{c||}{0.676} & \multicolumn{1}{c|}{0.182} & 0.333 & 0.462 & 0.501 & 0.573 & 0.671 \\ \hline
        \multicolumn{1}{|c||}{\cellcolor{lightgray}$4$}            & 0.224 & 0.315 & 0.456 & 0.497 & 0.574 & \multicolumn{1}{c||}{0.678} & \multicolumn{1}{c|}{0.183} & 0.334 & 0.464 & 0.503 & 0.575 & 0.671 \\ \hline
        \multicolumn{1}{|c||}{\cellcolor{lightgray}$5$}            & 0.195 & 0.330 & 0.463 & 0.503 & 0.577 & \multicolumn{1}{c||}{0.677} & \multicolumn{1}{c|}{0.183} & 0.335 & 0.464 & 0.503 & 0.575 & 0.672 \\ \hline
        \multicolumn{1}{|c||}{\cellcolor{lightgray}$6$}            & 0.187 & 0.334 & 0.465 & 0.504 & 0.577 & \multicolumn{1}{c||}{0.675} & \multicolumn{1}{c|}{0.183} & 0.335 & 0.464 & 0.503 & 0.575 & 0.672 \\ \toprule\toprule
        \multicolumn{1}{|c||}{\cellcolor{lightgray}\small theory}  & 0.181 & 0.333 & 0.462 & 0.500 & 0.571 & \multicolumn{1}{c||}{0.667} & \multicolumn{1}{c|}{0.181} & 0.333 & 0.462 & 0.500 & 0.571 & 0.667 \\ \toprule
        \end{tabular}\vspace{-2mm}
    \caption{Experimental order of convergence: $\texttt{EOC}_i(\smash{e_{q,i}^{L^{p'}}})$, ${i=1,\dots,6}$.}
    \label{tab:2}

[H]
    \setlength\tabcolsep{4.5pt}
    \centering
    \begin{tabular}{c |c|c|c|c|c|c|c|c|c|c|c|c|c|c|} 
        \cmidrule[\heavyrulewidth](){2-13}
        & \multicolumn{6}{c||}{\cellcolor{lightgray}first-order Bernardi--Raugel}   & \multicolumn{6}{c|}{\cellcolor{lightgray}conforming Crouzeix--Raviart}\\ 
        \cmidrule(){2-13}
        & \multicolumn{3}{c|}{\cellcolor{lightgray}using \eqref{eq:fem:defi95conv95recon}} & \multicolumn{3}{c||}{\cellcolor{lightgray}using \eqref{eq:fem:defi95conv95temam}}
        & \multicolumn{3}{c|}{\cellcolor{lightgray}using \eqref{eq:fem:defi95conv95recon}} & \multicolumn{3}{c|}{\cellcolor{lightgray}using \eqref{eq:fem:defi95conv95temam}}\\
        \hline 
        \multicolumn{1}{|c||}{\cellcolor{lightgray}\diagbox[height=1.1\line,width=0.11\dimexpr\textwidth]{\vspace{-0.6mm}$i$}{\\[-5mm] $p$}}
        & \cellcolor{lightgray}1.1 & \cellcolor{lightgray}1.2 & \cellcolor{lightgray}1.3 & \cellcolor{lightgray}4/3  & \cellcolor{lightgray}1.4  &  \multicolumn{1}{c||}{\cellcolor{lightgray}1.5} &
        \multicolumn{1}{c|}{\cellcolor{lightgray}1.1}  & \cellcolor{lightgray}1.2 &\cellcolor{lightgray}1.3 & \cellcolor{lightgray}4/3 &\cellcolor{lightgray}1.4 & \cellcolor{lightgray}1.5   \\ \toprule\toprule
        \multicolumn{1}{|c||}{\cellcolor{lightgray}$1$}            & 0.571 & 0.092 & -0.02 & 0.045 & 0.191 & \multicolumn{1}{c||}{0.364} & \multicolumn{1}{c|}{0.899} & 0.905 & 0.918 & 0.923 & 0.933 & 0.946 \\ \hline
        \multicolumn{1}{|c||}{\cellcolor{lightgray}$2$}            & 1.197 & 1.092 & 1.016 & 1.005 & 0.998 & \multicolumn{1}{c||}{1.007} & \multicolumn{1}{c|}{1.016} & 1.011 & 1.009 & 1.009 & 1.010 & 1.011 \\ \hline
        \multicolumn{1}{|c||}{\cellcolor{lightgray}$3$}            & 1.129 & 1.047 & 1.013 & 1.009 & 1.009 & \multicolumn{1}{c||}{1.014} & \multicolumn{1}{c|}{1.002} & 1.000 & 1.001 & 1.001 & 1.002 & 1.004 \\ \hline
        \multicolumn{1}{|c||}{\cellcolor{lightgray}$4$}            & 1.059 & 1.019 & 1.007 & 1.006 & 1.006 & \multicolumn{1}{c||}{1.008} & \multicolumn{1}{c|}{1.001} & 1.000 & 1.001 & 1.001 & 1.002 & 1.003 \\ \hline
        \multicolumn{1}{|c||}{\cellcolor{lightgray}$5$}            & 1.029 & 1.008 & 1.004 & 1.003 & 1.004 & \multicolumn{1}{c||}{1.005} & \multicolumn{1}{c|}{1.001} & 1.001 & 1.002 & 1.002 & 1.003 & 1.004 \\ \hline
        \multicolumn{1}{|c||}{\cellcolor{lightgray}$6$}            & 1.015 & 1.004 & 1.003 & 1.003 & 1.003 & \multicolumn{1}{c||}{1.004} & \multicolumn{1}{c|}{1.001} & 1.002 & 1.002 & 1.003 & 1.003 & 1.004 \\ \toprule\toprule
        \multicolumn{1}{|c||}{\cellcolor{lightgray}\small theory}  & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 & \multicolumn{1}{c||}{1.000} & \multicolumn{1}{c|}{1.000} & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 \\ \toprule
    \end{tabular}\vspace{-2mm}
    \caption{Experimental order of convergence: $\texttt{EOC}_i(\smash{e_{q,i}^{L^2}})$, ${i=1,\dots,6}$.}
    \label{tab:3}

References↩︎

---

1. Email: kaltenbach@math.tu-berlin.de↩︎

2. Email: julius.jessberger@mathematik.uni-freiburg.de↩︎

3. \(H(\operatorname{div};\Omega)\mathrel{\vcenter{:}}=\{\mathbf{z}\in (L^2(\Omega))^d\mid \operatorname{div}\mathbf{z}\in L^2(\Omega)\}\).↩︎

4. For a set \(M\subseteq \mathbb{R}^d\), by \(\mathrm{dim}_{\mathscr{H}}(M)\mathrel{\vcenter{:}}= \inf\{s>0\mid \mathscr{H}^s(M)=0\}\) we denote the Hausdorff dimension.↩︎