category: literaturenote citekey: badianonaugmentedvelocityvorticity2025 title: Non-augmented velocity-vorticity-pressure formulation for the Navier--Stokes--Brinkman--Forchheimer problem authors: "Badia, Santiago; Carstensen, Carsten; Martin, Alberto F.; Ruiz-Baier, Ricardo; Villa-Fuentes, Segundo" year: 2025 date: 2025-00-00 2025 url: "https://arxiv.org/abs/2506.10533v1" zotero_key: JRICDH58 zotero_storage: WP7XUC7N collections: transpordi_planeerimine folder: 001_artiklid firstAuthor: "Badia, Santiago"

status: converted

Non-augmented velocity-vorticity-pressure formulation for the Navier–Stokes–Brinkman–Forchheimer problem^∗

Santiago Badia† Carsten Carstensen‡ Alberto F. Mart´ın§ Ricardo Ruiz-Baier¶ Segundo Villa-Fuentes‖

June 13, 2025

Abstract

The flow of incompressible fluid in highly permeable porous media in vorticity - velocity - Bernoulli pressure form leads to a double saddle-point problem in the Navier–Stokes–Brinkman–Forchheimer equations. The paper establishes, for small sources, the existence of solutions on the continuous and discrete level of lowest-order piecewise divergence-free Crouzeix–Raviart finite elements. The vorticity employs a vector version of the pressure space with normal and tangential velocity jump penalisation terms. A simple Raviart–Thomas interpolant leads to pressure-robust a priori error estimates. An explicit residual-based a posteriori error estimate allows for efficient and reliable a posteriori error control. The efficiency for the Forchheimer nonlinearity requires a novel discrete inequality of independent interest. The implementation is based upon a light-weight forest-of-trees data structure handled by a highly parallel set of adaptive mesh refining algorithms. Numerical simulations reveal robustness of the a posteriori error estimates and improved convergence rates by adaptive mesh-refining.

Keywords: Navier–Stokes–Brinkman–Forchheimer equations; pressure robustness; nonconforming finite element methods; Banach fixed-point theory; a priori and a posteriori error estimates.

Mathematics Subject Classifications (2000): 65N15, 65N30, 76D05, 76M10.

1 Introduction

Scope. Equations of incompressible flow in dimensionless rotational or Lamb form (using vorticity) are of high importance in a number of applications. See, for example, the following non-exhaustive list of recent contributions analysing numerical methods based on different formulations [2, 14, 20, 32, 38, 46, 47, 52] (and see also the references therein). It was observed in [5, 4, 6] that, in order to control

^∗This work has been supported by Monash Mathematics through a Gordon Preston Sabbatical Fellowship; by the Australian Research Council through the Future Fellowship grant FT220100496 and Discovery Project grant DP22010316; and by the National Research and Development Agency (ANID) of the Ministry of Science, Technology, Knowledge and Innovation of Chile through the postdoctoral program Becas Chile grant 74220026. Computational resources were provided by the Australian Government through NCI under the National Computational Merit Allocation Scheme (NCMAS) and the ANU Merit Allocation Scheme (ANUMAS).

^†School of Mathematics, Monash University, 9 Rainforest Walk, Melbourne 3800 VIC, Australia. Email: santiago.badia@monash.edu.

Institut f¨ur Mathematik, Humboldt-Universit¨at zu Berlin, 10117 Berlin, Germany. Email: cc@math.hu-berlin.de. §School of Computing, Australian National University, Acton ACT 2601, Australia. Email: alberto.f.martin@anu.edu.au.

^¶School of Mathematics, Monash University, 9 Rainforest Walk, Melbourne 3800 VIC, Australia. Email: ricardo.ruizbaier@monash.edu.

^‖School of Mathematics, Monash University, 9 Rainforest Walk, Melbourne 3800 VIC, Australia. Present address: Departamento de Matem´atica, Universidad del B´ıo-B´ıo, Concepci´on, Chile. Email: svilla@ubiobio.cl.

the full H¹ (Ω) norm of the velocity and maintain optimal convergence, vorticity-based formulations for incompressible flow (with vorticity sought in L 2 (Ω) instead of the more common H(curl , Ω) case, velocity in H(curl , Ω) ∩ H(div, Ω), and Bernoulli pressure in L² (Ω)) required augmentation leastsquares terms coming from the incompressibility condition and constitutive equation for vorticity (the latter resembling also the vorticity-stabilisation from, e.g., [1, 13]). The aforementioned works [5, 4, 6] (which tackle Oseen, Navier–Stokes, and Forchheimer equations) also contain numerical evidence that either the vorticity or divergence stabilisation parameters could be zero for some finite element pairs approximating velocity and Bernoulli pressure.

The general Navier–Stokes and Forchheimer nonlinearities as well as the Brinkman drag will be included in this paper using fixed-point arguments. In contrast to the works above, we do not use augmentation techniques and treat the problem as a (perturbation of a) perturbed saddle-point problem embedded in another saddle-point problem. The analysis hinges on working on the kernel of the divergence operator. At the discrete level we can perform a fairly similar analysis as long as we use kernel-characterising spaces, such that the divergence of the discrete velocity is zero locally in each cell. For this we can take for instance the nonconforming Crouzeix–Raviart finite element pair [28], which is Stokes inf-sup stable and satisfies the required local kernel property. Pressure robust discretisations achieve velocity errors in the broken H¹ (Ω) norm which are proportional to the best approximation error in the velocity, without dependence on the velocity error. In order to do so in the present setting, we include a modification in the discrete right-hand side functional and in the nonlinear variational forms using a lowest-order H(div, Ω)-conforming interpolate of the velocity test function. This approach has been used for Stokes, Navier–Stokes, and many other variants in [1, 15, 40, 41, 43, 51], see also the numerous references therein. This property is also closely related to the L 2 (Ω) orthogonality of divergence-free functions onto the space of gradients of functions in H¹ (Ω). This approach comes at the price a small consistency error of optimal order, which is independent of the kinematic viscosity.

It is important to mention that the nonconformity of the method in combination with the need to control the H(curl , Ω) ∩ H(div, Ω) part of the velocity norm, imply that we need to stabilise the velocity-velocity bilinear form with tangential and normal jump terms to control the divergence and curl part of the velocity norm on the discrete level. This was done in, e.g., [17, 16] for the grad-div, curl-curl, and reduced Maxwell problems, and we recall that such stabilisation is not required for, e.g., Stokes equations with velocity in H¹ (Ω).

In many problems, solution singularities can cause suboptimal convergence, and adaptive mesh refinement is essential for recovering optimal rates. A number of variants of residual-based a posteriori estimators are available for Crouzeix–Raviart schemes applied to Stokes equations [25, 23]. In addition, the literature of pressure-robust methods also has works designing a posteriori error estimators [41, 39], including the use of divergence-free reconstruction operators and techniques that are commonly encountered in stream function-vorticity formulations of incompressible flow. The explicit residualbased a posteriori error estimate follows the overall frame of [24] (that includes historical remarks as well) with several additions for robustness (e.g. Lemma 5.3) required for the discrete norms. The efficiency result, however, was surprisingly subtle for some novel inverse estimate of independent innovation, we therefore establish in an appendix. The a posteriori error estimator contains a residual contribution and a non-conformity contribution. We show that it is reliable and efficient, where the efficiency proof relies upon a novel inverse estimate associated with element bubble functions. We use the estimator as an indicator for adaptive mesh refinement, and this restores optimal convergence in the case of singular solutions (which under uniform mesh refinement yield suboptimal convergence).

The particular adaptive meshes that we leverage in this work are hierarchically-adapted (i.e., nested) non-conforming octree-based meshes endowed with Morton (a.k.a., Z-shaped) space-fillingcurves for storage and data partitioning; see, e.g., [8]. This family of meshes can be very efficiently handled (i.e., refined, coarsened, re-partitioned, etc.) using high-performance and low-memory footprint algorithms [18]. While these are n-cube meshes (e.g., made of quadrilaterals or cubes in 2D and 3D, resp.), we split their elements (e.g., into 2 triangles or 6 tetrahedra, resp.) to obtain the simplicial meshes required by our finite element formulation. However, as these meshes are non-conforming (in particular they have hanging faces at cell interfaces between cells located at different levels of refinement), one needs to add additional multi-point constraints to the Crouzeix–Raviart finite element spaces in order to have optimal approximability properties. In particular, following [11], these additional constraints impose the trace average on a parent coarse face to be equivalent to the average of the trace averages on the children faces, and we show via numerical experiments that this approach recovers optimal convergence in the case of singular solutions.

Main contributions. In summary, to the best of our knowledge, the combination of the contributions addressed above (residual a posteriori error estimators for Navier–Stokes–Brinkman–Forchheimer equations in vorticity form using non-conforming methods) is novel. In particular this work features

a new non-augmented vorticity-based weak formulation for the Navier–Stokes equations with Brinkman and Forchheimer effects, generalising the recent work [6],
a rigorous solvability analysis and discrete problem for vorticity that attains pressure-robustness, complementing the works [40, 43],
a continuous and discrete analysis valid in 3D, extending the similar works for linear curl-curl type problems that address the 2D case [17, 16, 12],
a relatively simple residual-based a posteriori estimator for the pressure-robust scheme (compared to those from, e.g., [25, 39]),
the efficiency of the residual a posteriori error estimator requires an interesting novel result regarding inverse estimates for the Forchheimer nonlinearity,
efficient and reliable a posteriori estimators, which we prove theoretically and also confirmed numerically (noting that the effectivity index –the ratio between the total error and the global a posteriori error estimator– remains bounded between 1.6 and 1.9 for the tested cases),
new handling of multipoint constraints needed for implementation of Crouzeix–Raviart elements with hanging nodes in adaptive meshes constructed with octrees. This generalises the results from [8].

Outline. This article is organised as follows: in the remainder of this section we provide notational conventions and main assumptions to use throughout the paper. In Section 2 we give a brief overview on the governing equations and their statement in weak perturbed saddle-point form. Section 3 deals with the well-posedness analysis using Banach fixed-point theorem under small data assumptions and a global inf-sup argument. The definition of the discrete problem and the analysis of its unique solvability are addressed in Section 4. In Section 5 we derive C´ea estimates and error bounds for the specific finite element subspaces mentioned above. The definition of a posteriori error estimators and their robustness analysis is given in Section 6. We continue in Section 7 describing the benchmark setups we used in the numerical experiments, showcasing the properties of the proposed schemes, and confirming numerically the predicted a priori convergence estimates and robustness of the a posteriori error estimators. In Section 8 we give some concluding remarks, and Appendix A provides an inverse estimate for the efficiency proof of its own interest.

Preliminaries and notation. Let us denote by Ω ⊂ R ³ a bounded polyhedral Lipschitz domain with boundary Γ. Standard notations will be adopted for Lebesgue spaces L^p (Ω), with p ∈ [1, ∞] and Sobolev spaces W r,p(Ω) with r ≥ 0, endowed with the norms ∥ • ∥Lp(Ω) and ∥ • ∥^W r,p(Ω). Note that W ⁰,p(Ω) = L^p (Ω) and if p = 2, we write H^r (Ω) in place of W r,² (Ω), with the corresponding Lebesgue and Sobolev norms denoted by ∥ • ∥0,^Ω and ∥ • ∥r,Ω. The notation L² 0 (Ω) means the space of functions in L² (Ω) with zero mean value over Ω. We also write | • |r,^Ω for the H^r -seminorm. The bracket ⟨•, •⟩^Γ denotes duality that extends the L² (Γ) scalar product for smooth functions in the trace space H1/² (Γ) of H¹ (Ω) and its dual H−1/² (Γ). By S we will denote the corresponding vectorial counterpart of the generic scalar functional space S. The gradient, symmetric gradient, divergence and curl of a generic vector field v = (vi) are defined as

$$\nabla \boldsymbol{v} := \left(\frac{\partial v_i}{\partial xj}\right){i,j=1,3}, \quad \mathbf{D} \boldsymbol{v} := \frac{1}{2}(\boldsymbol{\nabla} \boldsymbol{v} + \boldsymbol{\nabla} \boldsymbol{v}^T), \quad \mathrm{div} \, \boldsymbol{v} := \sum_{j=1}^3 \frac{\partial v_j}{\partial x_j}, \quad \mathrm{and} \quad \mathbf{curl} \, \boldsymbol{v} := \boldsymbol{\nabla} \times \lambda \, \boldsymbol{v}.$$

In addition, we recall that the spaces

$$\begin{aligned} \mathbf{H}0(\mathrm{div};\Omega) &:= \big{ \boldsymbol{v} \in \mathbf{L}^2(\Omega) : & \mathrm{div}\, \boldsymbol{v} \in \mathbf{L}^2(\Omega) & \mathrm{and} & \gamma{\mathbf{n}} = 0 & \mathrm{on} & \partial \Omega \big}, \ \mathbf{H}0(\mathbf{curl};\Omega) &:= \big{ \boldsymbol{v} \in \mathbf{L}^2(\Omega) : & \mathbf{curl}\, \boldsymbol{v} \in \mathbf{L}^2(\Omega) & \mathrm{and} & \gamma{\boldsymbol{\tau}} = 0 & \mathrm{on} & \partial \Omega \big}, \end{aligned}$$

where γⁿ and γ^τ represent the normal and tangential trace, respectively, and Hilbert when equipped with the norms ∥v∥ 2 div,Ω := ∥v∥ 2 ⁰,^Ω + ∥div v∥ 2 0,Ω and ∥v∥ 2 curl ,Ω := ∥v∥ 2 ⁰,^Ω + ∥curl v∥ 2 0,Ω , respectively. Then, we define the following space

$$\mathbf{V} := \mathbf{H}_0(\operatorname{div}; \Omega) \cap \mathbf{H}_0(\mathbf{curl}; \Omega),$$

endowed with the norm

$$|v|{\mathbf{V}}^2 := |v|{0,\Omega}^2 + |\operatorname{div} v|{0,\Omega}^2 + |\operatorname{\mathbf{curl}} v|{0,\Omega}^2.$$

Finally, the notation A ≲ B abbreviates A ≤ CB with a generic h (mesh size)-independent constant C, while some of the constants are still written explicit to emphasise and quantify particular assumptions.

2 Model problem and its weak formulation

2.1 The governing equations

We start with the steady Navier–Stokes–Brinkman–Forchheimer equations in their usual velocity– pressure form. They consist in finding velocity u and kinematic pressure P such that

$$\kappa^{-1} \boldsymbol{u} - \nu \Delta \boldsymbol{u} + (\boldsymbol{u} \cdot \nabla) \boldsymbol{u} + \mathsf{F} |\boldsymbol{u}| \, \boldsymbol{u} + \frac{1}{\rho} \nabla P = \boldsymbol{f} \quad \text{in} \quad \Omega, \qquad \text{div } \boldsymbol{u} = 0 \quad \text{in} \quad \Omega$$ (2.1)

with κ > 0 the permeability of the porous media (assumed constant), ν the kinematic fluid viscosity, F > 0 the Forchheimer coefficient, ρ the fluid mass density, and f a given external force. Problem (2.1) can be equivalently set in terms of vorticity, velocity and pressure (similarly as done in, e.g., [7] for Brinkman and in [3] for Oseen equations). For this, we introduce the rescaled vorticity vector

$$\boldsymbol{\omega} := \sqrt{\nu} \operatorname{\mathbf{curl}} \boldsymbol{u}$$

and use the identity u · ∇u = curlu × u + ²∇(u · u). Then we introduce the rescaled Bernoulli pressure

$$p := \frac{1}{\rho}P + \frac{1}{2}\boldsymbol{u} \cdot \boldsymbol{u} - \lambda$$

for λ defined as the mean value of ¹ 2 u · u; and employ the following vector identity

$$\operatorname{\mathbf{curl}} \operatorname{\mathbf{curl}} \boldsymbol{u} = -\boldsymbol{\Delta} \boldsymbol{u} + \nabla (\operatorname{div} \boldsymbol{u})$$

together with the incompressibility constraint. These steps lead to the following equations

$$\kappa^{-1} \boldsymbol{u} + \sqrt{\nu} \operatorname{curl} \boldsymbol{\omega} + \mathsf{F} |\boldsymbol{u}| \, \boldsymbol{u} + \nabla p + \frac{1}{\sqrt{\nu}} \, \boldsymbol{\omega} \times \boldsymbol{u} = \boldsymbol{f}, \qquad \boldsymbol{\omega} - \sqrt{\nu} \operatorname{curl} \boldsymbol{u} = \boldsymbol{0}, \qquad \operatorname{div} \boldsymbol{u} = 0. \quad (2.2)$$

Furthermore, we focus on homogeneous Dirichlet boundary conditions for velocity and therefore an additional condition is required to enforce the uniqueness of the Bernoulli pressure. This gives

$$u = 0$$ on $\Gamma$ and $\int_{\Omega} p = 0$ . (2.3)

However, similar results as those shown below are also valid for other types of boundary conditions.

2.2 Mixed weak formulation

First, multiplying the first, second and third equations of (2.2) by v ∈ V, θ ∈ L 2 (Ω) and q ∈ L 2 0 (Ω), respectively, integrating by parts and utilising the boundary condition, we obtain the problem: Find ((u, ω), p) ∈ [V × L 2 (Ω)] × L 0 (Ω) such that

$$\kappa^{-1} \int{\Omega} \boldsymbol{u} \cdot \boldsymbol{v} + \sqrt{\nu} \int{\Omega} \boldsymbol{\omega} \cdot \operatorname{curl} \boldsymbol{v} - \int{\Omega} p \operatorname{div} \boldsymbol{v} - \frac{1}{\sqrt{\nu}} \int{\Omega} (\boldsymbol{u} \times \boldsymbol{\omega}) \cdot \boldsymbol{v} + \mathsf{F} \int{\Omega} |\boldsymbol{u}| \, \boldsymbol{u} \cdot \boldsymbol{v} = \int{\Omega} \boldsymbol{f} \cdot \boldsymbol{v},$$

$$\sqrt{\nu} \int{\Omega} \boldsymbol{\theta} \cdot \operatorname{curl} \boldsymbol{u} - \int{\Omega} \boldsymbol{\omega} \cdot \boldsymbol{\theta} = 0,$$

$$- \int_{\Omega} q \operatorname{div} \boldsymbol{u} = 0.$$ (2.4)

We introduce the bounded bilinear forms a : [V×L 2 (Ω)]×[V×L 2 (Ω)] → R, b : [V×L 2 (Ω)]×L 2 (Ω) → ^R, and, for each ^u^b ^∈ ^V, the bilinear form ^c ub : [V × L 2 (Ω)] × [V × L 2 (Ω)] → R as

$$a(\boldsymbol{u}, \boldsymbol{\omega}; \boldsymbol{v}, \boldsymbol{\theta}) := \kappa^{-1} \int{\Omega} \boldsymbol{u} \cdot \boldsymbol{v} + \sqrt{\nu} \int{\Omega} \boldsymbol{\omega} \cdot \operatorname{\mathbf{curl}} \boldsymbol{v} + \sqrt{\nu} \int{\Omega} \boldsymbol{\theta} \cdot \operatorname{\mathbf{curl}} \boldsymbol{u} - \int{\Omega} \boldsymbol{\omega} \cdot \boldsymbol{\theta}, \qquad (2.5a)$$

$$b(\boldsymbol{v},\boldsymbol{\theta};q) := -\int_{\Omega} q \operatorname{div} \boldsymbol{v}, \tag{2.5b}$$

$$c^{\widehat{\boldsymbol{u}}}(\boldsymbol{u},\boldsymbol{\omega};\boldsymbol{v},\boldsymbol{\theta}) := -\frac{1}{\sqrt{\nu}} \int{\Omega} (\widehat{\boldsymbol{u}} \times \boldsymbol{\omega}) \cdot \boldsymbol{v} + \mathsf{F} \int{\Omega} |\widehat{\boldsymbol{u}}| \, \boldsymbol{u} \cdot \boldsymbol{v}. \tag{2.5c}$$

On the other hand, we define the functional F ∈ [V × L 2 (Ω)]′ as

$$F(\boldsymbol{v}, \boldsymbol{\theta}) := \int_{\Omega} \boldsymbol{f} \cdot \boldsymbol{v}. \tag{2.6}$$

Then, the formulation consists in finding ((u, ω), p) ∈ [V × L 2 (Ω)] × L 2 0 (Ω), such that:

$$a(\boldsymbol{u}, \boldsymbol{\omega}; \boldsymbol{v}, \boldsymbol{\theta}) + b(\boldsymbol{v}, \boldsymbol{\theta}; p) + c^{\boldsymbol{u}}(\boldsymbol{u}, \boldsymbol{\omega}; \boldsymbol{v}, \boldsymbol{\theta}) = F(\boldsymbol{v}, \boldsymbol{\theta}),$$

$$b(\boldsymbol{u}, \boldsymbol{\omega}; q) = 0,$$ (2.7)

for all ((v, θ), q) ∈ [V × L (Ω)] × L 0 (Ω).

3 Analysis of the coupled problem

The following well-known symmetric and non-symmetric versions of the generalised Lax–Milgram lemma will be used in the forthcoming analysis (for a proof see, e.g., [30, Theorems 1.3 & 1.2]).

Lemma 3.1 Let H be a real Hilbert space, and let A : H×H → R be a symmetric and bounded bilinear form. Assume that

$$\sup{0 \neq v \in \mathcal{H}} \frac{A(u, v)}{|v|{\mathcal{H}}} \ge \widetilde{\alpha} |u|_{\mathcal{H}} \quad \forall u \in \mathcal{H}. \tag{3.1}$$

Then, for each F ∈ H′ there exists a unique u ∈ H such that

$$A(u,v) = F(v)$$ $\forall v \in \mathcal{H}$ , and $||u||{\mathcal{H}} \leq \frac{1}{\widetilde{\alpha}} ||F||{\mathcal{H}'}$ .

Lemma 3.2 Let H1, H² be real Hilbert spaces, and let B : H¹ × H² → R be a bounded bilinear form. Assume that

$$\sup_{0 \neq v \in \mathcal{H}2} \frac{B(u,v)}{|v|{\mathcal{H}2}} \, \geq \, \widehat{\alpha} \, |u|{\mathcal{H}_1} \quad \forall \, u \in \mathcal{H}_1, \tag{3.2a}$$

$$\sup_{u \in \mathcal{H}_1} B(u, v) > 0 \quad \forall v \in \mathcal{H}_2, v \neq 0.$$ (3.2b)

Then, for each F ∈ H′ 2 there exists a unique u ∈ H¹ such that

$$B(u,v) = F(v)$$ $\forall v \in \mathcal{H}2$ , and $|u|{\mathcal{H}1} \le \frac{1}{\widehat{\alpha}} |F|{\mathcal{H}'_2}$ .

We will combine these results with the Banach fixed-point theorem to demonstrate the wellposedness of (2.7) under a small data assumption.

3.1 Stability properties of a linear problem

First we recall the continuous embedding from H¹ (Ω) into L^p (Ω), for all p ∈ [1, 6]:

$$||w||{\mathrm{L}^{\mathbf{p}}(\Omega)} \le C{\mathrm{S}} ||w||_{1,\Omega} \quad \forall w \in \mathrm{H}^{1}(\Omega)$$

$$(3.3)$$

with C^S > 0 depending only on |Ω| and p (see [45, Theorem 1.3.4]).

Next we easily deduce from from the Cauchy–Schwarz inequality, the continuity of a(•, •), b(•, •):

$$\left| a(\boldsymbol{u}, \boldsymbol{\omega}; \boldsymbol{v}, \boldsymbol{\theta}) \right| \le (\kappa^{-1} + \sqrt{\nu} + 1) (|\boldsymbol{u}|{\mathbf{V}} + |\boldsymbol{\omega}|{0,\Omega}) (|\boldsymbol{v}|{\mathbf{V}} + |\boldsymbol{\theta}|{0,\Omega}), \tag{3.4a}$$

$$\left|b(\boldsymbol{v},\boldsymbol{\theta};q)\right| \le |q|{0,\Omega}(|\boldsymbol{v}|{\mathbf{V}} + |\boldsymbol{\theta}|_{0,\Omega}). \tag{3.4b}$$

In turn, using H¨older's inequality together with (3.3), we readily deduce that

$$\left| c^{\widehat{\boldsymbol{u}}}(\boldsymbol{u}, \boldsymbol{\omega}; \boldsymbol{v}, \boldsymbol{\theta}) \right| \le C{\mathrm{S}}^{2} \left( \frac{1}{\sqrt{\nu}} + \mathsf{F} \right) |\widehat{\boldsymbol{u}}|{\mathbf{V}} (|\boldsymbol{u}|{\mathbf{V}} + |\boldsymbol{\omega}|{0,\Omega}) (|\boldsymbol{v}|{\mathbf{V}} + |\boldsymbol{\theta}|{0,\Omega}). \tag{3.5}$$

Similarly, the linear functional F(•) is bounded

$$|F(\boldsymbol{v},\boldsymbol{\theta})| \le ||\boldsymbol{f}||{0,\Omega}(||\boldsymbol{v}||{\mathbf{V}} + ||\boldsymbol{\theta}||_{0,\Omega}). \tag{3.6}$$

Now, it is straightforward to see that the kernel of the bilinear form b(•, •) is a closed subspace of V × L 2 (Ω). It is denoted as V⁰ × L 2 (Ω), and the first component admits the characterisation

$$\mathbf{V}_0 := { \boldsymbol{v} \in \mathbf{V} : \operatorname{div} \boldsymbol{v} = 0 \text{ in } \Omega }.$$ (3.7)

Lemma 3.3 The bilinear form a(•, •) induces an invertible operator on V⁰ × L 2 (Ω).

Proof. We proceed using Lemma 3.1. First, from (3.4a) we observe that a(•, •) is bounded. To show that it also satisfies the inf-sup condition (3.1), we proceed as in, e.g., [3, Section 2.3]. For all (z, ζ) ∈ V⁰ × L 2 (Ω) (see (3.7)), we can define ^v^b := 2^z and ^θ^b := ^√ νcurl z − ζ, and then immediately assert that

$$a(\mathbf{z}, \boldsymbol{\zeta}; \widehat{\boldsymbol{v}}, \widehat{\boldsymbol{\theta}}) = 2\kappa^{-1} |\mathbf{z}|{0,\Omega}^2 + \nu |\mathbf{curl}\,\mathbf{z}|{0,\Omega}^2 + |\boldsymbol{\zeta}|{0,\Omega}^2 \ge \min{2\kappa^{-1}, \nu} |\mathbf{z}|{\mathbf{V}}^2 + |\boldsymbol{\zeta}|_{0,\Omega}^2.$$

Furthermore, it is clear that ^∥vb∥^V = 2∥z∥^V and ^∥θb∥0,^Ω ^≤ (1 + ^√ ν)(∥z∥^V + ∥ζ∥0,Ω), and from this, we can conclude that

$$\sup_{\mathbf{0} \neq (\boldsymbol{v}, \boldsymbol{\theta}) \in \mathbf{V}0 \times \mathbf{L}^2(\Omega)} \frac{a(\mathbf{z}, \boldsymbol{\zeta}; \boldsymbol{v}, \boldsymbol{\theta})}{|\boldsymbol{v}|{\mathbf{V}} + |\boldsymbol{\theta}|{0,\Omega}} \ge \frac{a(\mathbf{z}, \boldsymbol{\zeta}; \widehat{\boldsymbol{v}}, \widehat{\boldsymbol{\theta}})}{|\widehat{\boldsymbol{v}}|{\mathbf{V}} + |\widehat{\boldsymbol{\theta}}|{0,\Omega}} \ge \alpha \left( |\mathbf{z}|{\mathbf{V}} + |\boldsymbol{\zeta}|_{0,\Omega} \right) \quad \forall \left(\mathbf{z}, \boldsymbol{\zeta}\right) \in \mathbf{V}_0 \times \mathbf{L}^2(\Omega) \quad (3.8)$$

with $$\alpha := \frac{\min{2\kappa^{-1}, \nu, 1}}{2(3 + \sqrt{\nu})}$$ . Thus, the result follows.

On the other hand, from the equivalence H¹ 0 (Ω) = V (see [31, Lemma 2.5]), we have that b(•, •) satisfies the following inf-sup condition (see [31, Section 5.1])

$$\sup{\mathbf{0}\neq(\mathbf{v},\boldsymbol{\theta})\in\mathbf{V}\times\mathbf{L}^{2}(\Omega)}\frac{b(\mathbf{v},\boldsymbol{\theta};q)}{|\mathbf{v}|{\mathbf{V}}+|\boldsymbol{\theta}|{0,\Omega}}\geq\beta|q|{0,\Omega}\quad\forall\,q\in\mathrm{L}_{0}^{2}(\Omega).\tag{3.9}$$

Let us now define the bilinear form A : [V × L 2 (Ω)] × L 2 0 (Ω) × [V × L 2 (Ω)] × L 2 0 (Ω) → R as

$$\mathcal{A}(\mathbf{z}, \boldsymbol{\zeta}, r; \boldsymbol{v}, \boldsymbol{\theta}, q) := a(\mathbf{z}, \boldsymbol{\zeta}; \boldsymbol{v}, \boldsymbol{\theta}) + b(\mathbf{z}, \boldsymbol{\zeta}; q) + b(\boldsymbol{v}, \boldsymbol{\theta}; r). \tag{3.10}$$

Owing to (3.4a) and (3.4b), it is clear that A(•, •) is bounded. Moreover, from (3.8), (3.9) and [29, Proposition 2.36] it is not difficult to see that the following inf-sup condition holds:

$$\sup{\mathbf{0} \neq ((\boldsymbol{v}, \boldsymbol{\theta}), q) \in [\mathbf{V} \times \mathbf{L}^{2}(\Omega)] \times \mathbf{L}^{2}{0}(\Omega)} \frac{\mathcal{A}(\mathbf{z}, \boldsymbol{\zeta}, r; \boldsymbol{v}, \boldsymbol{\theta}, q)}{|((\boldsymbol{v}, \boldsymbol{\theta}), q)|} \ge \gamma |((\mathbf{z}, \boldsymbol{\zeta}), r)|$$ (3.11)

for all ((z, ζ), r) ∈ [V × L 2 (Ω)] × L 2 0 (Ω), where ∥((z, ζ), r)∥ := ∥z∥^V + ∥ζ∥0,^Ω + ∥r∥0,Ω, and

$$\gamma = \frac{\min{2\kappa^{-1}, \nu, 1}\beta^2}{(\beta + \kappa^{-1} + \sqrt{\nu} + 2)^2}.$$ (3.12)

3.2 Well-posedness analysis via Banach fixed-point

We proceed similarly to [27] using a fixed-point strategy to prove the well-posedness of (2.7). Let us introduce the bounded set

$$\mathbf{K} := \left{ \widehat{\boldsymbol{u}} \in \mathbf{V} : \quad |\widehat{\boldsymbol{u}}|{\mathbf{V}} \le \frac{2}{\gamma} |\boldsymbol{f}|{0,\Omega} \right}$$ (3.13)

with γ the constant defined in (3.12). Then, we define a fixed-point operator as

$$\mathscr{F}: \mathbf{K} \to \mathbf{K}, \quad \widehat{\mathbf{u}} \to \mathscr{F}(\widehat{\mathbf{u}}) = \mathbf{u},$$ (3.14)

where, given ^u^b ^∈ ^K, ^u is the first component of (u, ^ω), where ((u, ^ω), p) ^∈ [^V ^× ^L 2 (Ω)] × L 2 0 (Ω) is the solution of the linearised version of problem (2.7): Find ((u, ω), p) ∈ [V × L 2 (Ω)] × L 2 0 (Ω) such that

$$a(\boldsymbol{u}, \boldsymbol{\omega}; \boldsymbol{v}, \boldsymbol{\theta}) + b(\boldsymbol{v}, \boldsymbol{\theta}; p) + c^{\hat{\boldsymbol{u}}}(\boldsymbol{u}, \boldsymbol{\omega}; \boldsymbol{v}, \boldsymbol{\theta}) = F(\boldsymbol{v}, \boldsymbol{\theta}),$$

$b(\boldsymbol{u}, \boldsymbol{\omega}; q) = 0$

$$(3.15)$$

for all ((v, θ), q) ∈ [V × L 2 (Ω)] × L 2 0 (Ω). It is clear that ((u, ω), p) is a solution to (2.7) if and only if u satisfies F(u) = u, and consequently, the well-posedness of (2.7) is equivalent to the unique solvability of the fixed-point problem: Find u ∈ K such that

$$\mathscr{F}(\boldsymbol{u}) = \boldsymbol{u}.\tag{3.16}$$

In this way, in what follows we focus on proving the unique solvability of (3.16).

3.3 Well-definiteness of the fixed-point map

Let us first provide sufficient conditions under which the operator F (cf. (3.14)) is well-defined, or equivalently, the problem (3.15) is well-posed.

Lemma 3.4 (Unique solvability of the linearised problem) Let ^u^b ^∈ ^K and assume that

$$\frac{4}{\gamma^2}C{\rm S}^2 \left(\frac{1}{\sqrt{\nu}} + \mathsf{F}\right) |\boldsymbol{f}|{0,\Omega} \le 1 \tag{3.17}$$

with γ the positive constant in (3.12). Then, there exists a unique ((u, ω), p) ∈ [V × L 2 (Ω)] × L 2 0 (Ω) solution to (3.15). In addition, there holds

$$|((\boldsymbol{u},\boldsymbol{\omega}),p)| \le \frac{2}{\gamma} |\boldsymbol{f}|_{0,\Omega}.$$ (3.18)

Proof. We proceed similarly as in the proof of [19, Theorem 3.6]. In fact, given ^u^b ^∈ ^K, we begin by defining the bilinear form:

$$\mathcal{B}^{\widehat{u}}(\mathbf{z}, \boldsymbol{\zeta}, r; \boldsymbol{v}, \boldsymbol{\theta}, q) := \mathcal{A}(\mathbf{z}, \boldsymbol{\zeta}, r; \boldsymbol{v}, \boldsymbol{\theta}, q) + c^{\widehat{u}}(\mathbf{z}, \boldsymbol{\zeta}; \boldsymbol{v}, \theta)$$ (3.19)

with A(•, •) and c ub (•, •) the forms defined in (3.10) and (2.5c). Then, problem (3.15) can be rewritten equivalently as: Find ((u, ω), p) ∈ [V × L (Ω)] × L 2 0 (Ω), such that

$$\mathcal{B}^{\widehat{\boldsymbol{u}}}(\boldsymbol{u},\boldsymbol{\omega},p;\boldsymbol{v},\boldsymbol{\theta},q) = F(\boldsymbol{\theta},\boldsymbol{v}) \quad \forall ((\boldsymbol{v},\boldsymbol{\theta}),q) \in [\mathbf{V} \times \mathbf{L}^2(\Omega)] \times \mathbf{L}_0^2(\Omega). \tag{3.20}$$

Therefore, to prove the well-definiteness of F, in the sequel we equivalently prove that problem (3.20) is well-posed by means of Lemma 3.2. First, given ((z, ^ζ), r),((ve, ^θe), ^qe) ^∈ [^V ^× ^L 2 (Ω)] × L 2 0 (Ω) with ((ve, ^θe), ^qe) ̸⁼ ⁰, from (3.5) we observe that

$$\sup{\mathbf{0} \neq ((\boldsymbol{v},\boldsymbol{\theta}),q) \in [\mathbf{V} \times \mathbf{L}^{2}(\Omega)] \times \mathbf{L}{0}^{2}(\Omega)} \frac{\mathcal{B}^{\widehat{\boldsymbol{u}}}(\mathbf{z},\boldsymbol{\zeta},r;\boldsymbol{v},\boldsymbol{\theta},q)}{|((\boldsymbol{v},\boldsymbol{\theta}),q)|} \ge \frac{\left|\mathcal{A}(\mathbf{z},\boldsymbol{\zeta},r;\widetilde{\boldsymbol{v}},\widetilde{\boldsymbol{\theta}},\widetilde{q})\right|}{|((\widetilde{\boldsymbol{v}},\widetilde{\boldsymbol{\theta}}),\widetilde{q})|} - \frac{\left|c^{\widehat{\boldsymbol{u}}}(\mathbf{z},\boldsymbol{\zeta};\widetilde{\boldsymbol{v}},\widetilde{\boldsymbol{\theta}})\right|}{|((\widetilde{\boldsymbol{v}},\widetilde{\boldsymbol{\theta}}),\widetilde{q})|}$$ $$\ge \frac{\left|\mathcal{A}(\mathbf{z},\boldsymbol{\zeta},r;\widetilde{\boldsymbol{v}},\widetilde{\boldsymbol{\theta}},\widetilde{q})\right|}{|((\widetilde{\boldsymbol{v}},\widetilde{\boldsymbol{\theta}}),\widetilde{q})|} - C{\mathrm{S}}^{2}\left(\frac{1}{\sqrt{\nu}} + \mathsf{F}\right)|\widehat{\boldsymbol{u}}|{\mathbf{V}}(|\mathbf{z}|{\mathbf{V}} + |\boldsymbol{\zeta}|{0,\Omega}).$$

Together with the global inf-sup condition (3.11) and the fact that ((vb, ^θb), ^qe) is arbitrary, this implies

Hence, from the definition of the set K (cf. (3.13)), and assumption (3.17), we easily get

$$C{\mathrm{S}}^{2}\left(\frac{1}{\sqrt{\nu}} + \mathsf{F}\right) |\widehat{\boldsymbol{u}}|{\mathbf{V}} \leq \frac{2}{\gamma} C{\mathrm{S}}^{2}\left(\frac{1}{\sqrt{\nu}} + \mathsf{F}\right) |\boldsymbol{f}|{0,\Omega} \leq \frac{\gamma}{2} \tag{3.22}$$

and then, combining (3.21) and (3.22), we obtain

On the other hand, for a given ((z, ζ), r) ∈ [V × L 2 (Ω)] × L 2 0 (Ω), we observe that

$$\sup{\mathbf{0} \neq ((\boldsymbol{v}, \boldsymbol{\theta}), q) \in [\mathbf{V} \times \mathbf{L}^{2}(\Omega)] \times \mathbf{L}^{2}{0}(\Omega)} \mathcal{B}^{\widehat{\boldsymbol{u}}}(\boldsymbol{v}, \boldsymbol{\theta}, q; \mathbf{z}, \boldsymbol{\zeta}, r)$$

$$\geq \sup{\mathbf{0} \neq ((\boldsymbol{v}, \boldsymbol{\theta}), q) \in [\mathbf{V} \times \mathbf{L}^{2}(\Omega)] \times \mathbf{L}^{2}{0}(\Omega)} \frac{\mathcal{B}^{\widehat{\boldsymbol{u}}}(\boldsymbol{v}, \boldsymbol{\theta}, q; \mathbf{z}, \boldsymbol{\zeta}, r)}{|((\boldsymbol{v}, \boldsymbol{\theta}), q)|}$$

$$= \sup{\mathbf{0} \neq ((\boldsymbol{v}, \boldsymbol{\theta}), q) \in [\mathbf{V} \times \mathbf{L}^{2}(\Omega)] \times \mathbf{L}^{2}{0}(\Omega)} \frac{\mathcal{A}(\boldsymbol{v}, \boldsymbol{\theta}, q; \mathbf{z}, \boldsymbol{\zeta}, r) + c^{\widehat{\boldsymbol{u}}}(\boldsymbol{v}, \boldsymbol{\theta}; \mathbf{z}, \boldsymbol{\zeta})}{|((\boldsymbol{v}, \boldsymbol{\theta}), q)|}$$

with the problem definition in the last step. Putting this together with (3.5) implies

$$\sup{\mathbf{0} \neq ((\boldsymbol{v}, \boldsymbol{\theta}), q) \in [\mathbf{V} \times \mathbf{L}^{2}(\Omega)] \times \mathcal{L}{0}^{2}(\Omega)} \mathcal{B}^{\widehat{\boldsymbol{u}}}(\boldsymbol{v}, \boldsymbol{\theta}, q; \mathbf{z}, \boldsymbol{\zeta}, r)$$

$$\geq \sup{\mathbf{0} \neq ((\boldsymbol{v}, \boldsymbol{\theta}), q) \in [\mathbf{V} \times \mathbf{L}^{2}(\Omega)] \times \mathcal{L}{0}^{2}(\Omega)} \frac{\mathcal{A}(\boldsymbol{v}, \boldsymbol{\theta}, q; \mathbf{z}, \boldsymbol{\zeta}, r)}{|((\boldsymbol{v}, \boldsymbol{\theta}), q)|} - C{\mathbf{S}}^{2} \left(\frac{1}{\sqrt{\nu}} + \mathsf{F}\right) |\widehat{\boldsymbol{u}}|{\mathbf{V}} (|\mathbf{z}|{\mathbf{V}} + |\boldsymbol{\zeta}|{0,\Omega}). \tag{3.24}$$

Therefore, using the fact that B ub (•, •) is symmetric, from (3.11) and (3.24) we obtain

$$\sup_{\mathbf{0}\neq ((\boldsymbol{v},\boldsymbol{\theta}),q)\in [\mathbf{V}\times\mathbf{L}^2(\Omega)]\times\mathbf{L}^20(\Omega)}\mathcal{B}^{\widehat{\boldsymbol{u}}}(\boldsymbol{v},\boldsymbol{\theta},q;\mathbf{z},\boldsymbol{\zeta},r)\geq \left(\gamma-C{\mathrm{S}}^2\Big(\frac{1}{\sqrt{\nu}}+\mathsf{F}\Big)|\widehat{\boldsymbol{u}}|_{\mathbf{V}}\right)|((\mathbf{z},\boldsymbol{\zeta}),r)|.$$

Using also (3.22), yields

$$\sup{\mathbf{0} \neq ((\boldsymbol{v}, \boldsymbol{\theta}), q) \in [\mathbf{V} \times \mathbf{L}^{2}(\Omega)] \times \mathbf{L}{0}^{2}(\Omega)} \mathcal{B}^{\widehat{\boldsymbol{u}}}(\boldsymbol{v}, \boldsymbol{\theta}, q; \mathbf{z}, \boldsymbol{\zeta}, r) \ge \frac{\gamma}{2} |((\mathbf{z}, \boldsymbol{\zeta}), r)| > 0$$ (3.25)

for all ((z, ζ), r) ∈ [V × L 2 (Ω)] × L 2 0 (Ω).

In this way, from (3.23) and (3.25) we obtain that B ub (•, •) satisfies the hypotheses of Lemma 3.2, which allows us to conclude the unique solvability of (3.15), or equivalently, the existence of a unique ((u, ω), p) ∈ [V×L 2 (Ω)]×L 2 0 (Ω) such that F(u) = u. Finally, from (3.23), with ((z, ζ), r) = ((u, ω), p) and (3.20), we readily obtain that

$$|\boldsymbol{u}|{\mathbf{V}} \le |((\boldsymbol{u}, \boldsymbol{\omega}), p)| \le \frac{2}{\gamma} |\boldsymbol{f}|{0,\Omega}$$ (3.26)

implying that u belongs to K and concludes the proof. □

3.4 Well-posedness of the continuous problem

Now we provide the main result of this section, namely, the existence and uniqueness of solution of problem (2.7). This result is established in the following theorem.

Theorem 3.5 (Unique solvability) Let f ∈ L 2 (Ω) such that

$$\frac{4}{\gamma^2}C{\rm S}^2 \left(\frac{1}{\sqrt{\nu}} + \mathsf{F}\right) |f|{0,\Omega} < 1 \tag{3.27}$$

with γ the positive constant in (3.12). Then, F (cf. (3.14)) has a unique fixed-point u ∈ K. Equivalently, problem (2.7) has a unique solution ((u, ω), p) ∈ [K × L 2 (Ω)] × L 2 0 (Ω). Moreover, there holds

$$|((\boldsymbol{u},\boldsymbol{\omega}),p)| \le \frac{2}{\gamma} |\boldsymbol{f}|_{0,\Omega}. \tag{3.28}$$

Proof. Recall that (3.27) ensures the well-definiteness of ^F. Now, let ^ub1, ^ub2, ^u1, ^u² ^∈ ^K, be such that ^F(ub1) = ^u¹ and ^F(ub2) = ^u2. According to (3.14), it follows that there exist unique (ω1, p1), (ω2, p2) ∈ L 2 (Ω) × L 2 0 (Ω), such that for all ((v, θ), q) ∈ [V × L 2 (Ω)] × L 2 0 (Ω), there hold

$$\mathcal{B}^{\widehat{\boldsymbol{u}}_1}(\boldsymbol{u}_1,\boldsymbol{\omega}_1,p_1;\boldsymbol{v},\boldsymbol{\theta},q) = F(\boldsymbol{v},\boldsymbol{\theta})\,,\quad\text{and}\quad \mathcal{B}^{\widehat{\boldsymbol{u}}_2}(\boldsymbol{u}_2,\boldsymbol{\omega}_2,p_2;\boldsymbol{v},\boldsymbol{\theta},q) = F(\boldsymbol{v},\boldsymbol{\theta})\,.$$

Then, subtracting both equations, adding ±c ^ub¹ (u2, ω2; v, θ), and recalling the definition of B ub in (3.19), we easily arrive at

$$\mathcal{B}^{\widehat{\boldsymbol{u}}_1}(\boldsymbol{u}_1-\boldsymbol{u}_2,\boldsymbol{\omega}_1-\boldsymbol{\omega}_2,p_1-p_2;\boldsymbol{v},\boldsymbol{\theta},q) = c^{\widehat{\boldsymbol{u}}_2}(\boldsymbol{u}_2,\boldsymbol{\omega}_2;\boldsymbol{v},\boldsymbol{\theta}) - c^{\widehat{\boldsymbol{u}}_1}(\boldsymbol{u}_2,\boldsymbol{\omega}_2;\boldsymbol{v},\boldsymbol{\theta})\,.$$

Therefore, recalling that ^ub¹ ^∈ ^K from the latter identity, together with (3.23), the inequality ^|ub1| − ^|ub2| ≤ |ub¹ ⁻ ^ub2|, and simple computations, we obtain

$$\frac{\gamma}{2} | \boldsymbol{u}{1} - \boldsymbol{u}{2} |{\mathbf{V}} \leq \sup{\boldsymbol{0} \neq ((\boldsymbol{v},\boldsymbol{\theta}),q) \in [\mathbf{V} \times \mathbf{L}^{2}(\Omega)] \times \mathbf{L}{0}^{2}(\Omega)} \frac{\mathcal{B}^{\widehat{\boldsymbol{u}}{1}}(\boldsymbol{u}{1} - \boldsymbol{u}{2}, \boldsymbol{\omega}{1} - \boldsymbol{\omega}{2}, p{1} - p{2}; \boldsymbol{v}, \boldsymbol{\theta}, q)}{|((\boldsymbol{v},\boldsymbol{\theta}),q)|} \ = \sup{\boldsymbol{0} \neq ((\boldsymbol{v},\boldsymbol{\theta}),q) \in [\mathbf{V} \times \mathbf{L}^{2}(\Omega)] \times \mathbf{L}{0}^{2}(\Omega)} \frac{c^{\widehat{\boldsymbol{u}}{2}}(\boldsymbol{u}{2}, \boldsymbol{\omega}{2}; \boldsymbol{v}, \boldsymbol{\theta}) - c^{\widehat{\boldsymbol{u}}{1}}(\boldsymbol{u}{2}, \boldsymbol{\omega}{2}; \boldsymbol{v}, \boldsymbol{\theta})}{|((\boldsymbol{v},\boldsymbol{\theta}),q)|} \ \leq C{\mathrm{S}}^{2} \left(\frac{1}{\sqrt{\nu}} + \mathsf{F}\right) |\widehat{\boldsymbol{u}}{1} - \widehat{\boldsymbol{u}}{2}|{\mathbf{V}} (|\boldsymbol{u}{2}|{\mathbf{V}} + |\boldsymbol{\omega}{2}|{0,\Omega}).$$

Together with the fact that (u2, ω2) satisfy (3.18), this yields

$$|\mathscr{F}(\widehat{\boldsymbol{u}}_1) - \mathscr{F}(\widehat{\boldsymbol{u}}2)|{\mathbf{V}} = |\boldsymbol{u}_1 - \boldsymbol{u}2|{\mathbf{V}} \leq \frac{2}{\gamma} C{\mathrm{S}}^2 \Big(\frac{1}{\sqrt{\nu}} + \mathsf{F}\Big) \frac{2}{\gamma} |\boldsymbol{f}|{0,\Omega} |\widehat{\boldsymbol{u}}_1 - \widehat{\boldsymbol{u}}2|{\mathbf{V}}.$$

Combining the previous estimate with (3.27) and the Banach fixed-point theorem, readily implies that F has a unique fixed-point in K, and so there exists a unique ((u, ω), p) ∈ [V × L 2 (Ω)] × L 2 0 (Ω) solution to (2.7). Finally, the estimate (3.28) is obtained analogously to (3.26). □

Note that the formulation analysed above can also be defined in 2D. The vorticity is then the scalar ω = √ νcurlu, the operator curl is to be replaced by rot, and the weak convective term is now written as − √ ν R Ω ωu· v. The space for vorticity is then L² (Ω). At the discrete level these considerations also hold, but for sake of conciseness of the presentation we only discuss the 3D case.

4 Galerkin scheme

In this section we introduce the Galerkin scheme associated with problem (2.4), and show using Banach's fixed-point arguments that it admits a unique discrete solution.

4.1 Definition of the non-conforming method

First, let us denote by {Th} a family of non-degenerate simplicial meshes on Ω ⊂ R d (we simply assume that the domain is polytopal, so that no special treatment of the boundary is needed), and denote by E^h the set of all facets (edges in 2D) in the mesh, distinguishing between inner facets E int h and the set of facets lying on Γ, E Γ h . By h^K we denote the diameter of the element K and by h^F we denote the length/area of the facet F. As usual, by h we denote the maximum of the diameters of elements in Th. For a smooth vector, scalar, or tensor field ζ defined on Th, ζ ^± denote its traces taken from the interior of K⁺ and K−, respectively. We also denote by n ^± the outward unit normal vector to K^± (and for boundary faces it points outward of the domain Ω). For any inner facet F we define the normal and tangential jumps of any element-wise defined vector function v ∈ L 2 (Ω) across F by

$$[![\boldsymbol{v}\cdot\mathbf{n}]!]_F:=\boldsymbol{v}^+\cdot\mathbf{n}^++\boldsymbol{v}^-\cdot\mathbf{n}^-,\quad [![\boldsymbol{v}\times\mathbf{n}]!]_F:=\boldsymbol{v}^+\times\mathbf{n}^++\boldsymbol{v}^-\times\mathbf{n}^-$$

with K⁺ and K⁻ the two elements adjacent to F, and use the convention that [[v · n]]^F := v · n and [[v ×n]]^F := v ×n if F ∈ E^Γ h . For all meshes we assume that they are sufficiently regular (there exists a uniform positive constant η¹ such that each element K is star-shaped with respect to a ball of radius greater than η1hK). It is also assumed that there exists η² > 0 such that for each element and every facet F ∈ ∂K, we have that h^F ≥ η2hK, see, e.g., [29]). For ℓ ≥ 0, by Pℓ(K) we denote the space of polynomials of total degree at most ℓ defined locally on the generic element K ∈ Th.

For the approximation of velocity and pressure we use the nonconforming Crouzeix–Raviart Stokes inf-sup stable element (see [28]) where the velocity space consists of piecewise vector-valued d-linear polynomials on each dimension and continuous at the barycentre of the intra-element facets, the discrete pressure consist of piecewise constant functions, and for sake of inf-sup stability we also need that the curl of the discrete velocity lives in the space of vorticity and so we take piecewise vector-valued constants. This gives

$$\mathbf{V}{h} := { \boldsymbol{v}{h} \in \mathbf{L}^{2}(\Omega) : \boldsymbol{v}{h} \in \mathbb{P}{1}(K)^{d} \ \forall K \in \mathcal{T}{h}, \quad J{F}(\llbracket \boldsymbol{v}{h} \rrbracket{F}) = \mathbf{0} \ \forall F \in \mathcal{E}{h}^{\text{int}}, \quad J{F}(\boldsymbol{v}{h}|{F}) = \mathbf{0} \ \forall F \in \mathcal{E}_{h}^{\Gamma} },$$

$$\mathbf{W}{h} := { \boldsymbol{\theta}{h} \in \mathbf{L}^{2}(\Omega) : \ \boldsymbol{\theta}{h}|{K} \in \mathbb{P}{0}(K)^{d(d-1)/2} \ \ \forall K \in \mathcal{T}{h} },$$

$$\mathbf{Q}{h} := { q{h} \in \mathbf{L}{0}^{2}(\Omega) : \ q{h}|{K} \in \mathbb{P}{0}(K) \ \ \forall K \in \mathcal{T}_{h} }.$$

$$(4.1)$$

For any facet F ∈ E^h with barycentre C^F , the nodal functional J^F is defined by

$$J_F(\mathbf{v}) = \mathbf{v}(C_F)$$ or $J_F(\mathbf{v}) = \frac{1}{h_F} \int_F \mathbf{v} \, \mathrm{d}s$

and the degrees of freedom associated with facets on E Γ h vanish for any v^h ∈ Vh. We recall the definition of the Crouzeix–Raviart interpolation I CR : V → V^h as

$$\mathcal{I}^{\mathsf{CR}} v(C_F) = J_F(v) \qquad \forall F \in \mathcal{E}_h.$$

We also recall that element-wise integration by parts on a given K ∈ T^h readily gives that I CR preserves the averages of first derivatives.

Next, defining the lowest-order Raviart–Thomas space

$$\mathbf{RT}_0(\mathcal{T}_h) := { \boldsymbol{v}_h \in \mathbf{H}(\operatorname{div},\Omega) : \forall K \in \mathcal{T}_h, \ \exists \ \mathbf{c}_K \in \mathbb{R}^d, a_K \in \mathbb{R} : \ \boldsymbol{v}_h|_K(\boldsymbol{x}) = \mathbf{c}_K + a_K \boldsymbol{x} }$$

we recall the Raviart–Thomas interpolation I RT : V ⊕ V^h → RT0(Th) as

$$\mathbf{n}_F \cdot \mathcal{I}^{\mathsf{RT}} \boldsymbol{v} = \frac{1}{h_F} \int_F \boldsymbol{v} \cdot \mathbf{n}_F \qquad \forall F \in \mathcal{E}_h. \tag{4.2}$$

Note that even if RT0(Th) ̸⊂ V^h (since the tangential components of Raviart–Thomas functions are not necessarily continuous at each C^F ), the interpolation is well-defined for v ∈ V^h and we have that (see, e.g., [43])

I RTI CRv = I RTv ∀v ∈ V.

For the subsequent analysis we consider the following broken norm for the space V^h

$$|\boldsymbol{v}{h}|{h}^{2} := \sum{K \in \mathcal{T}{h}} \left(\frac{1}{\kappa} |\boldsymbol{v}{h}|{0,K}^{2} + \nu |\mathbf{curl}\,\boldsymbol{v}{h}|{0,K}^{2} + |\operatorname{div}\,\boldsymbol{v}{h}|{0,K}^{2}\right)$$

$$+ \sum{F \in \mathcal{E}{int}^{int}} \frac{1}{h{F}} \left(\nu |[![\boldsymbol{v}{h} \times \mathbf{n}]!]{F}|{0,F}^{2} + |[![\boldsymbol{v}{h} \cdot \mathbf{n}]!]{F}|_{0,F}^{2}\right)$$

$$(4.3)$$

as well as the piecewise H¹ (Th)-seminorm

$$|!|!| \boldsymbol{v} |!|!|^2{\text{pw}} := \sum{K \in \mathcal{T}h} |\boldsymbol{\nabla} \boldsymbol{v}|^2{1,K}.$$

We also have that the Raviart–Thomas interpolator is stable on V and also on V^h

$$|\mathcal{I}^{\mathsf{RT}}\boldsymbol{v}|{h} \le C{\mathsf{RT}}|\boldsymbol{v}|{h} \quad \forall \, \boldsymbol{v} \in \mathbf{V} \cup \mathbf{V}{h}. \tag{4.4}$$

The following approximability bounds are known for Crouzeix–Raviart and Raviart–Thomas interpolants

$$|\boldsymbol{v} - \boldsymbol{\mathcal{I}}^{\mathsf{CR}} \boldsymbol{v}|{h} \le C{\mathsf{CR}} h |\boldsymbol{v}|{2,\Omega} \qquad \forall \, \boldsymbol{v} \in \mathbf{H}{0}^{2}(\Omega),$$ (4.5a)

$$|\boldsymbol{v} - \mathcal{I}^{\mathsf{RT}}\boldsymbol{v}|_{0,\Omega} \le C_F h |\boldsymbol{v}|_h \qquad \forall \boldsymbol{v} \in \mathbf{V} \cup \mathbf{V}_h$$ (4.5b)

with (4.5a) stated in [15, Section 2.3], and the constant CCR depends only on the mesh regularity. In addition, the constant C^F only depends on the shape of the triangles/tetrahedra (maximum angle) but not on their size (see, e.g., [23, 35]). Let P^h denote the L² projection operator, which satisfies the following approximation property (see [30, Theorem 3.6]):

$$|\mathcal{P}h q - q|{0,\Omega} \le CP h|q|{1,\Omega} \qquad \forall q \in H^1(\Omega).$$ (4.6)

Since the method is nonconforming in the velocity space, for the discrete setting we will require the broken curl and broken divergence operators (associated with the non-diagonal part of the bilinear form ah(•, •) and the discrete bilinear form bh(•, •))

$$\operatorname{\mathbf{curl}}_h: \mathbf{V} \oplus \mathbf{V}_h \to \mathbf{L}^2(\Omega), \quad \operatorname{div}_h: \mathbf{V} \oplus \mathbf{V}_h \to \mathbf{L}^2(\Omega)$$

in the following sense

$$(\operatorname{\mathbf{curl}}_h v_h)|_K := \operatorname{\mathbf{curl}}(v_h|_K) \quad \text{and} \quad (\operatorname{\mathrm{div}}_h v_h)|_K := \operatorname{\mathrm{div}}(v_h|_K) \qquad \forall K \in \mathcal{T}_h.$$

With these ingredients, we can define element-wise variational forms. The forms that require modification are as follows

$$a{h}(\boldsymbol{u}{h}, \boldsymbol{\omega}{h}; \boldsymbol{v}{h}, \boldsymbol{\theta}{h}) := \frac{1}{\kappa} \int{\Omega} \boldsymbol{u}{h} \cdot \mathcal{I}^{\mathsf{RT}} \boldsymbol{v}{h} + \sum{F \in \mathcal{E}{h}^{\mathsf{int}}} \frac{\vartheta}{h{F}} \int{F} (\nu \llbracket \boldsymbol{u}{h} \times \mathbf{n} \rrbracket{F} \cdot \llbracket \boldsymbol{v}{h} \times \mathbf{n} \rrbracket{F} + \llbracket \boldsymbol{u}{h} \cdot \mathbf{n} \rrbracket{F} \llbracket \boldsymbol{v}{h} \cdot \mathbf{n} \rrbracket{F}$$

$$+ \sqrt{\nu} \int{\Omega} \boldsymbol{\omega}{h} \cdot \mathbf{curl}{h} \boldsymbol{v}{h} + \sqrt{\nu} \int{\Omega} \boldsymbol{\theta}{h} \cdot \mathbf{curl}{h} \boldsymbol{u}{h} - \int{\Omega} \boldsymbol{\omega}{h} \cdot \boldsymbol{\theta}_{h},$$

$$= \frac{1}{\kappa} \sum{K \in \mathcal{T}{h}} \int{K} \boldsymbol{u}{h} \cdot \mathcal{I}^{\mathsf{RT}} \boldsymbol{v}_{h}$$

$$+ \sum{F \in \mathcal{E}{h}^{\mathsf{int}}} \frac{\vartheta}{h{F}} \int{F} (\nu \llbracket \boldsymbol{u}{h} \times \mathbf{n} \rrbracket{F} \cdot \llbracket \boldsymbol{v}{h} \times \mathbf{n} \rrbracket{F} + \llbracket \boldsymbol{u}{h} \cdot \mathbf{n} \rrbracket{F} \llbracket \boldsymbol{v}{h} \cdot \mathbf{n} \rrbracket{F})$$

$$+ \sqrt{\nu} \sum{K \in \mathcal{T}{h}} \int{K} \boldsymbol{\omega}{h} \cdot \mathbf{curl} \boldsymbol{v}{h} + \sqrt{\nu} \sum{K \in \mathcal{T}{h}} \int{K} \boldsymbol{\theta}{h} \cdot \mathbf{curl} \boldsymbol{u}{h} - \int{\Omega} \boldsymbol{\omega}{h} \cdot \boldsymbol{\theta}_{h}, \qquad (4.7a)$$

$$b_h(\boldsymbol{v}_h, \boldsymbol{\theta}_h; qh) := -\int{\Omega} q_h \operatorname{div}_h \boldsymbol{v}h = -\sum{K \in \mathcal{T}_h} \int_K q_h \operatorname{div} \boldsymbol{v}_h, \tag{4.7b}$$

$$c_h^{\widehat{\boldsymbol{u}}_h}(\boldsymbol{u}_h, \boldsymbol{\omega}_h; \boldsymbol{v}_h, \boldsymbol{\theta}h) := -\frac{1}{\sqrt{\nu}} \int{\Omega} (\widehat{\boldsymbol{u}}_h \times \boldsymbol{\omega}_h) \cdot \mathcal{I}^{\mathsf{RT}} \boldsymbol{v}h + \mathsf{F} \int{\Omega} |\widehat{\boldsymbol{u}}_h| \, \boldsymbol{u}_h \cdot \mathcal{I}^{\mathsf{RT}} \boldsymbol{v}_h, \tag{4.7c}$$

$$F_h(\boldsymbol{v}_h, \boldsymbol{\theta}h) := \int{\Omega} \boldsymbol{f} \cdot \mathcal{I}^{\mathsf{RT}} \boldsymbol{v}_h \tag{4.7d}$$

with ϑ > 0 a sufficiently large, user specified penalty parameter. The stabilisation in ah(•, •) uses normal and tangential jumps across inter-element boundaries, which are needed for controlling the consistency error and in general for the convergence of the scheme, as discussed in, e.g., [33, 36] for elasticity equations (see also for example [37, 48] for the case of nonconforming schemes on quadrilaterals). We also provide numerical evidence in Section 7 that if ϑ = 0 then the method does not converge. Note also that, in [16] the jumps do not require a penalisation parameter since in that formulation the curl-curl and div-div terms are explicitly present in the continuous and discrete bilinear form (and the jump terms only contribute to maintain consistency). Also, note that the 2D Crouzeix–Raviart space used in [16] is also element-wise divergence-free, but the underlying continuous space only sets tangential components on boundary edges. Another variant in [34] imposes continuity only of the tangential components at the edges' midpoints.

The interpolation of the test velocity in the right-hand side functional (4.7d) follows the definition proposed in [40], but we stress that one could use any smoother operator such that the velocity error (in the broken norm (4.3)) is proportional to the corresponding best approximation error [51]. We proceed similarly for the convective and Forchheimer nonlinearities in (4.7c), as well as for the Brinkman term.

Having introduced the additional notations described above, the nonlinear discrete problem consists in finding ((uh, ωh), ph) ∈ (V^h × Wh) × Qh, such that:

$$a_h(\boldsymbol{u}_h, \boldsymbol{\omega}_h; \boldsymbol{v}_h, \boldsymbol{\theta}_h) + b_h(\boldsymbol{v}_h, \boldsymbol{\theta}_h; p_h) + c_h^{\boldsymbol{u}_h}(\boldsymbol{u}_h, \boldsymbol{\omega}_h; \boldsymbol{v}_h, \boldsymbol{\theta}_h) = F_h(\boldsymbol{v}_h, \boldsymbol{\theta}_h),$$

$$b_h(\boldsymbol{u}_h, \boldsymbol{\omega}_h; q_h) = 0$$

$$(4.8)$$

for all ((vh, θh), qh) ∈ (V^h × Wh) × Qh. Note that the interpolated test discrete velocity on the right-hand side functional induce a variational crime approach that maps discretely divergence-free test functions to divergence-free functions in H(div, Ω). This can also be regarded as a smoothing approach that permits to have a discrete load F^h well defined for all continuous functionals on V. This setting has been used extensively in, e.g., [1, 51, 40, 43], with the additional aim of achieving pressure robustness of the formulation. We also recall that interpolated test velocities are used in the convective nonlinearity. Finally, using the Cauchy–Schwarz and H¨older inequalities, it is clear that the bilinear forms a^h and b^h are bounded

$$a_h(\boldsymbol{u}_h, \boldsymbol{\omega}_h; \boldsymbol{v}_h, \boldsymbol{\theta}_h) \leq (|\boldsymbol{u}_h|_h + |\boldsymbol{\omega}h|{0,\Omega})(|\boldsymbol{v}_h|_h + |\boldsymbol{\theta}h|{0,\Omega}),$$

$$b_h(v_h, \theta_h; q_h) \le ||qh||{0,\Omega}(||v_h||_h + ||\thetah||{0,\Omega})$$

as well as c ubh h and Fh:

$$|\widehat{c}{h}^{\widehat{\boldsymbol{u}}{h}}(\boldsymbol{u}{h},\boldsymbol{\omega}{h};\boldsymbol{v}{h},\boldsymbol{\theta}{h})| \leq C{c{h}} |\widehat{\boldsymbol{u}}{h}|{h} (|\boldsymbol{u}{h}|{h} + |\boldsymbol{\omega}{h}|{0,\Omega}) (|\boldsymbol{v}{h}|{h} + |\boldsymbol{\theta}{h}|{0,\Omega}), \tag{4.9a}$$

$$\left| F_h(\boldsymbol{v}_h, \boldsymbol{\theta}h) \right| \le C{Fh} |\boldsymbol{f}|{0,\Omega} (|\boldsymbol{v}_h|_h + |\boldsymbol{\theta}h|{0,\Omega}) \tag{4.9b}$$

with Cc^h > 0 and CF^h > 0 depending on the boundedness constant of the operator I RT denoted by CRT, as well as on the penalty parameter ϑ.

4.2 Further properties of the discrete problem

Discrete kernel properties. First, we denote the kernel of the bilinear form bh(•, •) as V⁰ ^h × W^h (noting that the discrete vorticity space does not play an active role), and from [35, Lemma 4.62] we can see that, since the broken divergence of an element-wise affine function is element-wise constant, we can readily choose as test function q^h = divhvh, yielding the characterisation

$$\mathbf{V}_h^0 := \left{ \mathbf{v}_h \in \mathbf{V}_h : \operatorname{div}_h \mathbf{v}_h = 0 \right}. \tag{4.10}$$

Similarly, we stress that

$$\operatorname{curl}{h} \boldsymbol{v}{h} \in \mathbf{W}{h} \qquad \forall \boldsymbol{v}{h} \in \mathbf{V}_{h}.$$ (4.11)

The following lemma corresponds to the discrete version of Lemma 3.3. It depends on a mesh size smallness assumption, which can easily be avoided – and the proof thus further simplified – either if we have a discrete K¨orn-type inequality using the broken curl part of the discrete velocity norm (which is indeed valid trivially in the 2D case thanks to a discrete Poincar´e inequality for Crouzeix–Raviart elements [16]), or if the first term in the definition of ah(•, •) is symmetric (for example, if it has also the Raviart–Thomas interpolation applied to the trial discrete function). We opt to keep the present form as it makes the a posteriori analysis more straightforward.

Lemma 4.1 (Invertibility on the kernel) The restriction of ah(•, •) to the kernel of bh(•, •) induces an invertible operator, provided that the mesh size h is sufficiently small:

$$h \le \frac{1}{2C_F\sqrt{\nu}}.\tag{4.12}$$

Proof. First, we note that for any z^h ∈ Vh, from the definition (4.3), it readily holds

$$|\mathbf{z}h|{0,\Omega} \le \sqrt{\kappa} |\mathbf{z}_h|_h. \tag{4.13}$$

Then, we can assert that

$$(\mathbf{z}{h}, \mathcal{I}^{\mathsf{RT}}\mathbf{z}{h}){0,\Omega} = \frac{1}{\kappa} (\mathbf{z}{h}, \mathcal{I}^{\mathsf{RT}}\mathbf{z}{h} - \mathbf{z}{h}){0,\Omega} + (\mathbf{z}{h}, \mathbf{z}{h}){0,\Omega}$$

$$\geq -|\mathbf{z}{h} - \mathcal{I}^{\mathsf{RT}}\mathbf{z}{h}|{0,\Omega} |\mathbf{z}{h}|{0,\Omega} + |\mathbf{z}{h}|_{0,\Omega}^{2}$$

$$\geq -C{F}h\sqrt{\kappa}|\mathbf{z}{h}|{h}^{2} + |\mathbf{z}{h}|_{0,\Omega}^{2}$$

$$(4.14)$$

having used the Cauchy–Schwarz inequality, as well as (4.5b) and (4.13) in the last step.

Next, and similarly to the continuous case, it is clear that ah(•, •) is bounded. In addition, for all (zh, ζ^h ) ∈ V⁰ ^h×W^h (see (4.10)), and owing to (4.11), we can define ^vb^h := 2z^h and ^θ^b ^h := ^√ νcurl zh−ζ^h , and invoke (4.14), from which we obtain

$$\begin{split} &a_h(\mathbf{z}_h, \boldsymbol{\zeta}_h; \widehat{\boldsymbol{v}}_h, \widehat{\boldsymbol{\theta}}_h) \ &= \frac{2}{\kappa} (\mathbf{z}_h, \mathcal{I}^{\mathsf{RT}} \mathbf{z}h){0,\Omega} + \sum_{F \in \mathcal{E}_h^{\mathrm{int}}} \frac{2\vartheta}{h_F} \int_F \left( \nu \llbracket \boldsymbol{u}_h \times \mathbf{n} \rrbracket_F^2 + \llbracket \boldsymbol{u}_h \cdot \mathbf{n} \rrbracketF^2 \right) + \nu \sum{K \in \mathcal{T}_h} |\mathbf{curl} \, \mathbf{z}h|{0,K}^2 + |\boldsymbol{\zeta}h|{0,\Omega}^2 \end{split}$$

$$\begin{split} & \geq \frac{2}{\kappa} |\mathbf{z}h|{0,\Omega}^2 - \frac{2C_F h}{\sqrt{\kappa}} |\mathbf{z}_h|h^2 + \sum{F \in \mathcal{E}_h^{\mathrm{int}}} \frac{2\vartheta}{h_F} \int_F \left(\nu \llbracket \boldsymbol{u}_h \times \mathbf{n} \rrbracket_F^2 + \llbracket \boldsymbol{u}_h \cdot \mathbf{n} \rrbracketF^2 \right) \ & + \nu \sum{K \in \mathcal{T}_h} |\mathbf{curl} \, \mathbf{z}h|{0,K}^2 + |\boldsymbol{\zeta}h|{0,\Omega}^2 \ & \geq 2 \bigg( 1 - \frac{C_F h}{\sqrt{\kappa}} \bigg) |\mathbf{z}_h|_h^2 + |\boldsymbol{\zeta}h|{0,\Omega}^2. \end{split}$$

Thus, using (4.12) together with the fact that ^∥vbh∥^h = 2∥zh∥^h and ^∥θ^b ^h∥0,^Ω ^≤ (1 + ^√ ν)(∥zh∥^h + ∥ζh∥0,Ω), we can conclude that

$$\sup_{\mathbf{0}\neq(\boldsymbol{v}_h,\boldsymbol{\theta}_h)\in\mathbf{V}_h^0\times\mathbf{Q}}\frac{a_h(\mathbf{z}_h,\boldsymbol{\zeta}_h;\boldsymbol{v}_h,\boldsymbol{\theta}_h)}{|\boldsymbol{v}_h|_h+|\boldsymbol{\theta}h|{0,\Omega}}\geq \frac{a_h(\mathbf{z}_h,\boldsymbol{\zeta}_h;\widehat{\boldsymbol{v}}_h,\widehat{\boldsymbol{\theta}}_h)}{|\widehat{\boldsymbol{v}}_h|_h+|\widehat{\boldsymbol{\theta}}h|{0,\Omega}}\geq \alpha_h\left(|\mathbf{z}_h|_h+|\boldsymbol{\zeta}h|{0,\Omega}\right)$$ (4.15)

for all (zh, ζ^h ) ∈ V⁰ ^h ^× Q, with ^α^h := ¹ 2(3+^√ ν) . Thus, the result follows directly from Lemma 3.1. □

Discrete inf-sup conditions. The discrete inf-sup condition for the discrete divergence operator is satisfied for Crouzeix–Raviart elements [28]. This is recalled in the following result.

Lemma 4.2 The pair (Vh, Qh) is inf-sup stable with constant β^h > 0 independent of the mesh size

$$0 < \betah := \inf{q_h \in \mathbf{Q}h \setminus {0}} \sup{\boldsymbol{v}_h \in \mathbf{V}_h \setminus {\mathbf{0}}} \frac{b_h(\boldsymbol{v}_h, \boldsymbol{\theta}_h; q_h)}{|\boldsymbol{v}_h|_h |qh|{0,\Omega}}.$$ (4.16)

In addition, we have the following properties, shown in [43, Lemma 3.1].

Lemma 4.3 For all v^h ∈ V^h there holds

$$b_h(\boldsymbol{v}_h, \bullet; q) = b\left(\mathcal{I}^{\mathsf{RT}}\boldsymbol{v}_h, \bullet; q\right) \qquad \forall q \in L^2(\Omega),$$ (4.17a)

$$b\left(\mathcal{I}^{\mathsf{RT}}\boldsymbol{v}{h},\bullet;q\right) = \int{\Omega} \nabla q \cdot \left(\mathcal{I}^{\mathsf{RT}}\boldsymbol{v}_{h}\right) \qquad \forall q \in \mathcal{H}^{1}(\Omega). \tag{4.17b}$$

4.3 Discrete fixed-point arguments

Here we proceed similarly as in the continuous setting. To begin with, we define the bilinear form A^h : [V^h × Wh] × Q^h × [V^h × Wh] × Qh) → R as

$$\mathcal{A}_h(\mathbf{z}_h, \boldsymbol{\zeta}_h, r_h; \boldsymbol{v}_h, \boldsymbol{\theta}_h, q_h) := a_h(\boldsymbol{u}_h, \boldsymbol{\omega}_h; \boldsymbol{v}_h, \boldsymbol{\theta}_h) + b_h(\boldsymbol{v}_h, \boldsymbol{\theta}_h; p_h) + b_h(\boldsymbol{u}_h, \boldsymbol{\omega}_h; q_h). \tag{4.18}$$

It is easy to see that A^h is bounded (since a^h and b^h are), and furthermore, using (4.15), (4.16), and [29, Proposition 2.36], we have that A^h satisfies the following inf-sup condition

$$\sup_{\mathbf{0} \neq ((\boldsymbol{v}_h, \boldsymbol{\theta}_h), q_h) \in [\mathbf{V}_h \times \mathbf{W}_h] \times \mathbf{Q}_h} \frac{\mathcal{A}_h(\mathbf{z}_h, \boldsymbol{\zeta}_h, r_h; \boldsymbol{v}_h, \boldsymbol{\theta}_h, q_h)}{|((\boldsymbol{v}_h, \boldsymbol{\theta}_h), q_h)|_h} \ge \widetilde{\gamma} |((\mathbf{z}_h, \boldsymbol{\zeta}_h), r_h)|_h$$ (4.19)

for all ((zh, ζ^h ), rh) ∈ [V^h ×Wh]×Qh, where ∥((zh, ζ^h ), rh)∥^h := ^∥zh∥^h ⁺∥ζh∥0,^Ω ⁺∥rh∥0,Ω, and γ >^e ⁰ is the discrete version of γ (cf. (3.12)).

Let us introduce the following set

$$\mathbf{K}{h} := \left{ \widehat{\boldsymbol{u}} \in \mathbf{V}{h} : |\widehat{\boldsymbol{u}}{h}|{h} \leq \frac{2}{\widetilde{\gamma}} C{F{h}} |\boldsymbol{f}|_{0,\Omega} \right}$$

$$(4.20)$$

with ^γ^e the global inf-sup constant defined in (4.19). Then, and again analogously to the continuous case, we define the following fixed-point operator

$$\mathscr{F}_h: \mathbf{K}_h \to \mathbf{K}_h, \quad \widehat{\boldsymbol{u}}_h \to \mathscr{F}_h(\widehat{\boldsymbol{u}}_h) = \boldsymbol{u}_h,$$ (4.21)

where, given ^ub^h ^∈ ^Kh, ^u^h is the first component of ((uh, ^ωh), ph) ^∈ [V^h ^× ^Wh] ^× ^Qh, the solution of the linearised version of problem (4.8): Find ((uh, ωh), ph) ∈ [V^h × Wh] × Q^h such that

$$a{h}(\boldsymbol{u}{h},\boldsymbol{\omega}{h};\boldsymbol{v}{h},\boldsymbol{\theta}{h}) + b{h}(\boldsymbol{v}{h},\boldsymbol{\theta}{h};p{h}) + c{h}^{\hat{\boldsymbol{u}}{h}}(\boldsymbol{u}{h},\boldsymbol{\omega}{h};\boldsymbol{v}{h},\boldsymbol{\theta}{h}) = F{h}(\mathcal{I}^{\mathsf{RT}}\boldsymbol{v}{h},\boldsymbol{\theta}{h}),$$

$$b{h}(\boldsymbol{u}{h},\boldsymbol{\omega}{h};q{h}) = 0,$$

$$(4.22)$$

for all ((vh, θh), qh) ∈ [V^h × Wh] × Qh.

It is clear that ((uh, ωh), ph) is a solution to (4.8) if and only if u^h satisfies Fh(uh) = uh, and consequently, the well-posedness of (4.8) is equivalent to the unique solvability of the fixed-point problem: Find u^h ∈ K^h such that

$$\mathscr{F}_h(\boldsymbol{u}_h) = \boldsymbol{u}_h. \tag{4.23}$$

In what follows we focus on (4.23). We start by establishing that F^h is well-defined.

Lemma 4.4 (Wellposedness of the discrete linearised problem) Let ^ub^h ^∈ ^K^h and assume that

$$\frac{4}{\tilde{\gamma}^2} C_{ch} C{Fh} | \boldsymbol{f} |{0,\Omega} \le 1 \tag{4.24}$$

with the positive constant ^γ^e in (4.19). Then, there exists a unique ((uh, ^ωh), ph) ^∈ [V^h ^× ^Wh] ^× ^Q^h solution to (4.22). In addition, there holds

$$|((\boldsymbol{u}_h, \boldsymbol{\omega}_h), p_h)|h \le \frac{2}{\tilde{\gamma}} C{Fh} |\boldsymbol{f}|{0,\Omega}. \tag{4.25}$$

Proof. Given ^ub^h ^∈ ^Kh, we proceed as in the proof of Lemma 3.4 and define the bilinear form

$$\mathcal{B}{h}^{\widehat{\boldsymbol{u}}{h}}(\mathbf{z}{h},\boldsymbol{\zeta}{h},r{h};\boldsymbol{v}{h},\boldsymbol{\theta}{h},q{h}):=\mathcal{A}{h}(\mathbf{z}{h},\boldsymbol{\zeta}{h},r{h};\boldsymbol{v}{h},\boldsymbol{\theta}{h},q{h})\,+\,c{h}^{\widehat{\boldsymbol{u}}{h}}(\mathbf{z}{h},\boldsymbol{\zeta}{h};\boldsymbol{v}{h},\boldsymbol{\theta}_{h}). \tag{4.26}$$

Using (4.9a), (4.19), (4.24) and [29, Proposition 2.36] we obtain the inf-sup condition

$$\sup_{\mathbf{0} \neq ((\mathbf{v}_h, \boldsymbol{\theta}_h), q_h) \in [\mathbf{V}_h \times \mathbf{W}_h] \times \mathbf{Q}_h} \frac{\mathcal{B}_h^{\widehat{\mathbf{u}}_h}(\mathbf{z}_h, \boldsymbol{\zeta}_h, r_h; \mathbf{v}_h, \boldsymbol{\theta}_h, q_h)}{|((\mathbf{v}_h, \boldsymbol{\theta}_h), q_h)|_h} \ge \frac{\widetilde{\gamma}}{2} |((\boldsymbol{\zeta}_h, \mathbf{z}_h), r_h)|_h$$ (4.27)

for all ((vh, θh), qh) ∈ [V^h ×Wh] × Qh. Therefore, owing to the fact that for finite dimensional linear problems surjectivity and injectivity are equivalent, from (4.27) and Lemma 3.2 we obtain that there exists a unique ((uh, ωh), ph) ∈ [V^h × Wh] × Q^h satisfying (4.22) with u^h ∈ Kh. □

The following theorem establishes the well-posedness of the nonlinear discrete problem (4.8).

Theorem 4.5 (Unique solvability of the discrete nonlinear problem) Let f ∈ L 2 (Ω) such that

$$\frac{4}{\tilde{\gamma}^2} C_{ch} C{Fh} | \mathbf{f} |{0,\Omega} < 1 \tag{4.28}$$

with ^γ^e the positive constant in (4.19). Then, ^F^h (cf. (4.21)) has a unique fixed-point ^u^h ^∈ ^Kh. Equivalently, problem (4.8) has a unique solution ((uh, ωh), ph) ∈ [K^h × Wh] × Qh. This discrete solution satisfies

$$|((\boldsymbol{u}_h, \boldsymbol{\omega}_h), p_h)|h \le \frac{2}{\tilde{\gamma}} C{Fh} |\boldsymbol{f}|{0,\Omega}. \tag{4.29}$$

Proof. Employing (4.27) and (4.25), along with (4.28), the proof follows adapting the steps developed in the proof of Theorem 3.5. Further details are omitted. □

5 A priori error bounds

Now we turn to the error analysis. First we derive a Strang-type estimate. Then, under a small data assumption we show linear convergence of the method in the energy norm. Finally, we show that the velocity-vorticity error is independent of the pressure error.

Lemma 5.1 (C´ea estimate) Let ((u, ω), p) ∈ [V×L 2 (Ω)]×L 2 0 (Ω) and ((uh, ωh), ph) ∈ [Vh×Wh]× Qh, the solution of (2.7) and (4.8), respectively. Then there hold the error estimate

$$|((\boldsymbol{u},\boldsymbol{\omega}),p) - ((\boldsymbol{u}{h},\boldsymbol{\omega}{h}),p{h})|{h} \leq \left(1 + \frac{2}{\widetilde{\gamma}}\right) \inf{((\mathbf{z}{h},\boldsymbol{\zeta}{h}),r{h})\in[\mathbf{V}{h}\times\mathbf{W}{h}]\times\mathbf{Q}{h}} |((\boldsymbol{u},\boldsymbol{\omega}),p) - ((\mathbf{z}{h},\boldsymbol{\zeta}{h}),r{h})|_{h}$$

$$+ \frac{2}{\widetilde{\gamma}} \sup{\mathbf{0}\neq((\mathbf{v}{h},\boldsymbol{\theta}{h}),q{h})\in[\mathbf{V}{h}\times\mathbf{W}{h}]\times\mathbf{Q}{h}} \frac{\mathcal{B}{h}^{\boldsymbol{u}{h}}(\boldsymbol{u},\boldsymbol{\omega},p;\boldsymbol{v}{h},\boldsymbol{\theta}{h},q{h}) - F{h}(\boldsymbol{v}{h},\boldsymbol{\theta}{h})}{|((\boldsymbol{v}{h},\boldsymbol{\theta}{h}),q{h})|_{h}}.$$ (5.1)

Proof. Let ((zh, ζ^h ), rh) ∈ [V^h × Wh] × Q^h be arbitrary, we have the decomposition

$$\mathbf{u} - \mathbf{u}_h = \mathbf{u} - \mathbf{z}_h + (\mathbf{z}_h - \mathbf{u}_h) = \mathbf{u} - \mathbf{z}h + \chi{\mathbf{u}}, \quad \boldsymbol{\omega} - \boldsymbol{\omega}_h = \boldsymbol{\omega} - \zeta_h + (\zeta_h - \boldsymbol{\omega}_h) = \boldsymbol{\omega} - \boldsymbol{\omega}h + \chi{\boldsymbol{\omega}},$$

$$p - p_h = p - r_h + (r_h - p_h) = p - r_h + \chi_p.$$ (5.2)

Then, from (4.27), using (4.8) and the Cauchy–Schwarz inequality, we have

$$\frac{\widetilde{\gamma}}{2} | ((\boldsymbol{\chi}{\boldsymbol{u}}, \boldsymbol{\chi}{\boldsymbol{\omega}}), \chi{p}) |{h} \leq \sup{\boldsymbol{0} \neq ((\boldsymbol{v}{h}, \boldsymbol{\theta}{h}), q{h}) \in [\mathbf{V}{h} \times \mathbf{W}{h}] \times \mathbf{Q}{h}} \frac{\mathcal{B}{h}^{\boldsymbol{u}{h}}(\boldsymbol{\chi}{\boldsymbol{u}}, \boldsymbol{\chi}{\boldsymbol{\omega}}, \chi{p}; \boldsymbol{v}{h}, \boldsymbol{\theta}{h}, q{h})}{| ((\boldsymbol{v}{h}, \boldsymbol{\theta}{h}), q{h}) |{h}} \ \leq | ((\boldsymbol{u}, \boldsymbol{\omega}), p) - ((\mathbf{z}{h}, \boldsymbol{\zeta}{h}), r{h}) |_{h} \

\sup{\boldsymbol{0} \neq ((\boldsymbol{v}{h}, \boldsymbol{\theta}{h}), q{h}) \in [\mathbf{V}{h} \times \mathbf{W}{h}] \times \mathbf{Q}{h}} \frac{\mathcal{B}{h}^{\boldsymbol{u}{h}}(\boldsymbol{u} - \boldsymbol{u}{h}, \boldsymbol{\omega} - \boldsymbol{\omega}{h}, p - p{h}; \boldsymbol{v}{h}, \boldsymbol{\theta}{h}, q{h})}{| ((\boldsymbol{v}{h}, \boldsymbol{\theta}{h}), q{h}) |{h}} \ = | ((\boldsymbol{u}, \boldsymbol{\omega}), p) - ((\mathbf{z}{h}, \boldsymbol{\zeta}{h}), r{h}) |_{h} \
\sup{\boldsymbol{0} \neq ((\boldsymbol{v}{h}, \boldsymbol{\theta}{h}), q{h}) \in [\mathbf{V}{h} \times \mathbf{W}{h}] \times \mathbf{Q}{h}} \frac{\mathcal{B}{h}^{\boldsymbol{u}{h}}(\boldsymbol{u}, \boldsymbol{\omega}, p; \boldsymbol{v}{h}, \boldsymbol{\theta}{h}, q{h}) - F{h}(\boldsymbol{v}{h}, \boldsymbol{\theta}{h})}{| ((\boldsymbol{v}{h}, \boldsymbol{\theta}{h}), q{h}) |_{h}}. \tag{5.3}$$

Finally, (5.1) is obtained directly from estimate (5.3) along with the error decomposition (cf. (5.2)) and the triangle inequality. □

The first term on the right-hand side of (5.1) measures the approximation property of [V^h ×Wh]× Q^h with respect to the norm ∥((•, •), •)∥h, while the second term captures the consistency error arising from the nonconforming discretisation.

Let ((u, ω), p) be the solution of the continuous problem. Similarly to [15, Lemma 3.2], assuming that ((u, ω), p) ∈ (V × H¹ (Ω)) × H¹ (Ω), from (2.7), integrating by parts, using the strong form of the momentum balance (first equation in (2.2)), and applying some algebraic manipulations, we get

$$\begin{aligned} &|\mathcal{B}{h}^{\boldsymbol{u}{h}}(\boldsymbol{u},\boldsymbol{\omega},\boldsymbol{p};\boldsymbol{v}{h},\boldsymbol{\theta}{h},q{h}) - F{h}(\boldsymbol{v}{h},\boldsymbol{\theta}{h})| \ &= |\mathcal{B}{h}^{\boldsymbol{u}{h}}(\boldsymbol{u},\boldsymbol{\omega},\boldsymbol{p};\boldsymbol{v}{h},\boldsymbol{\theta}{h},q{h}) - \int{\Omega} (\kappa^{-1}\boldsymbol{u} + \sqrt{\nu}\mathbf{curl}\,\boldsymbol{\omega} + \mathsf{F}|\boldsymbol{u}|\,\boldsymbol{u} + \nabla\boldsymbol{p} + \frac{1}{\sqrt{\nu}}\,\boldsymbol{\omega} \times \boldsymbol{u}) \cdot \mathcal{I}^{\mathsf{RT}}\boldsymbol{v}{h}| \ &\leq \left| \sqrt{\nu} \sum{K \in \mathcal{T}{h}} \int{K} \mathbf{curl}\,\boldsymbol{\omega} \cdot (\boldsymbol{v}{h} - \mathcal{I}^{\mathsf{RT}}\boldsymbol{v}{h}) + \sum{K \in \mathcal{T}{h}} \int{K} \nabla\boldsymbol{p} \cdot (\boldsymbol{v}{h} - \mathcal{I}^{\mathsf{RT}}\boldsymbol{v}{h}) \right| \ &+ \left| \frac{1}{\sqrt{\nu}} \int{\Omega} ([\boldsymbol{u} - \boldsymbol{u}{h}] \times \boldsymbol{\omega}) \cdot \mathcal{I}^{\mathsf{RT}}\boldsymbol{v}{h} + \mathsf{F} \int{\Omega} [|\boldsymbol{u}{h}| - |\boldsymbol{u}|] \,\boldsymbol{u} \cdot \mathcal{I}^{\mathsf{RT}}\boldsymbol{v}_{h} \right|. \end{aligned} \tag{5.4}$$

Now, we are ready to determine the order of convergence of the proposed method.

Theorem 5.2 (Rate of convergence) Let ((u, ω), p) and ((uh, ωh), ph) solve the continuous and discrete problems (2.7) and (4.8), respectively. Assume that the data satisfies

$$C{\rm S}^2 C{\rm RT} \frac{8}{\gamma \tilde{\gamma}} \left( \frac{1}{\sqrt{\nu}} + \mathsf{F} \right) | \boldsymbol{f} |_{0,\Omega} \le 1 \tag{5.5}$$

and that ((u, ω), p) ∈ (H² 0 (Ω) × H¹ (Ω)) × H¹ (Ω). Then there exists Crate > 0, independent of h, such that

$$|((\boldsymbol{u},\boldsymbol{\omega}),p)-((\boldsymbol{u}_h,\boldsymbol{\omega}_h),p_h)|h \le C{\text{rate}} h { |\boldsymbol{u}|{2,\Omega} + |\boldsymbol{\omega}|{1,\Omega} + |p|_{1,\Omega} }.$$ (5.6)

Proof. To prove the result, we must estimate the two terms on the right-hand side of (5.1). For the first term, we can use the properties in (4.6) (to estimate the errors for ω and p) and (4.5a) (to estimate the error for u), while for the second term, we apply estimate (5.4), along with interpolation properties (4.4) and (4.5b), and the Cauchy–Schwarz inequality, to obtain

$$\begin{split} |((\boldsymbol{u},\boldsymbol{\omega}),p)-((\boldsymbol{u}_h,\boldsymbol{\omega}_h),p_h)|h &\leq (C{\mathrm{CR}}+2C{\mathrm{P}})\Big(1+\frac{2}{\widetilde{\gamma}}\Big)\,h\Big{|\boldsymbol{u}|{2,\Omega}\,+\,|\boldsymbol{\omega}|{1,\Omega}\,+\,|p|{1,\Omega}\Big}\ &\quad +CF(\sqrt{\nu}+1)\frac{2}{\widetilde{\gamma}}h\Big{|\mathbf{curl}\,\boldsymbol{\omega}|{0,\Omega}\,+\,|\nabla p|{0,\Omega}\Big}\ &\quad +C{\mathrm{S}}^2C{\mathrm{RT}}\frac{2}{\widetilde{\gamma}}(\frac{1}{\sqrt{\nu}}+\mathsf{F})(|\boldsymbol{u}|{0,\Omega}+|\boldsymbol{\omega}|_{0,\Omega})|\boldsymbol{u}-\boldsymbol{u}_h|_h. \end{split}$$

From the latest estimate, using the fact that (u, ω) satisfies estimate (3.26) and applying hypothesis (5.5), we obtain (5.6). □

The estimate above can be refined to reflect the pressure-robustness of the formulation. Consider the continuous and discrete problems in their reduced form (in the continuous and discrete kernels V⁰ and V⁰ h , respectively)

$$a(\boldsymbol{u}, \boldsymbol{\omega}; \boldsymbol{v}, \boldsymbol{\theta}) + c^{\boldsymbol{u}}(\boldsymbol{u}, \boldsymbol{\omega}; \boldsymbol{v}, \boldsymbol{\theta}) = F(\boldsymbol{v}, \boldsymbol{\theta}) \quad \forall (\boldsymbol{v}, \boldsymbol{\theta}) \in \mathbf{V}_0 \times \mathbf{W},$$ (5.7a)

and

$$a_h(\boldsymbol{u}_h, \boldsymbol{\omega}_h; \boldsymbol{v}_h, \boldsymbol{\theta}_h) + c_h^{\boldsymbol{u}_h}(\boldsymbol{u}_h, \boldsymbol{\omega}_h; \boldsymbol{v}_h, \boldsymbol{\theta}_h) = F_h(\boldsymbol{v}_h, \boldsymbol{\theta}_h) \quad \forall (\boldsymbol{v}_h, \boldsymbol{\theta}_h) \in \mathbf{V}_h^0 \times \mathbf{W}_h, \tag{5.7b}$$

which are equivalent to (2.7) and (4.8), respectively.

In order to show a pressure-robust refinement of the previous results, we require an auxiliary bound regarding the piecewise norm control in the discrete norm ∥ • ∥^h of vector fields from V + Vh. In turn, for this as well as for the a posteriori error estimation later on, we will employ the companion operator

$$\mathcal{J} \in \mathcal{L}(\mathbf{V}_h; \mathbf{V}) \tag{5.8}$$

that is a right-inverse of the Crouzeix–Raviart interpolation I CR ∈ L(V; Vh), satisfying

$$|h_{\mathcal{T}_h}(\boldsymbol{w}_h - \mathcal{J}\boldsymbol{w}h)|{0,\Omega} \lesssim |\boldsymbol{w}_h - \mathcal{J}\boldsymbol{w}h|{pw},$$

as well as other additional L ² orthogonality properties not needed herein (see the precise design for 2D and 3D in [22, 26]).

Lemma 5.3 For any v ∈ V and w^h ∈ V^h and 1 ≤ s ≤ 6 in 3D (and 1 ≤ s < ∞ in 2D), there holds

$$|\mathbf{v} + \mathbf{w}h|{\mathbf{L}^s(\Omega)} + |\mathbf{v} + \mathbf{w}h|{\text{pw}} \le C_{\sharp} |\mathbf{v} + \mathbf{w}_h|_h$$ (5.9)

with C^♯ depending on Ω and the shape-regularity of the mesh T^h (and as well on s in the 2D case).

Proof. We first add ±J w^h to the left-hand side of (5.9), use triangle inequality, and invoke the well-known discrete Sobolev embedding with constant CdS(s) > 0:

$$|\boldsymbol{w}|{\mathbf{L}^{s}(\Omega)} \le C{\mathrm{dS}}(s) |\boldsymbol{w}|_{\mathrm{pw}} \qquad \forall \boldsymbol{w} \in \mathbf{V}$$

applied on the term v − J w^h ∈ V. Then, we employ the continuous Sobolev embedding (3.3), and this gives

$$LHS_{(5.9)} := | \boldsymbol{v} + \boldsymbol{w}h |{\mathbf{L}^s(\Omega)} + | \boldsymbol{v} + \boldsymbol{w}h |{pw} \le C_{dS}(s) | \boldsymbol{v} + \mathcal{J} \boldsymbol{w}h |{pw} + | \boldsymbol{w}_h - \mathcal{J} \boldsymbol{w}h |{\mathbf{L}^s(\Omega)} + | \boldsymbol{v} + \boldsymbol{w}h |{pw}$$

$$\begin{split} & \leq (1 + C_{\mathrm{dS}}(s)) |\boldsymbol{v} + \mathcal{J}\boldsymbol{w}h|{\mathrm{pw}} + |\boldsymbol{w}_h - \mathcal{J}\boldsymbol{w}h|{\mathbf{L}^s(\Omega)} + |\boldsymbol{w}_h - \mathcal{J}\boldsymbol{w}h|{\mathrm{pw}} \ & \leq (1 + C_{\mathrm{dS}}(s)) |\boldsymbol{v} + \mathcal{J}\boldsymbol{w}h|{\mathrm{pw}} + (1 + C_{\mathrm{S}}) |\boldsymbol{w}_h - \mathcal{J}\boldsymbol{w}h|{\mathrm{pw}} \ & \leq (1 + C_{\mathrm{dS}}(s)) [|\mathbf{curl}(\boldsymbol{v} + \mathcal{J}\boldsymbol{w}h)|{0,\Omega} + |\mathrm{div}(\boldsymbol{v} + \mathcal{J}\boldsymbol{w}h)|{0,\Omega}] + (1 + C_{\mathrm{S}}) |\boldsymbol{w}_h - \mathcal{J}\boldsymbol{w}h|{\mathrm{pw}} \end{split}$$

with the equivalence — valid for v + J w^h ∈ V — between the piecewise norm ||| • |||pw and the semi-norm in H0(curl , Ω) ∩ H0(div, Ω) [31, Lemma 2.5 & Remark 2.7] in the last step.

Using next again triangle inequality and the definition of the broken curl, the broken divergence, and the discrete velocity norm ∥ • ∥h, from the bounds above we readily get

$$\begin{split} \text{LHS}{(5.9)} & \leq [(1 + C{\text{dS}}(s))C{\text{norm}} + (1 + C{\text{S}})] |\boldsymbol{w}_h - \mathcal{J}\boldsymbol{w}h|{\text{pw}} \ & + (1 + C_{\text{dS}}(s))[|\text{curl}_h(\boldsymbol{v} + \boldsymbol{w}h)|{0,\Omega} + |\text{div}_h(\boldsymbol{v} + \boldsymbol{w}h)|{0,\Omega}] \ & \leq C_{\sharp} |\boldsymbol{v} + \boldsymbol{w}_h|_h \end{split}$$

with the following estimate from [26] in the last step:

$$|\boldsymbol{w}_h - \mathcal{J}\boldsymbol{w}h|{\text{pw}} \lesssim |\boldsymbol{v} + \boldsymbol{w}h|{\text{pw}},$$

as well as the fact that

$$|\mathbf{curl}_h(\boldsymbol{v}+\boldsymbol{w}h)|{0,\Omega} + |\mathrm{div}_h(\boldsymbol{v}+\boldsymbol{w}h)|{0,\Omega} \lesssim |\boldsymbol{v}+\boldsymbol{w}h|{\mathrm{pw}}.$$

Therefore C^♯ > 0 depends on CS, CdS(s) and on the parameter-dependent constant Cnorm > 0. □ Lemma 5.3 implies, in particular, that

$$|\boldsymbol{u} - \boldsymbol{u}h|{\mathbf{L}^4(\Omega)} \le C_{\sharp} |\boldsymbol{u} - \boldsymbol{u}_h|_h. \tag{5.10}$$

Theorem 5.4 (Pressure-robust error bound) Assume that (u, ω),(uh, ωh) are the unique solutions to (5.7a) and (5.7b), respectively. Suppose further that the data satisfies (5.5) with r = 1−C 2 ^SCRT γγe √ 1 +F ∥f∥0,^Ω > 0, and that the continuous vorticity is more regular ω ∈ H¹ (Ω). Then

$$|\boldsymbol{u} - \boldsymbol{u}_h|_h + |\boldsymbol{\omega} - \boldsymbol{\omega}h|{0,\Omega} \le Ch|\boldsymbol{\omega}|{1,\Omega} + \inf{\substack{(\boldsymbol{v}_h,\boldsymbol{\theta}_h) \in \mathbf{V}_h^0 \times \mathbf{W}_h \ \alpha_h|\boldsymbol{v}_h|_h = 1}} \left\frac{1}{r} + \frac{1}{r\alpha_h}\right. \quad (5.11)$$

Proof. Let us adopt the notation ⃗u := (u, ω), ⃗v := (v, θ), ⃗z := (z, ζ) ∈ V⁰ × W and similarly for their discrete counterparts ⃗uh, ⃗vh, ⃗z^h ∈ V⁰ ^h × Wh, denoting the corresponding discrete norm as, e.g., ∥⃗zh∥ := ∥zh∥^h + ∥ζh∥0,Ω. First we decompose

$$\vec{u} - \vec{u}_h = \vec{u} - \vec{v}_h + \vec{v}_h - \vec{u}_h = \vec{u} - \vec{v}_h - \vec{y}_h$$

and note that ⃗y^h := ⃗u^h − ⃗v^h belongs to V⁰ ^h × Wh. Then, using the inf-sup condition of ah(•, •), the definition of ⃗yh, and the definition of the discrete problem (5.7b), we can write the following

$$\alpha{h} |\vec{z}{h}| |\vec{y}{h}| \leq \sup{0 \neq \vec{z}{h} \in \mathbf{V}{h}^{0} \times \mathbf{W}{h}} a{h}(\vec{y}{h}, \vec{z}{h}) = \sup{0 \neq \vec{z}{h} \in \mathbf{V}{h}^{0} \times \mathbf{W}{h}} a{h}(\vec{u}{h} - \vec{v}{h}, \vec{z}{h}) = \sup{0 \neq \vec{z}{h} \in \mathbf{V}{h}^{0} \times \mathbf{W}{h}} \left( a{h}(\vec{u} - \vec{v}{h}, \vec{z}{h}) - a{h}(\vec{u}, \vec{z}{h}) + a{h}(\vec{u}{h}, \vec{z}{h}) \right) = \sup{0 \neq \vec{z}{h} \in \mathbf{V}{h}^{0} \times \mathbf{W}{h}} \left( a{h}(\vec{u} - \vec{v}{h}, \vec{z}{h}) + \left[ F{h}(\vec{z}{h}) - c{h}^{u{h}}(\vec{u}{h}, \vec{z}{h}) - a{h}(\vec{u}, \vec{z}_{h}) \right] \right).$$ (5.12)

Consequently, dividing through αh∥⃗zh∥ and using the boundedness of ah(•, •) as well as the triangle inequality ∥⃗u − ⃗uh∥ ≤ ∥⃗u − ⃗vh∥ + ∥⃗yh∥, we obtain the following estimate

$$|\vec{u} - \vec{u}h| \le \inf{\vec{v}_h \in \mathbf{V}_h^0 \times \mathbf{W}_h} \left( 1 + \frac{1}{\alpha_h} \right) |\vec{u} - \vec{v}h| + \sup{\vec{z}_h \in \mathbf{V}_h^0 \times \mathbf{W}_h} \frac{\left| F_h(\vec{z}_h) - c_h^{\mathbf{u}_h}(\vec{u}_h, \vec{z}_h) - a_h(\vec{u}, \vec{z}_h) \right|}{\alpha_h |\vec{z}_h|}$$ (5.13)

composed by the best approximation in V⁰ ^h × W^h and the consistency error.

Note that even if I RTz^h is not in H¹ (Ω), we have

$$|\mathcal{I}^{\mathsf{RT}}\mathbf{z}h|{\mathbf{L}^4(\Omega)} \leq |\mathcal{I}^{\mathsf{RT}}\mathbf{z}_h - \mathbf{z}h|{\mathbf{L}^4(\Omega)} + |\mathbf{z}h|{\mathbf{L}^4(\Omega)} \leq C|\mathbf{z}h|{\mathbf{L}^4(\Omega)} \leq C|\nabla_h\mathbf{z}h|{0,\Omega} \leq C_*|\mathbf{z}_h|_h, \quad (5.14)$$

thanks to triangle inequality, inverse estimates, and discrete Sobolev properties [42, Theorem 4.12].

Next, we proceed to add ±c u h (⃗u, ⃗zh) to the numerator in the consistency error ah(⃗u, ⃗zh)+c u^h h (⃗uh, ⃗zh)− Fh(⃗zh). The + will contribute to complete a full residual ah(⃗u, ⃗zh) + c u h (⃗u, ⃗zh) − Fh(⃗zh), so we need to investigate first the remainder terms as follows, adding and subtracting appropriate terms. Applying H¨older's inequality, property (5.14), and reverse triangle inequality, gives

$$\begin{aligned} |c{h}^{\mathbf{u}{h}}(\vec{u}{h}, \vec{z}{h}) - c{h}^{\mathbf{u}}(\vec{u}, \vec{z}{h})| &= \left| -\frac{1}{\sqrt{\nu}} \int{\Omega} \mathbf{u}{h} \times (\boldsymbol{\omega}{h} - \boldsymbol{\omega}) \cdot \mathcal{I}^{\mathsf{RT}} \mathbf{z}{h} + \mathsf{F} \int{\Omega} |\mathbf{u}{h}| (\mathbf{u}{h} - \boldsymbol{u}) \cdot \mathcal{I}^{\mathsf{RT}} \mathbf{z}{h} \right| \ &+ \frac{1}{\sqrt{\nu}} \int{\Omega} (\mathbf{u} - \mathbf{u}{h}) \times \boldsymbol{\omega} \cdot \mathcal{I}^{\mathsf{RT}} \mathbf{z}{h} - \mathsf{F} \int{\Omega} [|\mathbf{u}| - |\mathbf{u}{h}|] \mathbf{u} \cdot \mathcal{I}^{\mathsf{RT}} \mathbf{z}{h} | \ &\leq \frac{C{*}^{2}}{\sqrt{\nu}} |\mathbf{u}{h}|{h} |\boldsymbol{\omega} - \boldsymbol{\omega}{h}|{0,\Omega} |\mathbf{z}{h}|{h} + \mathsf{F} C{}^{2} |\mathbf{u}{h}|{h} |\mathbf{u} - \mathbf{u}{h}|{0,\Omega} |\mathbf{z}{h}|{h} \ &+ \frac{C_{}}{\sqrt{\nu}} |\mathbf{u} - \mathbf{u}{h}|{\mathbf{L}^{4}(\Omega)} |\boldsymbol{\omega}|{0,\Omega} |\mathbf{z}{h}|{h} + \mathsf{F} C{} |\mathbf{u} - \mathbf{u}{h}|{\mathbf{L}^{4}(\Omega)} |\mathbf{u}|{0,\Omega} |\mathbf{z}{h}|{h} \ &\lesssim C{} [M{h}C{*} + MC{\sharp}] \left(\mathsf{F} + \frac{1}{\sqrt{\nu}}\right) (|\mathbf{u} - \mathbf{u}{h}|{h} + |\boldsymbol{\omega} - \boldsymbol{\omega}{h}|{0,\Omega}) |\mathbf{z}{h}|_{h} \end{aligned} \tag{5.15}$$

where for the last estimation we have used that ⃗u, ⃗u^h are solutions to the continuous and discrete problems featuring a continuous dependence on data that we denote here by M = 2 γ ∥f∥0,^Ω and M^h = 2 ^γeCF^h ∥f∥0,Ω, respectively (cf. (3.28) and (4.29), respectively); we have also used (5.10).

We now look again at the numerator of the consistency error and rewrite f in terms of the left-hand side of the momentum balance equation in (2.2), use the fact that

$$\int_{\Omega} \nabla p \cdot \mathcal{I}^{\mathsf{RT}} \mathbf{z}_h = 0$$

(see, e.g., [40]) and apply integration by parts on the term ^√ ν R Ω curl ^hz^h · ω, to get

$$\begin{split} a_h(\vec{u}, \vec{z}_h) + c_h^{\boldsymbol{u}_h}(\vec{u}_h, \vec{z}_h) - F_h(\vec{z}_h) + \left[c_h^{\boldsymbol{u}_h}(\vec{u}_h, \vec{z}_h) - c_h^{\boldsymbol{u}}(\vec{u}, \vec{z}h)\right] \ &= \sqrt{\nu} \int{\Omega} (\mathbf{curl} \, \boldsymbol{u} - \boldsymbol{\omega}) \cdot \boldsymbol{\zeta}h + \sqrt{\nu} \int{\Omega} \mathbf{curl} \, \boldsymbol{\omega} \cdot (\mathbf{z}_h - \mathcal{I}^{\mathsf{RT}} \mathbf{z}_h) + \left[c_h^{\boldsymbol{u}_h}(\vec{u}_h, \vec{z}_h) - c_h^{\boldsymbol{u}}(\vec{u}, \vec{z}h)\right] \ &\leq \sqrt{\nu} |\mathbf{curl} \, \boldsymbol{\omega}|{0,\Omega} |\mathbf{z}_h - \mathcal{I}^{\mathsf{RT}} \mathbf{z}h|{0,\Omega} \ &+ C_* [Mh C* + M C_\sharp] \left(\mathsf{F} + \frac{1}{\sqrt{\nu}}\right) (|\boldsymbol{u} - \boldsymbol{u}_h|_h + |\boldsymbol{\omega} - \boldsymbol{\omega}h|{0,\Omega}) |\mathbf{z}_h|h \ &\leq \left[ch|\boldsymbol{\omega}|{1,\Omega} + C_* [Mh C* + M C_\sharp] \left(\mathsf{F} + \frac{1}{\sqrt{\nu}}\right) |\vec{u} - \vec{u}h|{1,\Omega} \right] |\mathbf{z}_h|_h. \end{split}$$

Here we have used the estimate in (5.15) and the fact that if ⃗u is the exact smooth solution of (5.7a) then ^√ νcurlu = ω in Ω, and have also used approximation properties of I RT. This yields

$$\sup_{\vec{z}_h \in \mathbf{V}_h^0 \times \mathbf{W}_h} \frac{\left| F_h(\vec{z}_h) - c_h^{\mathbf{u}_h}(\vec{u}_h, \vec{z}_h) - a_h(\vec{u}, \vec{z}_h) \right|}{\alpha_h |\vec{z}h|} \leq Ch |\boldsymbol{\omega}|{1,\Omega} + C_[Mh C + M C_{\sharp}] \left(\mathsf{F} + \frac{1}{\sqrt{\nu}}\right) |\vec{u} - \vec{u}_h|.$$

And the proof is complete after combining this estimate, the small data assumption (5.5), and (5.14). □

6 A posteriori error analysis

This section is devoted to a reliable and efficient computable error control. The derivation of the reliable error bound departs from the assumption that (u, ω, p) ∈ V × L 2 (Ω) × L 2 0 (Ω) solves (2.4) and (uh, ωh, ph) ∈ V^h × W^h × Q^h solves (4.8) to define the errors e^u := u − J uh, e^ω := ω − ωh, and e^p := p − p^h as in (5.2) except that e^u := u − J u^h ∈ V is not equal to u − u^h ∈ V + Vh.

Recall from (3.26) that u belongs to K and thereafter define the bilinear form B u •, • as in (3.19) (with ^u^b replaced by ^u). Recall (3.23) and deduce that there exists some test function (v, ^θ, q) ^∈ V × L 2 (Ω) × L 2 0 (Ω) of norm ∥((v, θ), q)∥ ≤ 1 at most one and

$$\frac{\gamma}{2} | (\mathbf{e}{\boldsymbol{u}}, \mathbf{e}{\boldsymbol{\omega}}, ep) | = \mathcal{B}^{\boldsymbol{u}} (\mathbf{e}{\boldsymbol{u}}, \mathbf{e}_{\boldsymbol{\omega}}, e_p; \boldsymbol{v}, \boldsymbol{\theta}, q).$$

The solution (u, ω, p) ∈ V×L (Ω)×L 2 0 (Ω) to (2.4) also satisfies B u u, ω, p; v, θ, q = F(v, θ), whence

$$\frac{\gamma}{2}|(\mathbf{e}{\bm{u}},\mathbf{e}{\bm{\omega}},e{\bm{p}})| = \int{\Omega} \bm{f}\cdot\bm{v} - \mathcal{B}^{\bm{u}}\big(\mathcal{J}\bm{u}_h,\bm{\omega}_h,p_h;\bm{v},\bm{\theta},q\big).$$

The right-hand side in this identity defines a residual in terms of the test functions v, θ, and q. The discrete solution (uh, ωh, ph) ∈ Vh×Wh×Q^h to (4.8) involves the discrete right-hand side from (4.7d) and the operator I RT from (4.2). With the definition of the Crouzeix–Raviart interpolation v^h := I CRv and the piecewise constant integral means ω^h and p^h of ω and p, respectively, we investigate the first identity of (4.8) for θ^h = 0 and (4.7a)-(4.7c), namely

$$\begin{split} \int_{\Omega} \boldsymbol{f} \cdot \boldsymbol{\mathcal{I}}^{\mathsf{RT}} \boldsymbol{v}h &= \frac{1}{\kappa} \int{\Omega} \boldsymbol{u}h \cdot \boldsymbol{\mathcal{I}}^{\mathsf{RT}} \boldsymbol{v} + \sum{F \in \mathcal{E}_h^{\mathsf{int}}} \frac{\vartheta}{h_F} \int_F \left( \nu \llbracket \boldsymbol{u}_h \times \mathbf{n} \rrbracket_F \cdot \llbracket \boldsymbol{v}_h \times \mathbf{n} \rrbracket_F + \llbracket \boldsymbol{u}_h \cdot \mathbf{n} \rrbracket_F \llbracket \boldsymbol{v}_h \cdot \mathbf{n} \rrbracketF \right) \ &+ \int{\Omega} \left( \sqrt{\nu} \ \boldsymbol{\omega}_h \cdot \mathbf{curl} \, \boldsymbol{v} - ph \mathrm{div} \, \boldsymbol{v} \right) + \int{\Omega} \left( \mathsf{F} |\boldsymbol{u}_h| \, \boldsymbol{u}_h - \nu^{-1/2} \boldsymbol{u}_h \times \boldsymbol{\omega}_h \right) \cdot \boldsymbol{\mathcal{I}}^{\mathsf{RT}} \boldsymbol{v} \end{split}$$

with R K ( ν ω^h · curl v^h − phdiv vh) = R K ( ν ω^h · curl v − phdiv v) for all K ∈ T^h from the integral mean property of the gradients for the Crouzeix–Raviart interpolation in the last step. The combination of the last two identities and the definitions (3.19) and (3.10) reveal the key identity

$$\begin{split} \frac{\gamma}{2} | (\mathbf{e}{\boldsymbol{u}}, \mathbf{e}{\boldsymbol{\omega}}, ep) | &= \int{\Omega} \boldsymbol{f} \cdot (\boldsymbol{v} - \mathcal{I}^{\mathsf{RT}} \boldsymbol{v}) + \kappa^{-1} \int_{\Omega} (\boldsymbol{u}_h \cdot \mathcal{I}^{\mathsf{RT}} \boldsymbol{v} - (\mathcal{J} \boldsymbol{u}h) \cdot \boldsymbol{v}) \ &+ \sum{F \in \mathcal{E}_h^{\mathrm{int}}} \frac{\vartheta}{hF} \int{F} (\nu \llbracket \boldsymbol{u}_h \times \mathbf{n} \rrbracket_F \cdot \llbracket \boldsymbol{v}_h \times \mathbf{n} \rrbracket_F + \llbracket \boldsymbol{u}_h \cdot \mathbf{n} \rrbracket_F \llbracket \boldsymbol{v}_h \cdot \mathbf{n} \rrbracketF) \ &+ \int{\Omega} (\mathsf{F} |\boldsymbol{u}_h| \, \boldsymbol{u}_h - \nu^{-1/2} \boldsymbol{u}_h \times \boldsymbol{\omega}h) \cdot \mathcal{I}^{\mathsf{RT}} \boldsymbol{v} + \int{\Omega} (\nu^{-1/2} \boldsymbol{u} \times \boldsymbol{\omega}_h - \mathsf{F} |\boldsymbol{u}| \, \mathcal{J} \boldsymbol{u}h) \cdot \boldsymbol{v} \ &- \int{\Omega} \left( \nu^{1/2} \boldsymbol{\theta} \cdot \mathbf{curl} \, \mathcal{J} \boldsymbol{u}_h - \boldsymbol{\omega}h \cdot \boldsymbol{\theta} \right) + \int{\Omega} q \, \operatorname{div} \mathcal{J} \boldsymbol{u}_h. \end{split}$$

The jump terms with [[uh]]^F = [[uh−J uh]]^F and [[vh]]^F = [[vh−v]]^F allow for standard trace inequalities and interpolation local interpolation error estimates and eventually verify

$$\sum_{F \in \mathcal{E}_h^{\mathrm{int}}} \frac{\vartheta}{h_F} \int_F \left( \nu \llbracket \boldsymbol{u}_h \times \mathbf{n} \rrbracket_F \cdot \llbracket \boldsymbol{v}_h \times \mathbf{n} \rrbracket_F + \llbracket \boldsymbol{u}_h \cdot \mathbf{n} \rrbracket_F \llbracket \boldsymbol{v}_h \cdot \mathbf{n} \rrbracket_F \right) \lesssim (1 + \nu) | \boldsymbol{u}_h - \mathcal{J} \boldsymbol{u}h |{\mathrm{pw}}$$

with |||v − vh|||pw ≲ 1 in the last step. The remaining terms that involve a factor v, v^h = I CRv, or I RTv^h = I RTv combine to one residual term plus perturbations. The residual reads

$$\int_{\Omega} \left( \boldsymbol{f} - \kappa^{-1} \boldsymbol{u}_h + \nu^{-1/2} \boldsymbol{u}_h \times \boldsymbol{\omega}_h - \mathsf{F} |\boldsymbol{u}_h| \, \boldsymbol{u}h \right) \cdot \left( \boldsymbol{v} - \boldsymbol{\mathcal{I}}^{\mathsf{RT}} \boldsymbol{v} \right) + \kappa^{-1} \int{\Omega} (\boldsymbol{u}_h - \boldsymbol{\mathcal{J}} \boldsymbol{u}_h) \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot \boldsymbol{v}_h \cdot$$

and the remaining perturbations read

$$\nu^{-1/2}\int_{\Omega}\left(\left(\boldsymbol{u}-\boldsymbol{u}h\right)\times\boldsymbol{u}+\int{\Omega}\mathsf{F}\big(\left(\left|\boldsymbol{u}_h\right|-\left|\boldsymbol{u}\right|\big)\big)\biguplus_h+\left|\boldsymbol{u}\left|\left(\boldsymbol{u}_h-\mathcal{J}\boldsymbol{u}_h\right)\right)\cdot v$$

$$\leq \left[ \nu^{-1/2} |\boldsymbol{\omega}h|{0,\Omega} + \mathsf{F} |\boldsymbol{u}h|{0,\Omega} \right] |\boldsymbol{u} - \boldsymbol{u}h|{\mathbf{L}^4(\Omega)} |\boldsymbol{v}|{\mathbf{L}^4(\Omega)} + \mathsf{F} |\boldsymbol{u}|{\mathbf{L}^4(\Omega)} |\boldsymbol{v}|_{\mathbf{L}^4(\Omega)} |\boldsymbol{u}_h - \mathcal{J}\boldsymbol{u}h|{0,\Omega}$$

$$\leq C_{dS}(4)C_S(4) \left[ \nu^{-1/2} |\boldsymbol{\omega}h|{0,\Omega} + \mathsf{F} |\boldsymbol{u}h|{0,\Omega} \right] |\boldsymbol{u} - \boldsymbol{u}h|{\mathrm{pw}} + C_{dS}(4)CS(4)\mathsf{F} |\boldsymbol{u}|{1,\Omega} h_{\mathrm{max}} C_I |\boldsymbol{u}_h - \mathcal{J}\boldsymbol{u}h|{\mathrm{pw}}$$

with discrete Sobolev (resp. Sobolev) inequalities with constant CdS(4) ≲ 1 (resp. CS(4) ≲ 1), interpolation error estimates, and |v|1,^Ω ≤ 1 in the last step.

The discrete equations also reveal ah(uh, ωh; 0, θh) = 0 = bh(uh, ωh; qh) for all piecewise constant θ^h ∈ W^h and for all piecewise constant q^h ∈ Q^h with integral mean zero over the domain. Those identities localise (utilise R Ω divhu^h = 0 from a piecewise integration by parts to overcome the global constraint in Qh) and lead to the discrete identities

$$\sqrt{\nu}\operatorname{\mathbf{curl}}_h \boldsymbol{u}_h = \boldsymbol{\omega}_h \quad \text{and} \quad \operatorname{div}_h \boldsymbol{u}_h = 0$$

for the piecewise constant functions (from piecewise action of the differential operators). Hence

$$\begin{split} &-\int{\Omega}\left(\sqrt{\nu}\ \boldsymbol{\theta}\cdot\mathbf{curl}\ \mathcal{J}\boldsymbol{u}{h}-\boldsymbol{\omega}{h}\cdot\boldsymbol{\theta}\right)+\int{\Omega}q\ \mathrm{div}\mathcal{J}\boldsymbol{u}{h}\ &=\int{\Omega}\sqrt{\nu}\ \boldsymbol{\theta}\cdot\mathbf{curl}\ {h}(\boldsymbol{u}{h}-\mathcal{J}\boldsymbol{u}{h})-\int{\Omega}q\ \mathrm{div}{h}(\boldsymbol{u}{h}-\mathcal{J}\boldsymbol{u}{h})\ &\leq\left(\sqrt{\nu}|\boldsymbol{\theta}|+\sqrt{2}|q|\right)\ |!|\boldsymbol{u}{h}-\mathcal{J}\boldsymbol{u}{h}|{\mathrm{pw}}\leq\sqrt{2+\nu}\ |!|\boldsymbol{u}{h}-\mathcal{J}\boldsymbol{u}{h}|_{\mathrm{pw}} \end{split}$$

with Cauchy inequalities and ∥((v, θ), q)∥ ≤ 1 in the last steps.

The combination of all the above estimates and the standard error estimate ∥h −1 T (v − IRTv)∥0,^Ω ≲ |||vh|||pw ≤ 1 reveal that

$$|(\mathbf{e}{\boldsymbol{u}}, \mathbf{e}{\boldsymbol{\omega}}, e{p})| - 2C{dS}(4)C{S}(4)\gamma^{-1} \left(\nu^{-1/2} |\boldsymbol{\omega}{h}|{0,\Omega} + \mathsf{F}|\boldsymbol{u}{h}|{0,\Omega}\right) |\boldsymbol{u} - \boldsymbol{u}{h}|_{pw}$$

$$\lesssim (1 + \nu^{1/2} + (\mathsf{F} + \kappa^{-1})h{\max}) |\boldsymbol{u}{h} - \mathcal{J}\boldsymbol{u}{h}|{pw}$$

$$+ |h{\mathcal{T}}(\boldsymbol{f} - \kappa^{-1}\boldsymbol{u}{h} + \nu^{-1/2}\boldsymbol{u}{h} \times \boldsymbol{\omega}{h} - \mathsf{F}|\boldsymbol{u}{h}|\,\boldsymbol{u}{h})|_{0,\Omega}.$$

(The notation ≲ includes generic constants as well as γ ≈ 1, while we keep γ in the rather explicit negative term of the lower bound.) Since e^u = u − J uh, the triangle inequality |||u − uh|||pw ≤ ∥(eu, eω, ep)∥ + |||u^h − J uh|||pw, provides

$$| \boldsymbol{u} - \boldsymbol{u}h |{\text{pw}} + | \mathbf{e}{\boldsymbol{\omega}} |{0,\Omega} + | ep |{0,\Omega} - 2C_{dS}(4)C_S(4)\gamma^{-1} \left( \nu^{-1/2} | \boldsymbol{\omega}h |{0,\Omega} + \mathsf{F} | \boldsymbol{u}h |{0,\Omega} \right) | \boldsymbol{u} - \boldsymbol{u}h |{\text{pw}}$$

$$\lesssim (1 + \nu^{1/2} + (\mathsf{F} + \kappa^{-1})h_{\text{max}}) | \boldsymbol{u}_h - \mathcal{J} \boldsymbol{u}h |{\text{pw}}$$

$$+ | h_{\mathcal{T}} (\boldsymbol{f} - \kappa^{-1} \boldsymbol{u}_h + \nu^{-1/2} \boldsymbol{u}_h \times \boldsymbol{\omega}_h - \mathsf{F} | \boldsymbol{u}_h | \boldsymbol{u}h ) |{0,\Omega}.$$

$$(6.1)$$

The proof concludes with a discussion of the factor CdS(4)CS(4) ν ⁻1/² ∥ωh∥0,^Ω + F∥uh∥0,^Ω ≤ γ/4 which follows from (4.29) for small sources f in L 2 (Ω). The point is that the latter estimate allows us to absorb the negative term on the lower bound of (6.2) and this leads to the a posteriori error control. The discrete Friedrichs inequality ∥uh∥0,^Ω ≤ CdF |||uh|||pw for all Crouzeix–Raviart functions with homogeneous boundary conditions and |||uh|||pw ≤ C♯∥uh∥^h by (5.9) provide

$$\nu^{-1/2} |\boldsymbol{\omega}h|{0,\Omega} + \mathsf{F} |\boldsymbol{u}h|{0,\Omega} \leq \nu^{-1/2} |\boldsymbol{\omega}h|{0,\Omega} + C{dF} C{\sharp} \mathsf{F} |\boldsymbol{u}_h|h \leq \sqrt{\nu^{-1} + C{dF}^2 C{\sharp}^2 \mathsf{F}^2} \frac{2}{\widetilde{\gamma}} C{Fh} |\boldsymbol{f}|{0,\Omega}$$

with a Cauchy inequality and (4.29) in the last step. Hence in case that

$$|\mathbf{f}|{0,\Omega} \le \gamma \widetilde{\gamma} / \left( 8C{Fh} C{dS}(4) CS(4) \sqrt{\nu^{-1} + C{dF}^2 C_{\sharp}^2 \mathsf{F}^2} \right),$$ (6.2)

the lower bound in (6.2) provides error control over ¹ 2 |||u − uh|||pw + ∥eω∥0,^Ω + ∥ep∥0,Ω. This concludes the proof of the a posteriori error estimate. The final form of the explicit residual-based a posteriori error estimate employs the well-established formula

$$|!|!|!|!| \boldsymbol{u}_h - \mathcal{J}\boldsymbol{u}h |!|!|!|^2{\mathrm{pw}} \lesssim \sum_{F \in \mathcal{E}} h_F^{-1} |!|!| [\boldsymbol{u}_h]_F \times \mathbf{n}_F |!|!|^2F \approx \sum{F \in \mathcal{E}} h_F |!|!| [\mathbf{D}\boldsymbol{u}_h]_F \times \mathbf{n}_F |!|!|^2_F.$$

This gives rise to the explicit residual-based a posteriori error estimator with the contribution

$$\eta^{2}(K) := |K|^{2/d} |\boldsymbol{f} - \kappa^{-1} \boldsymbol{u}{h} + \nu^{-1/2} \boldsymbol{u}{h} \times \boldsymbol{\omega}{h} - \mathsf{F}|\boldsymbol{u}{h}| \, \boldsymbol{u}{h}|{0,K}^{2} + |K|^{1/d} \sum{F \in \mathcal{F}(K)} |[\mathbf{D}\boldsymbol{u}{h}]{F} \times \mathbf{n}{F}|_{0,F}^{2} \quad (6.3)$$

for each triangle ^K ∈ T^h and the global version ^η(Th) := qP K∈T^h η ²(K).

Theorem 6.1 (A posteriori error control) Provided the source is sufficiently small such that (6.2) holds, we have reliability

$$|\mathbf{u} - \mathbf{u}h|{\text{pw}} + |\mathbf{e}{\boldsymbol{\omega}}|{0,\Omega} + |ep|{0,\Omega} \lesssim \eta(\mathcal{T}_h).$$

Efficiency holds even in local form up to data oscillations: For any K ∈ T^h with neighbourhood Ω(K) covered by the union of all simplices in T^h with zero distance to K, we have

$$\eta(K) \lesssim |D_{pw}(u - uh)|{0,\Omega(K)} + |u - uh|{0,K} + |\omega - \omegah|{0,K} + |ph - p|{0,K} + osc_k(f, K).$$

The generic multiplicative constants behind the notation ≲ exclusively depend on the material constants, upper bounds of the solutions u and u^h in L 2 (Ω), and the shape regularity of Th.

Proof. Reliability follows from the analysis prior to the statement of Theorem 6.1. The remaining efficiency follows with Verf¨urth's bubble-function methodology. This is well established for the side contributions

$$|K|^{1/(2d)} |[\mathbf{D} \boldsymbol{u}_h]_F \times \mathbf{n}F|{L^2(\partial K)} \lesssim |D_{pw}(\boldsymbol{u} - \boldsymbol{u}h)|{0,\Omega(K)}$$

and follows verbatim [26, 50]. The volume contribution, however, challenges with a technical (possibly unexpected) trap. The overall idea is to design a local test function v^K that allows an evaluation of the residual functional and thereby involves the exact solution. The data oscillations arise from the approximation of the source f in the very first step. In order to estimate the volume contribution

$$\mu_K := |K|^{1/d} | \boldsymbol{f} - \kappa^{-1} \boldsymbol{u}_h + \nu^{-1/2} \boldsymbol{u}_h \times \boldsymbol{\omega}_h - \mathsf{F} |\boldsymbol{u}_h| \, \boldsymbol{u}h |{0,K}$$

of K ∈ Th, we consider the volume-bubble function b^K ∈ W 1,∞ 0 (K) on a simplex K ⊂ R d . The latter is the product of all d + 1 barycentric coordinates of the vertices of K times a factor (d + 1)d+1 for the normalisation 0 ≤ b^K ≤ 1 = max b^K in K. With the abbreviation

$$\boldsymbol{v}_K := \Pi_1 \boldsymbol{f} - \kappa^{-1} \boldsymbol{u}_h + \nu^{-1/2} \boldsymbol{u}_h \times \boldsymbol{\omega}_h - \mathsf{F} |\boldsymbol{u}_h| \, \boldsymbol{u}_h,$$

the admissible test function reads bKv^K and belongs to V (for b^K is extended by zero). The typical application in Verf¨urth's bubble-function methodology considers only polynomial test functions that allow a standard inverse estimate ∥g∥0,K ≤ Ceq∥b 1/2 ^K g∥0,K for all g ∈ P1(K) d . In the application of this paper v^K is polynomial up to this extra term bKF|uh|u^h and hence requires a new inverse estimate. After surprisingly large efforts, Appendix A provides

$$|||f|f+g||{0,K} \le C{\text{eq}}||bK^{1/2}(|f|f+g)||{0,K} \text{ for all } f,g \in P_1(K)^d.$$

The test function bKv^K ∈ V therefore satisfies the novel nonstandard inverse estimate

$$|\boldsymbol{v}K|{0,K} \le C_{\text{eq}} |b_K^{1/2} \boldsymbol{v}K|{0,K}.$$

The remaining arguments in the efficiency proof are standard nowadays and (partly) apply verbatim. This, a triangle inequality, and the definition of the oscillations reveal

$$\mu_K \le |K|^{1/d} |\boldsymbol{f} - \Pi1 \boldsymbol{f}|{0,K} + \boldsymbol{v}K|{0,K} \le \operatorname{osc}k(\boldsymbol{f}, K) + C{\text{eq}} |K|^{1/d} |\boldsymbol{b}_K^{1/2} \boldsymbol{v}K|{0,K}.$$ (6.4)

With v replaced by bKv^K ∈ V in the first equation in (2.4), the last term relates to

$$|\boldsymbol{b}{K}^{1/2}\boldsymbol{v}{K}|{0,K}^{2} = \int{K} b{K}\boldsymbol{v}{K} \cdot (\Pi{1}\boldsymbol{f} - \kappa^{-1}\boldsymbol{u}{h} + \nu^{-1/2}\boldsymbol{u}{h} \times \boldsymbol{\omega}{h} - \mathsf{F}|\boldsymbol{u}{h}|\,\boldsymbol{u}{h})$$

$$= \int{K} b{K} \boldsymbol{v}{K} \cdot (\Pi{1} - 1) \boldsymbol{f} + \int{K} \left( \nu^{1/2} \boldsymbol{\omega} \cdot \operatorname{\mathbf{curl}} \left( b{K} \boldsymbol{v}{K} \right) - p \operatorname{div} (b{K} \boldsymbol{v}_{K}) \right)$$

$$+ \int{K} b{K} \boldsymbol{v}{K} \cdot \left( \kappa^{-1} (\boldsymbol{u} - \boldsymbol{u}{h}) - \nu^{-1/2} (\boldsymbol{u} \times \boldsymbol{\omega} - \boldsymbol{u}{h} \times \boldsymbol{\omega}{h}) + \mathsf{F} (|\boldsymbol{u}| \, \boldsymbol{u} - |\boldsymbol{u}{h}| \, \boldsymbol{u}{h}) \right).$$

$$(6.5)$$

The last right-hand side (6.5) consists of six summands we enumerate S1, . . . , S⁶ in the order displayed. Notice that there is a factor |K| ¹/d omitted compared to (6.4) and hence we can afford a factor |K| −1/d in the upper bounds of S1, . . . , S6. A Cauchy inequality and the definition of the oscillations provide

$$S_1 = \int_K b_K \boldsymbol{v}_K \cdot (\Pi_1 - 1) \boldsymbol{f} \le |b_K^{1/2} \boldsymbol{v}K|{0,K} |K|^{-1/d} \operatorname{osc}_k(\boldsymbol{f}, K).$$

Recall that bKv^K has a support on the simplex K and vanishes on its boundary ∂K. Hence the Gauss and Stokes theorems show R K div(bKvK) = 0 and R ^K ω · curl(bKvK) = 0. Since p^h and ω^h are constant on K, we infer

$$S_2 + S_3 = \nu^{1/2} \int_K (\boldsymbol{\omega} - \boldsymbol{\omega}_h) \cdot \mathbf{curl} (b_K \boldsymbol{v}_K) + \int_K (p_h - p) \operatorname{div} (b_K \boldsymbol{v}_K)$$

$$\lesssim (|\boldsymbol{\omega} - \boldsymbol{\omega}h|{0,K} + |ph - p|{0,K}) |b_K^{1/2} \boldsymbol{v}K|{0,K} |K|^{-1/d}$$

with an inverse estimate |bKvK|1,K ≲ ∥b 1/2 ^K vK∥0,K|K| ⁻1/d in the last step. A Cauchy inequality controls the term

$$S_4 = \int_K b_K \boldsymbol{v}_K \cdot \kappa^{-1} (\boldsymbol{u} - \boldsymbol{u}_h) \le \kappa^{-1} |b_K^{1/2} \boldsymbol{v}K|{0,K} |\boldsymbol{u} - \boldsymbol{u}h|{0,K}.$$

The nonlinear term S⁵ is related to ν ¹/2S⁵ ≤ ∥b 1/2 ^K vK∥L∞(K)∥(u−uh)×ω+uh×(ω−ωh)∥L1(K) . Since ω and u^h are bounded in L 2 (Ω) by the a priori error analysis, a Cauchy inequality for the L 1 integrals may be written as ∥(u − uh) × ω∥L1(K) ≲ ∥u − uh∥0,K and ∥u^h × (ω − ωh)∥L1(K) ≲ ∥ω − ωh∥0,K. This and the inverse estimate ∥b 1/2 ^K ^vK∥L∞(K) ^≲ ^|K^| ⁻1/d∥b 1/2 ^K vK∥0,K establish

$$S_5 \lesssim \nu^{-1/2} |K|^{-1/d} ||b_K^{1/2} \boldsymbol{v}K||{0,K} \left( ||\boldsymbol{u} - \boldsymbol{u}h||{0,K} + ||\boldsymbol{\omega} - \boldsymbol{\omega}h||{0,K} \right).$$

We can afford the above inverse estimate ∥b 1/2 ^K ^vK∥L∞(K) ^≲ ^|K^| ⁻1/d∥b 1/2 ^K vK∥0,K in the nonlinear term

$$S_6 \lesssim |K|^{-1/d} |b_K^{1/2} \boldsymbol{v}K|{0,K} |\mathsf{F}| ||\boldsymbol{u}| \, \boldsymbol{u} - |\boldsymbol{u}_h| \, \boldsymbol{u}h|{\mathbf{L}^1(K)}.$$

The elementary estimate | |a|a−|b|b | ≤ (|a|+|b|)|a−b| for vectors a, b ∈ R ^d and a Cauchy and triangle inequality provide

$$||\boldsymbol{u}|\,\boldsymbol{u}-|\boldsymbol{u}_h|\,\boldsymbol{u}h|{\mathbf{L}^1(K)} \le (|\boldsymbol{u}|_{0,K}+|\boldsymbol{u}h|{0,K})|\boldsymbol{u}-\boldsymbol{u}h|{0,K}.$$

Let us hide several constants in the notation ≲ as before but also upper bounds of ∥uh∥0,K and ∥u∥0,K to infer

$$S_6 \lesssim |K|^{-1/d} |b_K^{1/2} \boldsymbol{v}K|{0,K} |\mathsf{F}| |\boldsymbol{u} - \boldsymbol{u}h|{0,K}.$$

Let us hide the material constants |F|, κ −1 , ν, and ν −1 in the notation ≲ for a summary of the above estimates of S1,. . . , S⁶ in (6.5). After a division of the factor ∥b 1/2 ^K vK∥0,K (if positive as else there is nothing left to prove) we infer

$$|K|^{1/d} |b_K^{1/2} \boldsymbol{v}K||{0,K} \lesssim ||\boldsymbol{u} - \boldsymbol{u}h||{0,K} + ||\boldsymbol{\omega} - \boldsymbol{\omega}h||{0,K} + ||ph - p||{0,K} + \operatorname{osc}_k(\boldsymbol{f}, K).$$

The combination with (6.4) concludes the proof of the local efficiency

$$\mu_K \lesssim |\boldsymbol{u} - \boldsymbol{u}h|{0,K} + |\boldsymbol{\omega} - \boldsymbol{\omega}h|{0,K} + |ph - p|{0,K} + \operatorname{osc}_k(\boldsymbol{f}, K)$$

of the volume contribution in lower order terms. The efficiency proof does not require any smallness assumption neither on the mesh-size nor on the closeness of the exact and discrete solution. □

Algorithm 7.1 Adaptive Mesh Refinement and coarsening algorithm

1: INPUT: coarse mesh Ch, θ
                              r ∈ (0, 1), θ
                                        c ∈ (0, 1), lmax
2: OUTPUT: solution of (4.8) on an adapted mesh
3: Set Th to be the result of (optionally) applying several levels of uniform refinement to Ch
4: SOLVE the discrete problem (4.8) on Th
5: for l = 1, . . . lmax do
6: ESTIMATE: for every cell K ∈ Th, compute the local error indicator η(K) from (6.3)
7: MARK for refinement a Mr ⊂ Th with largest η(K) such that |Mr
                                                                        | ≈ θ
                                                                            r
                                                                             |Th|
8: MARK for coarsening another set Mc ⊂ Th with smallest η(K) such that |Mc
                                                                                   | ≈ θ
                                                                                        c
                                                                                        |Th|
9: ADAPT: refine K ∈ Mr and coarsen K ∈ Mc
                                                     to construct a new Th for the next step
10: SOLVE the discrete problem (4.8) on Th
11: end for

7 Numerical results

In this section we report and discuss a number of numerical examples to illustrate the performance of the proposed mixed finite element schemes and a posteriori error estimators. The realisation of the numerical methods designed in the paper is conducted with a combination of open source software packages in the so-called Gridap ecosystem [10, 49, 9]. In all cases, we use a Newton method with exact Jacobian, setting a tolerance of 10−⁸ on the ℓ ² norm of the increment and 10−¹² on the ℓ^∞ norm of the residual. The linear systems were solved either with UMFPACK (2D cases) or the multifrontal massively parallel sparse direct solver MUMPS (3D cases).

The adaptive mesh refinement procedure follows Algorithm 7.1. It comprises an standard SOLVE → ESTIMATE → MARK → ADAPT loop. This algorithm can be in principle combined with any kind of adaptive mesh approach. Our implementation of Algorithm 7.1 particularly leverages hierarchicallyadapted (i.e., nested) non-conforming octree-based meshes; see, e.g., [8]. Forest-of-octrees meshes can be seen as a two-level decomposition of Ω, referred to as macro and micro level, respectively. The macro level is a suitable conforming partition C^h of Ω into quadrilateral (d = 2) or hexahedral cells (d = 3). This mesh, which may be generated, e.g., using an unstructured mesh generator, is referred to as the coarse mesh. At the micro level, each of the cells of C^h becomes the root of an adaptive octree with cells that can be recursively and dynamically refined or coarsened using the so-called 1 : 2^d uniform partition rule. If a cell is marked for refinement, then it is split into 2^d children cells by subdividing all parent cell edges. If all children cells of a parent cell are marked for coarsening, then they are collapsed into the parent cell. The union of all leaf cells in this hierarchy forms the decomposition of the domain at the micro level, i.e., Th. While these meshes are made of quadrilaterals or cubes, we split their elements into simplices (2 triangles per mesh quadriteral in 2D, 6 tetrahedra per mesh cube in 3D) to obtain the simplicial meshes for our formulation; see, e.g., Figure 7.3.

The adaptive meshes resulting from this approach are non-conforming. In particular, they have hanging faces at the interface of neighbouring cells at different levels of refinement. Mesh nonconformity requires special adaptations to the Crouzeix–Raviart finite element space used in our formulation. In particular, following [11], we impose a set of linear multi-point constraints for the velocity degrees of freedom located at these non-conforming interfaces such that the average of the averages on the finer children faces equals the average on the coarse parent face (so-called Option C in [11]). While this approach was mathematically proven in [11] to lead to optimal convergence rates in the case of the Douglas–Santos–Sheen–Ye finite element, our numerical results confirm that this is also the case for the Crouzeix–Raviart finite element. We used the GridapP4est.jl [44] Julia package in order to handle such kind of meshes (including facet integration on non-conforming interfaces as per required by the stabilization terms in the formulation and the computation of the a posteriori error estimator) and finite element space constraints. This package, built upon the p4est meshing engine [18], is endowed with the so-called Morton space-filling curves, and it provides high-performance and low-memory footprint algorithms to handle forest-of-octrees.

| DoFs | h |-------------------------------------------------------------------------------------------------|------- | P0
P0
standard
CR
−
−
scheme, with kinematic viscosity
ν
= 1 | | 33 | 0.7071 | 145 | 0.3536 | 609 | 0.1768 | 2497 | 0.0884 | 10113 | 0.0442 | 40705 | 0.0221 | 163329 | 0.0110 | P0
P0
CR
−
−
scheme with variational crime, with kinematic viscosity
ν
= 1 | | 33 | 0.7071 | 145 | 0.3536 | 609 | 0.1768 | 2497 | 0.0884 | 10113 | 0.0442 | 40705 | 0.0221 | 163329 | 0.0110 | | | 33 | 0.7071 | 145 | 0.3536 | 609 | 0.1768 | 2497 | 0.0884 | 10113 | 0.0442 | 40705 | 0.0221 | 163329 | 0.0110 | | 33 | 0.7071 | 145 | 0.3536 | 609 | 0.1768 | 2497 | 0.0884 | 10113 | 0.0442 | 40705 | 0.0221 | 163329 | 0.0110 | ∥u
−
uh∥h | rate | ∥ω
−
ωh∥0,Ω | rate | ∥p
−
ph∥0,Ω | rate | lossdiv | losscurl | | --------|-----------------|--------------|---------------------------------------------------------|--------|----------------------------------|-------|-------------|----------|--| | | | | | | | | | | | 6.09e-02 | ⋆ | 5.60e-02 | ⋆ | 1.93e-01 | ⋆ | 3.47e-18 | 1.04e-17 | | | 3.56e-02 | 0.774 | 3.26e-02 | 0.779 | 1.01e-01 | 0.925 | 6.94e-18 | 2.78e-17 | | | 1.82e-02 | 0.969 | 1.64e-02 | 0.994 | 5.27e-02 | 0.945 | 2.43e-17 | 2.78e-17 | | | 9.05e-03 | 1.007 | 8.04e-03 | 1.028 | 2.72e-02 | 0.952 | 1.38e-16 | 6.94e-17 | | | 4.50e-03 | 1.008 | 3.97e-03 | 1.019 | 1.38e-02 | 0.975 | 1.08e-15 | 2.06e-16 | | | 2.25e-03 | 1.003 | 1.97e-03 | 1.007 | 6.97e-03 | 0.989 | 1.80e-15 | 4.02e-16 | | | 1.12e-03 | 1.001 | 9.86e-04 | 1.002 | 3.50e-03 | 0.995 | 1.49e-13 | 8.33e-16 | | | | | | | | | | | | | 5.59e-02 | ⋆ | 5.38e-02 | ⋆ | 1.71e-01 | ⋆ | 1.73e-18 | 4.34e-18 | | | 3.43e-02 | 0.705 | 3.30e-02 | 0.705 | 9.31e-02 | 0.873 | 6.94e-18 | 2.78e-17 | | | 1.75e-02 | 0.971 | 1.66e-02 | 0.995 | 4.91e-02 | 0.923 | 2.78e-17 | 4.16e-17 | | | 8.66e-03 | 1.015 | 8.08e-03 | 1.036 | 2.57e-02 | 0.932 | 6.77e-17 | 8.33e-17 | | | 4.30e-03 | 1.011 | 3.97e-03 | 1.023 | 1.32e-02 | 0.964 | 1.19e-15 | 2.26e-16 | | | 2.14e-03 | 1.003 | 1.98e-03 | 1.008 | 6.67e-03 | 0.984 | 4.77e-15 | 4.16e-16 | | | 1.07e-03 | 1.001 | 9.86e-04 | 1.002 | 3.35e-03 | 0.993 | 5.55e-15 | 1.05e-15 | | | standard
CR | P0
−
− | P0 | | scheme, with kinematic viscosity | ν | = 10−4 | | | | 3.85e-02 | ⋆ | 1.31e-03 | ⋆ | 1.81e-01 | ⋆ | 1.39e-17 | 8.67e-19 | | | 3.78e-02 | 0.029 | 4.25e-03 | -1.699 | 9.30e-02 | 0.958 | 2.78e-17 | 1.73e-18 | | | 1.33e-02 | 1.503 | 2.43e-03 | 0.806 | 4.72e-02 | 0.977 | 8.33e-17 | 1.73e-18 | | | 5.68e-03 | 1.230 | 1.37e-03 | 0.823 | 2.37e-02 | 0.994 | 8.33e-17 | 1.30e-18 | | | 2.56e-03 | 1.151 | 6.28e-04 | 1.129 | 1.19e-02 | 0.999 | 1.18e-16 | 1.30e-18 | | | 1.21e-03 | 1.083 | 2.00e-04 | 1.653 | 5.94e-03 | 1.000 | 4.02e-16 | 3.79e-18 | | | 5.89e-04 | 1.034 | 5.50e-05 | 1.860 | 2.97e-03 | 1.000 | 9.99e-16 | 8.24e-18 | | | P0
CR
− | P0
− | | scheme with variational crime, with kinematic viscosity | | | | = 10−4
ν | | | | 5.79e-03 | ⋆ | 5.89e-04 | ⋆ | 1.70e-01 | ⋆ | 1.73e-18 | 2.71e-20 | | | 4.44e-03 | 0.384 | 4.31e-04 | 0.449 | 9.26e-02 | 0.880 | 6.94e-18 | 1.08e-19 | | | 1.20e-03 | 1.883 | 1.86e-04 | 1.216 | 4.72e-02 | 0.972 | 2.26e-17 | 4.34e-19 | | | 3.48e-04 | 1.789 | 8.82e-05 | 1.073 | 2.37e-02 | 0.993 | 6.42e-17 | 8.67e-19 | | | 1.03e-04 | 1.758 | 4.31e-05 | 1.033 | 1.19e-02 | 0.998 | 1.11e-16 | 1.36e-18 | | | 3.39e-05 | 1.602 | 2.12e-05 | 1.023 | 5.94e-03 | 1.000 | 2.22e-16 | 2.87e-18 | | | 1.35e-05 | 1.329 | 1.05e-05 | 1.016 | 2.97e-03 | 1.000 | 4.88e-15 | 9.81e-18 | |

Table 7.1: Error decay with respect to smooth manufactured solutions in 2D for two methods, and a large and a small viscosity. This illustrates the (Bernoulli) pressure robustness property, evidenced more clearly in the third and fourth blocks of the table for (4.8).

Verification of convergence to smooth solutions. First we exemplify the convergence of the method against smooth manufactured solutions and illustrate the pressure robustness property. For this we proceed similarly as in [40, Sect. 5.2] and consider the domain Ω = (0, 1)² with two values for the kinematic viscosity ν = 1 and ν = 10−⁴ , and with two discretisations (including or not the interpolation of suitable test functions in ah, F^h and c u^h h that lead to a variational crime). We also consider a manufactured stream function ξ = x 2 (1−x) 2y 2 (1−y 2 ) and a manufactured velocity, vorticity and Bernoulli pressure that solve the coupled system (2.2) and that satisfy also the homogeneous boundary condition and zero-mean Bernoulli pressure constraint (2.3), as follows

$$\boldsymbol{u} = \operatorname{\mathbf{curl}} \xi, \quad \boldsymbol{\omega} = \sqrt{\nu} \operatorname{\mathbf{curl}} \boldsymbol{u}, \quad p = x^3 + y^3 - \frac{1}{2}.$$

The remaining model parameters and stabilisation constant are taken as F = 1, κ = 1, and ϑ = 10, respectively.

We generate successively refined simplicial grids and compute errors between the approximate and exact solutions on each refinement level. For uniform mesh refinement, the experimental convergence

Figure 7.1: Approximate velocity (line integral contours and magnitude), vorticity, and Bernoulli pressure profiles computed with the modified CR−P0−P⁰ scheme, with kinematic viscosity ν = 10−⁴ .

order is computed as

$$\mathtt{rate} = \frac{\log(ei(\bullet)) - \log(e{i+1}(\bullet))}{\log(hi) - \log(h{i+1})}$$

with ei(•) denoting the error associated with the quantity • in its natural norm and hⁱ the mesh size corresponding to a refinement level i. The error history (errors and experimental convergence rates) shown in Table 7.1 confirms the optimal convergence of the non-conforming scheme with and without the variational crime approach, for all variables in their respective norms. In particular we confirm that for the larger viscosity the methods deliver comparable results. On the other hand, for the smaller viscosity, with the modified scheme the velocity error (measured in the broken norm) is almost two orders of magnitude smaller than that obtained with the standard scheme (the latter violates the invariance condition discussed in [43] since the applied gradient force on the right-hand side induces an incompatible velocity field). In the table we also report on the local kernel characterising properties divhu^h = 0 and curl ^hu^h − ω^h = 0, shown by projecting onto the pressure and vorticity spaces these residuals, and taking the ℓ^∞ norm of the corresponding vector of degrees of freedom. We obtain a machine precision accuracy for all cases. The approximate solutions computed with the modified scheme using the second last mesh refinement and the smaller viscosity are shown in Figure 7.1, indicating well-resolved patterns.

We note that for this relatively small problem we simply use a direct method for each tangent system in the linearisation process, and confirm that, in order to reach the stopping criterion on each refinement level, at most four and two Newton–Raphson iterations were needed for the small and large viscosity cases, respectively.

In order to illustrate the need for the stabilisation parameter θ as discussed in Section 4, we run again – only for four coarse mesh refinement levels – the test reported in the fourth block of Table 7.1 (CR − P⁰ − P⁰ modified scheme, with kinematic viscosity ν = 10−⁴ ) and compare that block with the results in Table 7.2, that use ϑ ∈ {0, 0.01, 1}. Without stabilisation the method does not converge in the velocity nor vorticity, whereas for not big enough stabilisation one observes a slightly suboptimal convergence in the velocity field (as well as a larger Newton iteration count). Again as in Table 7.1, the energy error in the modified method is not affected by the relatively larger Bernoulli pressure error.

We also showcase the convergence properties in 3D by considering the following smooth manufactured solutions in Ω = (0, 1)³

$$\mathbf{u} = \begin{pmatrix} \sin(\pi x)\cos(\pi y)\cos(\pi z) \ -2\cos(\pi x)\sin(\pi y)\cos(\pi z) \ \cos(\pi x)\cos(\pi y)\sin(\pi z) \end{pmatrix}, \quad \boldsymbol{\omega} = \sqrt{\nu} \operatorname{\mathbf{curl}} \mathbf{u}, \quad p = \sin(\pi x)\sin(\pi y)\sin(\pi z)$$

| DoFs | h | ∥u
−
uh∥h | rate | ∥ω
−
ωh∥0,Ω | rate | ∥p
−
ph∥0,Ω | rate | lossdiv | losscurl | | | |----------|-------------|-----------------|--------|-------------------|--------|-------------------|-------|----------|----------|--|--| | ϑ
= 0 | | | | | | | | | | | | | 33 | 0.7071 | 1.94e-02 | ⋆ | 5.67e-04 | ⋆ | 1.70e-01 | ⋆ | 6.94e-18 | 4.82e-20 | | | | 145 | 0.3536 | 5.81e-02 | -1.580 | 5.72e-04 | -0.012 | 9.26e-02 | 0.880 | 1.39e-17 | 2.22e-19 | | | | 609 | 0.1768 | 1.40e-01 | -1.273 | 5.72e-04 | 0.000 | 4.72e-02 | 0.972 | 5.55e-17 | 5.52e-19 | | | | 2497 | 0.0884 | 2.96e-01 | -1.077 | 5.72e-04 | 0.000 | 2.37e-02 | 0.993 | 2.22e-16 | 1.29e-18 | | | | 10113 | 0.0442 | 6.01e-01 | -1.020 | 5.71e-04 | 0.001 | 1.19e-02 | 0.998 | 4.44e-16 | 3.38e-18 | | | | | ϑ
= 0.01 | | | | | | | | | | | | 33 | 0.7071 | 1.76e-02 | ⋆ | 5.84e-04 | ⋆ | 1.70e-01 | ⋆ | 2.17e-19 | 2.37e-20 | | | | 145 | 0.3536 | 3.68e-02 | -1.067 | 5.61e-04 | 0.057 | 9.26e-02 | 0.880 | 2.08e-17 | 1.08e-19 | | | | 609 | 0.1768 | 3.56e-02 | 0.048 | 3.88e-04 | 0.534 | 4.72e-02 | 0.972 | 2.78e-17 | 3.79e-19 | | | | 2497 | 0.0884 | 2.21e-02 | 0.690 | 2.32e-04 | 0.742 | 2.37e-02 | 0.993 | 5.55e-17 | 5.42e-19 | | | | 10113 | 0.0442 | 1.18e-02 | 0.904 | 1.24e-04 | 0.902 | 1.19e-02 | 0.998 | 1.11e-16 | 1.63e-18 | | | | | ϑ
= 1 | | | | | | | | | | | | 33 | 0.7071 | 5.93e-03 | ⋆ | 6.46e-04 | ⋆ | 1.70e-01 | ⋆ | 6.94e-18 | 1.08e-19 | | | | 145 | 0.3536 | 3.97e-03 | 0.579 | 3.76e-04 | 0.780 | 9.26e-02 | 0.880 | 8.67e-18 | 2.17e-19 | | | | 609 | 0.1768 | 9.68e-04 | 2.037 | 1.70e-04 | 1.148 | 4.72e-02 | 0.972 | 2.78e-17 | 4.34e-19 | | | | 2497 | 0.0884 | 3.29e-04 | 1.555 | 8.42e-05 | 1.010 | 2.37e-02 | 0.993 | 5.55e-17 | 5.42e-19 | | | | 10113 | 0.0442 | 1.41e-04 | 1.222 | 4.20e-05 | 1.004 | 1.19e-02 | 0.998 | 1.73e-16 | 1.68e-18 | | |

Table 7.2: Error decay with respect to mesh refinement (only coarser levels are shown) using smooth manufactured solutions in 2D for the modified CR − P⁰ − P⁰ scheme, with small kinematic viscosity ν = 10−⁴ and three values of the penalisation parameter (compare also with the fourth block in Table 7.1).

DoFs	h	∥u − uh∥h	rate	∥ω − ωh∥0,Ω	rate	∥p − ph∥0,Ω	rate	lossdiv	losscurl
43	0.6124	1.94e+00	⋆	3.42e-01	⋆	3.15e-01	⋆	4.44e-16	5.55e-17
409	0.3062	1.33e+00	0.417	2.59e-01	0.402	2.24e-01	0.495	1.78e-15	2.22e-16
3553	0.1531	9.83e-01	0.432	1.64e-01	0.663	1.13e-01	0.987	4.44e-15	5.83e-16
29569	0.0765	5.16e-01	0.931	8.68e-02	0.915	5.21e-02	1.114	1.34e-14	1.78e-15
241153	0.0383	2.57e-01	1.004	4.39e-02	0.982	2.52e-02	1.046	3.60e-14	4.22e-15
1947649	0.0191	1.28e-01	1.002	2.21e-02	0.994	1.25e-02	1.015	7.90e-14	8.88e-15

Table 7.3: Error decay with respect to mesh refinement using smooth manufactured solutions in 3D for the modified (with variational crime) CR−P⁰ −P⁰ scheme, with ν = 0.01, κ = 100, F = 10, ϑ = 1.

and take model and stabilisation parameters as follows ν = 0.01, κ = 100, F = 10, ϑ = 1. The error history is presented in Table 7.3, where the results are consistent with the 2D case: they confirm the optimal order of convergence in the three unknowns, and give evidence of the expected conservation properties. For these runs the Newton–Raphson algorithm has taken up to six iterations to converge on each mesh refinement. The components of the approximate solution are portrayed in Figure 7.2.

Testing the robustness of the posteriori error estimator with non-smooth solutions. Next, we assess the convergence of the method when approximating a non-smooth exact solution under uniform mesh refinement, and check how the adaptive mesh refinement guided by the a posteriori error estimator defined in (6.3) is able to restore optimal convergence. For this we take unity parameters ν = F = κ = 1, and use the following manufactured solutions defined on the L-shaped domain Ω = (−1, 1)² \ ([0, 1) × (−1, 0]) (see, e.g., [21, 50]):

$$\mathbf{u} = r^{\lambda} \begin{pmatrix} (1+\lambda)\sin(\theta)\psi(\theta) + \cos(\theta)\psi'(\theta) \ \sin(\theta)\psi'(\theta) - (1+\lambda)\cos(\theta)\psi(\theta) \end{pmatrix}, \quad \boldsymbol{\omega} = \sqrt{\nu}\mathbf{curl}\,\mathbf{u},$$ $$p = -\nu \frac{r^{\lambda-1}}{1-\lambda} ((1+\lambda)^2 \psi'(\theta) + \psi'''(\theta))$$

Figure 7.2: Approximate velocity (streamlines), vorticity (streamlines), and Bernoulli pressure profiles on a coarse mesh, computed with the modified CR − P⁰ − P⁰ scheme, with kinematic viscosity ν = 1.

in polar coordinates centred at the origin (r, θ) ∈ (0, ∞) × (0, 3π 2 ), where

$$\psi(\theta) = \frac{\sin((1+\lambda)\theta)\cos(\lambda w)}{1+\lambda} - \cos((1+\lambda)\theta) - \frac{\sin((1-\lambda)\theta)\cos(\lambda w)}{1-\lambda} + \cos((1-\lambda)\theta).$$

Here λ = 856399 ^1572864 ≈ 0.5444837 is the smallest positive solution of sin(λw) + λ sin(w) = 0, and we take w = 2 . Second-order derivatives for velocity (and first order derivatives for pressure and vorticity) are not square integrable, and therefore these solutions do not have higher regularity (u ∈/ H¹ (Ω), ω, p /∈ H¹ (Ω)). Nevertheless, the exact boundary velocity is zero on the reentrant edges (at θ = 0 and θ = 2 ) and so the boundary data oscillation can be considered of high order. Note that the exact velocity and pressure above are such that ^√ νcurl ω + ∇p = 0 and so f = κ ⁻1u + ν ⁻1/2ω × u + F|u|u. For this test the stabilisation parameter is taken as ϑ = 10.

We show in Figure 7.3 the approximate solutions as well as a sample of adaptively refined meshes generated by Algorithm 7.1 (where it is evident that the solution singularity induces a concentration of elements near the reentrant corner). In Table 7.4 we display the difference in error decay between the uniform and adaptive mesh refinement as well as the effectivity index

$$\text{eff}(\eta) = \frac{|\boldsymbol{u} - \boldsymbol{u}_h|_h + |\boldsymbol{\omega} - \boldsymbol{\omega}h|{0,\Omega} + |p - ph|{0,\Omega}}{\eta(\mathcal{T}_h)}$$

in both cases. For adaptive mesh refinement, the experimental convergence order is computed as

$$\mathtt{rate} = \frac{\log(ei(\bullet)) - \log(e{i+1}(\bullet))}{-\frac{1}{2}[\log(\mathtt{DoF}i)/\log(\mathtt{DoF}{i+1})]}$$

with ei(•) the error associated with the quantity • in its natural norm and DoFⁱ the total number of degrees of freedom corresponding to a refinement level i. We used a coarse mesh C^h with three quadtrees to mesh the L-shaped domain, no initial uniform refinements, and θ ^r = 0.275 and θ ^c = 0.0 in Algorithm 7.1. We can observe that the error under uniform mesh refinement goes to zero with a suboptimal rate (the expected O(h λ )) while for the adaptive case we recover optimal linear convergence, necessitating much fewer degrees of freedom to achieve the same error. In addition, from the last column of the table we can see that the effectivity index remains bounded between 1.6 and 1.9 for all cases. Irrespective of the refinement strategy and refinement level, the Newton–Raphson scheme has taken no more than four iterations to reach the prescribed tolerance of 10−⁸

8 Concluding remarks

We have developed and analysed a novel non-augmented vorticity–velocity–Bernoulli formulation for incompressible flow in highly permeable porous media, governed by the Navier–Stokes–Brinkman– Forchheimer equations. The approach avoids least-squares or stabilisation-based augmentation by

Figure 7.3: Sample of approximate solutions for the convergence test on an L-shaped domain (top rows), and coarse meshes produced after three, six, and nine steps of the adaptive refinement algorithm guided by the a posteriori error estimator.

| DoFs | h | ∥u
−
uh∥h | rate | ∥ω
−
ωh∥0,Ω | rate | ∥p
−
ph∥0,Ω | rate | lossdiv | losscurl | eff(η) | | |--------|-------------------------|-----------------|-------|--------------------------|-------|-------------------|-------|----------|----------|--------|--| | | uniform mesh refinement | | | | | | | | | | | | 23 | 0.7071 | 2.64e+0 | ⋆ | 2.50e+0 | ⋆ | 2.57e+0 | ⋆ | 6.21e-11 | 8.88e-16 | 1.635 | | | 105 | 0.3536 | 1.91e+0 | 0.467 | 1.80e+0 | 0.473 | 2.00e+0 | 0.366 | 2.12e-11 | 1.78e-15 | 1.826 | | | 449 | 0.1768 | 1.36e+0 | 0.489 | 1.27e+0 | 0.496 | 1.69e+0 | 0.241 | 7.28e-12 | 4.44e-15 | 1.810 | | | 1857 | 0.0884 | 9.51e-01 | 0.516 | 8.91e-01 | 0.517 | 1.23e+0 | 0.459 | 2.50e-12 | 8.88e-15 | 1.797 | | | 7553 | 0.0442 | 6.59e-01 | 0.529 | 6.17e-01 | 0.530 | 8.61e-01 | 0.514 | 8.81e-13 | 3.73e-14 | 1.790 | | | 30465 | 0.0221 | 4.54e-01 | 0.537 | 4.25e-01 | 0.537 | 5.96e-01 | 0.530 | 3.41e-13 | 8.26e-14 | 1.786 | | | 122369 | 0.0110 | 3.12e-01 | 0.540 | 2.92e-01 | 0.541 | 4.11e-01 | 0.537 | 2.27e-13 | 1.55e-13 | 1.784 | | | | | | | adaptive mesh refinement | | | | | | | | | 23 | 0.7071 | 2.64e+0 | ⋆ | 2.50e+0 | ⋆ | 2.57e+0 | ⋆ | 6.21e-11 | 8.88e-16 | 1.635 | | | 79 | 0.3536 | 2.18e+0 | 0.305 | 2.05e+0 | 0.317 | 2.14e+0 | 0.155 | 1.27e-06 | 8.88e-16 | 1.876 | | | 169 | 0.1768 | 1.64e+0 | 0.757 | 1.53e+0 | 0.768 | 1.86e+0 | 0.609 | 2.12e-11 | 3.55e-15 | 1.811 | | | 405 | 0.0884 | 1.25e+0 | 0.625 | 1.16e+0 | 0.646 | 1.62e+0 | 0.212 | 7.42e-09 | 3.55e-15 | 1.782 | | | 921 | 0.0442 | 8.94e-01 | 0.808 | 8.26e-01 | 0.817 | 1.31e+0 | 0.572 | 2.50e-12 | 7.11e-15 | 1.759 | | | 1891 | 0.0221 | 6.54e-01 | 0.871 | 6.00e-01 | 0.886 | 1.03e+0 | 0.726 | 8.63e-13 | 1.24e-14 | 1.726 | | | 3833 | 0.0110 | 4.72e-01 | 0.920 | 4.33e-01 | 0.925 | 7.72e-01 | 0.815 | 3.13e-13 | 1.78e-14 | 1.698 | | | 7771 | 0.0055 | 3.37e-01 | 0.960 | 3.07e-01 | 0.973 | 5.72e-01 | 0.848 | 1.21e-13 | 2.84e-14 | 1.671 | | | 15153 | 0.0028 | 2.39e-01 | 1.020 | 2.18e-01 | 1.032 | 4.22e-01 | 0.911 | 6.39e-14 | 4.26e-14 | 1.649 | | | 29121 | 0.0014 | 1.72e-01 | 1.011 | 1.56e-01 | 1.024 | 3.13e-01 | 0.919 | 5.68e-14 | 7.82e-14 | 1.625 | |

Table 7.4: Error decay, divergence and curl loss, and effectivity index tabulated with respect to mesh refinement using non-smooth manufactured solutions in the L-shaped domain for the modified CR − P⁰ − P⁰ scheme, and using uniform vs. adaptive mesh refinement.

formulating the problem as a nested saddle-point system and working in divergence-free velocity spaces. This framework enables both continuous and discrete solvability analyses under small data assumptions, relying on fixed-point arguments and inf-sup conditions in divergence-free subspaces. At the discrete level, we employ a Crouzeix–Raviart finite element method enhanced with tangential and normal jump penalisation to control the full H(curl) ∩ H(div) velocity norm, and we ensure pressure-robustness through a modified variational form using a Raviart–Thomas interpolant. We introduced a fully computable residual-based a posteriori error estimator, which is shown to be reliable and efficient—the latter hinging on a new inverse inequality tailored to the Forchheimer nonlinearity. Numerical experiments validate all theoretical results and demonstrate that adaptive mesh refinement, driven by the proposed estimator, restores optimal convergence even in the presence of solution singularities. The practical implementation leverages a lightweight, parallel octree-based infrastructure capable of generating non-conforming adaptive meshes, where appropriate multipoint constraints allow the recovery of optimal approximation properties. The combined contributions in formulation, analysis, estimator design, and implementation provide a foundation for further development of adaptive solvers for nonlinear porous flow models and illustrate the benefits of structure-preserving discretisations in complex multiphysics regimes.

Possible directions of future work include the extension of the analysis to more general data regimes, particularly beyond the smallness assumptions required for fixed-point arguments, potentially through monotonicity or compactness-based techniques. It would also be valuable to investigate higher-order nonconforming discretisations that retain divergence-free properties and support pressure-robustness. Finally, we want to address the coupling with transport and reaction models within the same porous medium, where interface dynamics or multiphysics constraints may necessitate tailored trace operators and stabilisation.

References

[1] N. Ahmed, G. R. Barrenechea, E. Burman, J. Guzman, A. Linke, and C. Merdon ´ , A pressure-robust discretization of Oseen's equation using stabilization in the vorticity equation, SIAM J. Numer. Anal., 59 (2021), pp. 2746–2774.
[2] M. Amara, D. Capatina-Papaghiuc, and D. Trujillo, Stabilized finite element method for Navier-Stokes equations with physical boundary conditions, Math. Comp., 76 (2007), pp. 1195– 1217.
[3] V. Anaya, A. Bouharguane, D. Mora, C. Reales, R. Ruiz-Baier, N. Seloula, and H. Torres, Analysis and approximation of a vorticity–velocity–pressure formulation for the Oseen equations, J. Sci. Comput., 80 (2019), pp. 1577–1606.
[4] V. Anaya, R. Caraballo, S. Caucao, L. F. Gatica, R. Ruiz-Baier, and I. Yotov, A vorticity-based mixed formulation for the unsteady Brinkman–Forchheimer equations, Comput. Methods Appl. Mech. Engrg., 404 (2023), p. 115829.
[5] V. Anaya, R. Caraballo, B. Gomez-Vargas, D. Mora, and R. Ruiz-Baier, Velocityvorticity-pressure formulation for the Oseen problem with variable viscosity, Calcolo, 58 (2021), pp. 1–25.
[6] V. Anaya, R. Caraballo, R. Ruiz-Baier, and H. Torres, On augmented finite element formulations for the Navier–Stokes equations with vorticity and variable viscosity, Comput. Math. Appl., 143 (2023), pp. 397–416.
[7] V. Anaya, D. Mora, R. Oyarzua, and R. Ruiz-Baier ´ , A priori and a posteriori error analysis of a mixed scheme for the Brinkman problem, Numer. Math., 133 (2016), pp. 781–817.
[8] S. Badia, A. F. Mart´ın, E. Neiva, and F. Verdugo, A generic finite element framework on parallel tree-based adaptive meshes, SIAM J. Sci. Comput., 42 (2020), pp. C436–C468.
[9] S. Badia, A. F. Mart´ın, and F. Verdugo, GridapDistributed: a massively parallel finite element toolbox in Julia, J. Open Source Softw., 7 (2022), p. 4157.
[10] S. Badia and F. Verdugo, Gridap: An extensible finite element toolbox in julia, Journal of Open Source Software, 5 (2020), p. 2520.
[11] W. Bangerth, I. Kim, D. Sheen, and J. Yim, On hanging node constraints for nonconforming finite elements using the Douglas–Santos–Sheen–Ye element as an example, SIAM J. Numer. Anal., 55 (2017), pp. 1719–1739.
[12] M. Barker, S. Cao, and A. Stern, A nonconforming primal hybrid finite element method for the two-dimensional vector Laplacian, SMAI J. Comput. Math., 10 (2024), pp. 85–106.
[13] L. Beirao da Veiga, F. Dassi, and G. Vacca ˜ , Vorticity-stabilized virtual elements for the Oseen equation, Math. Models Methods Appl. Sci., 31 (2021), pp. 3009–3052.
[14] M. Benzi, M. A. Olshanskii, L. G. Rebholz, and Z. Wang, Assessment of a vorticity based solver for the Navier–Stokes equations, Comput. Methods Appl. Mech. Engrg., 247-248 (2012), pp. 216–225.
[15] C. Brennecke, A. Linke, C. Merdon, and J. Schoberl ¨ , Optimal and pressure-independent L 2 velocity error estimates for a modified Crouzeix-Raviart Stokes element with BDM reconstructions, J. Comput. Math., 33 (2015), pp. 191–208.
[16] S. C. Brenner, J. Cui, F. Li, and L.-Y. Sung, A nonconforming finite element method for a two-dimensional curl–curl and grad-div problem, Numer. Math., 109 (2008), pp. 509–533.
[17] S. C. Brenner, F. Li, and L.-Y. Sung, A locally divergence-free nonconforming finite element method for the time-harmonic Maxwell equations, Math. Comp., 76 (2007), pp. 573–595.
[18] C. Burstedde, L. C. Wilcox, and O. Ghattas, p4est: Scalable Algorithms for Parallel Adaptive Mesh Refinement on Forests of Octrees, SIAM J. Sci. Comput., 33 (2011), pp. 1103– 1133.
[19] J. Camano, C. Garc ˜ ´ıa, and R. Oyarzua´ , Analysis of a momentum conservative mixed-FEM for the stationary Navier–Stokes problem, Numer. Methods Partial Differential Equations, 37 (2021), pp. 2895–2923.
[20] R. Caraballo, C. W. In, A. F. Mart´ın, and R. Ruiz-Baier, Robust finite element methods and solvers for the Biot–Brinkman equations in vorticity form, Numer. Methods Partial Differential Equations, 40 (2024), p. e23083.
[21] C. Carstensen, P. Causin, and R. Sacco, A posteriori dual-mixed adaptive finite element error control for Lam´e and Stokes equations, Numer. Math., 101 (2005), pp. 309–332.
[22] C. Carstensen, D. Gallistl, and M. Schedensack, Adaptive nonconforming Crouzeix-Raviart FEM for eigenvalue problems, Math. Comp., 84 (2015), pp. 1061–1087.
[23] C. Carstensen, J. Gedicke, and D. Rim, Explicit error estimates for Courant, Crouzeix-Raviart and Raviart-Thomas finite element methods, J. Comput. Math., 30 (2012), pp. 337–353.
[24] C. Carstensen, B. Graßle, and N. Nataraj ¨ , Unifying a posteriori error analysis of five piecewise quadratic discretisations for the biharmonic equation, J. Numer. Math., 32 (2024), pp. 77–109.
[25] C. Carstensen and C. Merdon, Computational survey on a posteriori error estimators for the Crouzeix–Raviart nonconforming finite element method for the Stokes problem, Comput. Methods Appl. Math., 14 (2014), pp. 35–54.
[26] C. Carstensen and S. Puttkammer, Adaptive guaranteed lower eigenvalue bounds with optimal convergence rates, Numer. Math., 156 (2024), pp. 1–38.
[27] S. Caucao, R. Oyarzua, and S. Villa-Fuentes ´ , A new mixed-FEM for steady-state natural convection models allowing conservation of momentum and thermal energy, Calcolo, 57 (2020), p. 36.
[28] M. Crouzeix and P.-A. Raviart, Conforming and nonconforming finite element methods for solving the stationary Stokes equations I, RAIRO: Math., 7 (1973), pp. 33–75.
[29] A. Ern and J.-L. Guermond, Theory and Practice of Finite Elements, Applied Mathematical Sciences, 159. Springer-Verlag, New York, 2004.
[30] G. N. Gatica, A Simple Introduction to the Mixed Finite Element Method, Springer-Verlag, Berlin, 2014.
[31] V. Girault and P.-A. Raviart, Finite element approximation of the Navier-Stokes equations, vol. 749 of Lecture Notes in Mathematics, Springer-Verlag, Berlin-New York, 1979.
[32] M.-L. Hanot, An arbitrary order and pointwise divergence-free finite element scheme for the incompressible 3D Navier–Stokes equations, SIAM J. Numer. Anal., 61 (2023), pp. 784–811.
[33] P. Hansbo and M. G. Larson, Discontinuous Galerkin and the Crouzeix–Raviart element: application to elasticity, ESAIM: Math. Model. Numer. Anal., 37 (2003), pp. 63–72.
[34] P. Hansbo and T. Rylander, A linear nonconforming finite element method for Maxwell's equations in two dimensions. Part I: Frequency domain, J. Comput. Phys., 229 (2010), pp. 6534– 6547.
[35] V. John, Finite element methods for incompressible flow problems, vol. 51, Springer Verlag, Berlin-New York, 2016.
[36] P. Knobloch, On Korn's inequality for nonconforming finite elements, Tech. Mech., 20 (2000), pp. 205–214.
[37] P. Knobloch and L. Tobiska, On Korn's first inequality for quadrilateral nonconformimg finite elements of first order approximation properties, Int. J. Numer. Anal. Model., 2 (2005), pp. 439–458.
[38] W. Layton, C. C. Manica, M. Neda, M. Olshanskii, and L. G. Rebholz, On the accuracy of the rotation form in simulations of the Navier–Stokes equations, J. Comput. Phys., 228 (2009), pp. 3433–3447.
[39] P. L. Lederer, C. Merdon, and J. Schoberl ¨ , Refined a posteriori error estimation for classical and pressure-robust Stokes finite element methods, Numer. Math., 142 (2019), pp. 713– 748.
[40] A. 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), pp. 782–800.
[41] A. Linke and C. Merdon, Guaranteed energy error estimators for a modified robust Crouzeix– Raviart Stokes element, J. Sci. Comput., 64 (2015), pp. 541–558.
[42] , Pressure-robustness and discrete Helmholtz projectors in mixed finite element methods for the incompressible Navier–Stokes equations, Comput. Methods Appl. Mech. Engrg., 311 (2016), pp. 304–326.
[43] A. Linke, C. Merdon, and W. Wollner, Optimal L2 velocity error estimate for a modified pressure-robust Crouzeix–Raviart Stokes element, IMA J. Numer. Anal., 37 (2017), pp. 354–374.
[44] A. F. Mart´ın, Gridapp4est.jl, 2025.
[45] A. Quarteroni and A. Valli, Numerical Approximation of Partial Differential Equations, vol. 23 of Springer Series in Computational Mathematics, Springer-Verlag, Berlin, 1994.
[46] M. Salaun and S. Salmon ¨ , Low-order finite element method for the well-posed bidimensional Stokes problem, IMA J. Numer. Anal., 35 (2015), pp. 427–453.
[47] C.-C. Tsai and S.-Y. Yang, On the velocity–vorticity–pressure least-squares finite element method for the stationary incompressible Oseen problem, J. Comput. Appl. Math., 182 (2005), pp. 211–232.
[48] S. Turek and A. Ouazzi, Unified edge-oriented stabilization of nonconforming FEM for incompressible flow problems: Numerical investigations, J. Numer. Math., 15 (2007).
[49] F. Verdugo and S. Badia, The software design of Gridap: A Finite Element package based on the Julia JIT compiler, Computer Physics Communications, 276 (2022), p. 108341.
[50] R. Verfurth ¨ , A posteriori error estimators for the Stokes equations, Numer. Math., 55 (1989), pp. 309–325.
[51] R. Verfurth and P. Zanotti ¨ , A quasi-optimal Crouzeix–Raviart discretization of the Stokes equations, SIAM J. Numer. Anal., 57 (2019), pp. 1082–1099.
[52] Y. Zhang, A. Palha, M. Gerritsma, and Q. Yao, A MEEVC discretization for twodimensional incompressible Navier-Stokes equations with general boundary conditions, J. Comput. Phys., 510 (2024), p. 113080.

A A new inverse estimate

Let us first recall that the scalar volume-bubble function b^K ∈ W 1,∞ 0 (K) on a simplex K ⊂ R d is the product of all d + 1 barycentric coordinates of the vertices of K times a factor (d + 1)d+1 for the normalisation 0 ≤ b^K ≤ 1 = max b^K in K. Note that typical efficiency estimates rely upon Verf¨urth's bubble-function methodology [50]. For the present case, the Forchheimer nonlinearity requires a nontrivial modification of this approach, presented in the following key inverse estimate.

Lemma A.1 (New inverse estimate) There exists a universal positive constant Ceq, that exclusively depends on d = 2, 3, such that all affine vector-valued functions f, g ∈ P1(K) d satisfy

$$|||f|f+g||{0,K} \le C{\text{eq}} ||bK^{1/2}(|f|f+g)||{0,K}.$$

Proof of independence of K. At first glance, the most striking aspect of the lemma is the independence of K, that makes Ceq a constant that depends on the reference simplex Kref of volume |Kref| > 0. But that aspect follows from an elementary affine transformation Ψ from K to Kref: All integrals of both sides of the squared estimate are Lebesgue functions (and no derivative appears) so the transformation shows

$$||\boldsymbol{f}|\boldsymbol{f}+\boldsymbol{g}|{0,K} = \sqrt{\frac{|K|}{|K{\text{ref}}|}} \, ||\hat{\boldsymbol{f}}|\hat{\boldsymbol{f}}+\hat{\boldsymbol{g}}|{0,K{\text{ref}}} \quad \text{and} \quad |bK^{1/2}(|\boldsymbol{f}|\boldsymbol{f}+\boldsymbol{g})|{0,K} = \sqrt{\frac{|K|}{|K{\text{ref}}|}} \, |b{K{\text{ref}}}^{1/2}(|\hat{\boldsymbol{f}}|\hat{\boldsymbol{f}}+\hat{\boldsymbol{g}})|{0,K_{\text{ref}}}$$

for ^f^ˆ := ^f ◦ Ψ, ^g^ˆ := ^g ◦ Ψ, ^bKref ⁼ ^b^K ◦ Ψ, and ^Kref = Ψ(K). Thus, once the assertion holds on the reference simplex K = Kref, it holds on any simplex K with the same constant Ceq. □

Proof in case g = 0 or f ∈ P0(K) d . If g = 0 then only the term |f|f with Euclid length |f| ² arises and |f| 2 is a quadratic polynomial (for the affine f). Therefore a classical inverse estimate provides Ceq,2 with

$$||q2||{0,K} \le C_{\text{eq},2} ||b_K^{1/2} q2||{0,K} \quad \text{for all } q_2 \in \mathbb{P}_2(K)^d.$$

This concludes the proof of the assertion for f ∈ P1(K) d , g = 0, and Ceq = Ceq,2:

$$|||f|f||{0,K} = |||f|^2||{0,K} \le C_{\text{eq},2}||bK^{1/2}|f|^2||{0,K} = C_{\text{eq},2}||bK^{1/2}|f|f||{0,K}.$$

In case f ∈ P0(K) d , |f|f + g is a polynomial and the above arguments verify the assertion. □

Notation throughout the remaining parts of the proof. The following notation on the affine function f = a+Bξ for ξ = x−mid(K) applies throughout the proof. So a = f(mid(K)) and B = Df ∈ R d×d has the maximal singular value σ1(B) ≥ 0 (and σ1(B) = 0 iff B = 0).

The volume-bubble function b^K with maximal value 1 attained at x = mid(K) is continuous and vanishes outside K. The open set {x ∈ R d : 1/2 < bK(x)} lies compactly in the interior of the simplex K and contains its midpoint mid(K). Consequently there exists some positive r ^∗ < 1 such that the open ball ω around mid(K) with radius r ∗ is contained in {x ∈ R d : 1/2 < bK(x)}. Let r ^∗ denote the maximal radius with this property and let ω denote the associated open ball. Notice that r ∗ is a universal constant in the sense that it exclusively depends on the reference simplex K = Kref (since one we know the scaling invariance we can reduce the analysis to the reference element). Note that x ∈ Ω is equivalent to ξ = x − mid(K) ∈ R d satisfies |ξ| < r^∗ .

Proof of RHS = 0 ⇔ Df = 0 and g = Π0|f|f. (Note that we are proving that RHS = 0 ⇔ f, g ∈ P0(K) ^d and f = g). Let us denote by Π^k the L ² projection onto Pk(K) d . Suppose the upper bound in the assertion, called RHS, vanishes for the affine functions f and g; i.e. |f|f is affine in ω. If f(x) = 0 for all x ∈ ω then B = Df = 0. Otherwise |f| is positive at some x ∈ ω and so (by continuity) also in a small neighbourhood of x. We fix this point x throughout this proof, select a direction ζ ∈ R d , and consider

$$g(t) := |\boldsymbol{f}(\boldsymbol{x} + t\boldsymbol{\zeta})| \, \boldsymbol{f}(\boldsymbol{x} + t\boldsymbol{\zeta}) \cdot \boldsymbol{f}(\boldsymbol{x})$$

as a smooth function in the real parameter t near zero. Since |f|f is affine in ω, so is g(t) affine in t and its second derivative g ′′(0) = 0 vanishes. Abbreviate a := f(x) ̸= 0, b := Bζ, B = Df, and f := f(x + tζ) = a + tb and compute g(t) = |f|(a · f), g ′ (t) = (a · f) (b · f)/|f| + |f|(a · b), and eventually

$$g''(0) = |\boldsymbol{a}| |\boldsymbol{b}|^2 + (\boldsymbol{a} \cdot \boldsymbol{b})^2 / |\boldsymbol{a}|.$$

Since we have a ̸= 0, g ′′(0) = 0 implies b = Bζ = 0. Recall ζ is an arbitrary direction, Bζ = 0 follows for all ζ ∈ R d , whence B = 0. The remaining details are straightforward and omitted for brevity. □

Proof for ε|f(mid(K))| ≤ σ1(Df). Recall that the extra condition means that for some small and positive ε, the parameters a = f(mid(K)) and B = Df ∈ R d×d , that determine f ∈ P1(K) d , satisfy εα ≤ β for length α := |a| and the first singular value β := σ1(B). We may exclude the case β = 0 (as then 0 ≤ α ≤ β/ε = 0 means f = 0), when the assertion follows with Ceq,2 from classical inverse estimates mentioned above. In the remaining case β > 0, the assertion can be rewritten with the quotient

$$Q(\boldsymbol{f}) := \max_{\boldsymbol{g} \in \mathbb{P}1(K)^d} |||\boldsymbol{f}||\boldsymbol{f} + \boldsymbol{g}||{0,K} / ||bK^{1/2}(|\boldsymbol{f}||\boldsymbol{f} + \boldsymbol{g})||{0,K} \le C_{\text{eq}}.$$

The characterisation of RHS = 0 in the previous step reveals that Q(f) ≥ 1 (by b^K ≤ 1) is quotient of two positive terms (indeed, since β > 0 then B is not the zero matrix and f is not constant, whence the two terms in the quotient are positive). We rescale the data and observe Q(f) = Q(f/t) for any positive real t. The choice t := α + β > 0 reveals that we may and will assume α + β = 1 without loss of generality in the sequel and consider

$$\max{(\boldsymbol{a},\mathbf{B})\in\mathcal{C}(\varepsilon)}Q(\boldsymbol{f})=:C{\mathrm{eq}}(\varepsilon)$$

with the above notation f = a + Bξ and the compact parameter set

$$C(\varepsilon) = {(\boldsymbol{a}, \mathbf{B}) \in \mathbb{R}^{d \times (d+1)} : \alpha + \beta = 1 \text{ and } \varepsilon \alpha \leq \beta \text{ hold for } \alpha = |\boldsymbol{a}| \text{ and } \beta = \sigma_1(\mathbf{B})}.$$

Since C(ε) is compact and does not include B = 0 (because ε(1 − β) = εα ≤ β implies β ≥ ε/(1 + ε)), the quotient Q(f) = Q(a + Bξ) depends continuously on (a, B). Consequently, that maximum Ceq(ε) < ∞ of all Q(f) with f from the parameter set C(ε) is attained and in particular finite. (At this stage Ceq(ε) < ∞ might monotonically depend on ε > 0 with possibly C(ε) → ∞ as ε → 0.) □

Proof for σ1(Df) < ε0|f(mid(K))|. A rescaling shows that we may suppose a unit vector f(mid(K)) = a, |a| = 1, and we rename Df as t > 0 times B for a matrix B with largest singular value σ1(B) = 1. Hence the parameter regime translates into 0 < t < ε⁰ for a small ε⁰ < 1/2 defined below. The first ingredient in this proof is a Taylor expansion in the parameter t of |f|f. Since t is small and |a| = 1, |f| is differentiable and allows for a power series expansion. From f = a + tb and b = Bξ, we eventually infer |f| = 1 + (a · b)t + |b| ² − (a · b) 2 /2 + O(t 3 ) and thereafter

$$|\mathbf{f}|\mathbf{f} = \mathbf{a} + t(\mathbf{b} + (\mathbf{a} \cdot \mathbf{b}) \mathbf{a}) + t^2((|\mathbf{b}|^2 - (\mathbf{a} \cdot \mathbf{b})^2) \mathbf{a} + 2(\mathbf{a} \cdot \mathbf{b}) \mathbf{b})/2 + \delta$$

with some third-order remainder δ. This expansion holds for all x ∈ K and the difference δ ∈ L∞(K) d satisfies ∥δ∥L∞(K) ^d ≤ C0t 3 for all 0 < t < ε⁰ and ε⁰ < 1/2 sufficiently small. After abbreviating the terms in the displayed expansion of |f|f by a ∈ P0(K) d , g¹ := b + (a · b) a ∈ P1(K) d , and g2 := (|b| ² − (a · b) 2 ) a + 2(a · b) b ∈ P2(K) d , we have

$$|||\mathbf{f}||\mathbf{f} - \mathbf{a} - t\mathbf{g}_1 - t^2\mathbf{g}2/2||{L^{\infty}(K)^d} \le C_0 t^3.$$ (A.1)

The second ingredient is the b^K projection P : L 2 (K) ^d → P1(K) ^d defined, for any q ∈ L 2 (K) d , by

$$P\boldsymbol{q} \in \mathbb{P}_1(K)^d \quad \text{satisfies} \quad \int_K b_K(\boldsymbol{q} - P\boldsymbol{q}) \cdot \boldsymbol{\phi}_1 \; \mathrm{d}\boldsymbol{x} = 0 \quad \text{for all } \boldsymbol{\phi}_1 \in \mathbb{P}_1(K)^d.$$

Notice that P ∈ L(L 2 (K) d ) is a projection and the orthogonal projection in the b^K weighted L 2 (K) d scalar product (•, •)b^K := (bK•, •)K; P is merely an oblique projection in L 2 (K) d . Let ∥ • ∥b^K := (•, •) 1/2 b^K denote the induced b^K weighted L 2 (K) norm. The stability ∥Pq∥0,K ≤ Ceq,2∥q∥b^K is immediate from the above classical inverse estimate: C −2 eq,2∥Pq∥ 2 ⁰,K ≤ ∥Pq∥ 2 b^K = (Pq, q)b^K ≤ ∥Pq∥b^K ∥q∥b^K with a Cauchy inequality in the Hilbert space (L 2 (K) d ,(•, •)b^K ) in the last step, whence ∥Pq∥0,K ≤ Ceq,2∥q∥b^K .

In the second step of the proof we establish ∥(1 − P)(|f|f)∥0,K ≤ C1t 2 for 0 < t < ε0. A triangle inequality, the elementary estimate ∥ • ∥0,K ≤ |K| ¹/2∥ • ∥L∞(K) for the reference simplex K = Kref of volume |K| ≤ 1, and (A.1) result in

$$||(1-P)(|\boldsymbol{f}|\boldsymbol{f})||_{0,K} \le ||(1-P)(\boldsymbol{a}-t\boldsymbol{g}_1-t^2\boldsymbol{g}2/2)||{0,K} + C_0t^3.$$

Since (1 − P)(a − tg¹ ) = 0, the above stability ∥P • ∥0,K ≤ Ceq,2∥ • ∥b^K provides

$$2t^{-2}|(1-P)(\boldsymbol{a}-t\boldsymbol{g}_1-t^2\boldsymbol{g}2/2)|{0,K}=|(1-P)\boldsymbol{g}2|{0,K}\leq (1+C_{\text{eq},2})|\boldsymbol{g}2|{0,K}.$$

The point is that |a| = 1 = σ1(B) is bounded and so is b = Bξ with |ξ| < 1 (recall K = Kref is the reference simplex): Thus |b| = |Bξ| ≤ σ1(B) = 1 and eventually

$$|g_2| := |(|b|^2 - (a \cdot b)^2) a + 2(a \cdot b) b| \le ||b|^2 - (a \cdot b)^2| + 2|a \cdot b| \le 3$$ a.e. in $K$ .

The combination with the above estimates provides C¹ := 8(1 + Ceq,2) + C0/2 in

$$||(1-P)(|\mathbf{f}|\mathbf{f})||{0,K} \le 8t^2(1+C{\text{eq},2})+C_0t^3 \le C_1t^2.$$

The third step of the proof establishes C2t ² ≤ ∥b 1/2 ^K (1 − P)(|f|f)∥0,K. Recall b^K > 1/2 on the ball ω around mid(K) of radius r ∗ . Since σ1(B) = 1 we find some unit vector ζ with |Bζ| = 1 and hence at x := mid(K) + r ^∗ζ we have b = Bξ for ξ = x − mid(K) = r ^∗ζ and infer |b| = r ∗ |Bζ| = r ∗ . At other points in ω, |b| ≤ r ∗ , for instance b vanishes at mid(K). Altogether r ^∗ = ∥b∥L∞(ω) ^d . This and the calculation a · g² = |b| ² + (a · b) 2 imply

$$(r^*)^2 = |\boldsymbol{b}|_{L^{\infty}(\omega)^d}^2 \le |\boldsymbol{a} \cdot \boldsymbol{g}2|{L^{\infty}(\omega)} \le 2 \min_{g_0 \in \mathbb{R}} |\boldsymbol{a} \cdot \boldsymbol{g}_2 - g0|{L^{\infty}(\omega)}$$

with optimal value g⁰ = ∥a · g2∥L∞(ω)/2 in the last step (recall that a · g² ≥ 0 and it vanishes at mid(K)). The choice of g as the integral mean of a · g² over ω is therefore an upper bound,

$$\min_{g_0 \in \mathbb{R}} |\boldsymbol{a} \cdot \boldsymbol{g}_2 - g0|{L^{\infty}(\omega)} \leq |\boldsymbol{a} \cdot \boldsymbol{g}2 - \overline{g}|{L^{\infty}(\omega)} \leq C_{\text{inv}} |\boldsymbol{a} \cdot \boldsymbol{g}2 - \overline{g}|{L^2(\omega)}$$

from another inverse estimate (from equivalence of norms in P2(Ω)) in the last step. The key insight starts with the quadratic polynomial a · g² − g with integral mean zero over the symmetric (fixed) domain ω such that a·g2−g is symmetric along all straight lines through mid(K), while ξ is asymmetric (with respect to the centre mid(K)). Thus the product R ω (a · g² − g0)ξ dx = 0 vanishes. This means that a· g² −g⁰ has vanishing moments up to degree one: a· g² −g⁰ is L 2 (ω) perpendicular to all affine functions. Consequently

$$|\bm{a}\cdot\bm{g}2 - \overline{g}|{L^2(\omega)} = \min_{q_1 \in \mathbb{P}_1(\omega)} |\bm{a}\cdot\bm{g}_2 - q1|{L^2(\omega)} \leq 2 \min_{q_1 \in \mathbb{P}_1(\omega)} |\bm{a}\cdot\bm{g}_2 - q1|{bK} \leq 2 \min{q_1 \in \mathbb{P}_1(\omega)^d} |\bm{g}_2 - q1|{b_K}$$

with b^K > 1/2 a.e. in ω ⊂ K and |a| = 1 in the last steps. Notice that the last term minq1∈P1(ω) ^d ∥g2− q1∥b^K is equal to ∥(1 − P)g2∥b^K . The combination of the aforementioned inequalities reads

$$\frac{t^2(r^*)^2}{8C_{\text{inv}}} \le |t^2(1-P)\boldsymbol{g}2/2|{bK} \le |(1-P)(|\boldsymbol{f}|\boldsymbol{f})|{b_K} + C_0 t^3$$

with (A.1) and a triangle inequality in the last step. Under the condition

$$0 < t \le \varepsilon_0 := \min{1/2, (r^*)^2/(16C0C{\text{inv}})},\,$$

we conclude the proof of the assertion with C² = (r ∗ ) ²/(16C0Cinv) in the third step.

Step four concludes the proof for t ≤ ε0. The previous two steps provide

$$C_1^{-1}C2||(1-P)(|\mathbf{f}|\mathbf{f})||{0,K} \le t^2 C2 \le ||(1-P)(|\mathbf{f}|\mathbf{f})||{b_K}$$

Given any g ∈ P1(K) ^d and f as before, a triangle and the last inequality lead to

$$|||f|f + g||{0,K} \le ||(1 - P)(|f|f)||{0,K} + ||g + P(|f|f)||_{0,K}$$

$$\le C_1 C2^{-1} ||(1 - P)(|f|f)||{bK} + C{\text{eq},2} ||g + P(|f|f)||_{b_K}$$

with the aforementioned classical inverse estimate in the last step. A Cauchy inequality in R ^d and a Pythagoras identity in the b^K weighted Lebesgue norms based on the orthogonality ((1−P)(|f|f), g+ P(|f|f))b^K = 0 leads to

$$|||f|f + g||_{0,K} \le \sqrt{C_1^2 C2^{-2} + C{eq,2}^2} |||f|f + g||_{b_K}.$$

This concludes the proof of the assertion with C 2 eq = C 2 ¹C −2 ² + C 2 eq,2. □

Conclusion of the proof of the lemma. We depart from the previous definition of ε⁰ as a universal constant (it exclusively depends on K = Kref, whence solely on d) and see that the assertion follows with a constant Ceq = q C² 1C −2 ² + C² eq,2 (that exclusively depends on ^d) in case ^f ^∈ ^P1(K) d satisfies σ1(Df) < ε0|f(mid(K))|. The other regime is covered by the choice ε = ε⁰ and leads to the assertion with a constant Ceq(ε0) in case f ∈ P1(K) d satisfies ε0|f(mid(K))| ≤ σ1(Df). This concludes the proof. □

badianonaugmentedvelocityvorticity2025.md 179 KB Histórico Raw