Incorporating variable viscosity in vorticity-based formulations for Brinkman equations

TL;DR

This paper presents a novel mixed finite element formulation for Brinkman equations with variable viscosity, providing stability analysis and optimal error estimates, validated by computational experiments.

Contribution

It introduces a non-symmetric mixed finite element approach for Brinkman equations with variable viscosity, extending stability and error analysis to arbitrary order vorticity discretizations.

Findings

Stable finite element coupling for velocity, pressure, and vorticity.

Optimal a priori error estimates established.

Computational examples confirm theoretical results.

Abstract

In this brief note, we introduce a non-symmetric mixed finite element formulation for Brinkman equations written in terms of velocity, vorticity and pressure with non-constant viscosity. The analysis is performed by the classical Babu\v{s}ka-Brezzi theory, and we state that any inf-sup stable finite element pair for Stokes approximating velocity and pressure can be coupled with a generic discrete space of arbitrary order for the vorticity. We establish optimal a priori error estimates which are further confirmed through computational examples

Tables

Table 1. Table 1: Error history associated to the augmented scheme using ( 3.2 )-( 3.4 ) with k = ℓ = 1 𝑘 ℓ 1 k=\ell=1 .

DoF	$h$	${\Vert 𝒖-{𝒖}_h\Vert }_{1,\Omega }$	rate	${\Vert 𝝎-{𝝎}_h\Vert }_{0,\Omega }$	rate	${\Vert p-p_h\Vert }_{0,\Omega }$	rate
smooth viscosity ${\nu }_a$
84	0.7071	11.233	–	10.580	–	2126	–
284	0.3536	4.4150	1.347	3.6531	1.524	1194	0.832
1044	0.1768	1.2351	1.838	1.0024	1.863	271.24	2.136
4004	0.0884	0.3092	1.999	0.2482	2.016	44.490	2.609
15684	0.0442	0.0767	2.011	0.0609	2.027	6.2553	2.732
62084	0.0221	0.0191	2.005	0.0150	2.015	0.8594	2.525
247044	0.0111	0.0047	1.999	0.0037	2.008	0.2503	2.318
steeper viscosity ${\nu }_b$
84	0.7071	11.233	–	10.581	–	2125	–
284	0.3536	4.4150	1.347	3.6528	1.524	1193	0.832
1044	0.1768	1.2350	1.837	1.0024	1.862	271.25	2.136
4004	0.0884	0.3093	1.998	0.2484	2.016	44.491	2.609
15684	0.0442	0.0767	2.011	0.0609	2.027	6.2553	2.731
62084	0.0221	0.0191	2.005	0.0151	2.016	0.8603	2.437
247044	0.0111	0.0048	1.999	0.0037	2.008	0.2487	2.290

Equations

$$\nu {\mathbf{K}}^{-1}\mathbf{u}+\nu \mathrm{curl}\mathbf{\omega }-2\mathbf{\epsilon }(\mathbf{u})\nabla \nu +\nabla p$$

$$\mathbf{\omega }-\mathrm{curl}\mathbf{u}$$

$$\mathrm{div}\mathbf{u}$$

$$\mathbf{u}$$

$$(p,1{)}_{0,\Omega }$$

$$0<{\nu }_0\le \nu \le {\nu }_1.$$

$${\sigma }_{\min }|\mathbf{v}{|}^2\le {\mathbf{v}}^t{\mathbf{K}}^{-1}\mathbf{v}\le {\sigma }_{\max }|\mathbf{v}{|}^2\forall \mathbf{v}\in {\mathbb{R}}^3.$$

$$-2\mathrm{div}(\nu \mathbf{\epsilon }(\mathbf{u}))=-2\nu \mathrm{div}(\mathbf{\epsilon }(\mathbf{u}))-2\mathbf{\epsilon }(\mathbf{u})\nabla \nu =-\nu \Delta \mathbf{u}-2\mathbf{\epsilon }(\mathbf{u})\nabla \nu =\nu \mathrm{curl}(\mathrm{curl}\mathbf{u})-2\mathbf{\epsilon }(\mathbf{u})\nabla \nu ,$$

$$\mathrm{curl}(\mathrm{curl}\mathbf{v})=-\Delta \mathbf{v}+\nabla (\mathrm{div}\mathbf{v}).$$

$$|||\mathbf{v}||{|}_{1,\Omega }^2:=\Vert \mathbf{v}{\Vert }_{0,\Omega }^2+\Vert \mathrm{curl}\mathbf{v}{\Vert }_{0,\Omega }^2+\Vert \mathrm{div}\mathbf{v}{\Vert }_{0,\Omega }^2.$$

$$\Vert \mathbf{v}{\Vert }_{1,\Omega }^2\le C_{pf}(\Vert \mathrm{curl}\mathbf{v}{\Vert }_{0,\Omega }^2+\Vert \mathrm{div}\mathbf{v}{\Vert }_{0,\Omega }^2)\forall \mathbf{v}\in {\mathrm{H}}_0^1(\Omega {)}^3,$$

$${\int }_{\Omega }\mathrm{curl}\mathbf{\omega }\cdot \mathbf{v}={\int }_{\Omega }\mathbf{\omega }\cdot \mathrm{curl}\mathbf{v}+⟨\mathbf{\omega }\times \mathbf{n},\mathbf{v}{⟩}_{\partial \Omega },$$

$${\int }_{\Omega }\nu {\mathbf{K}}^{-1}\mathbf{u}\cdot \mathbf{v}-2{\int }_{\Omega }\mathbf{\epsilon }(\mathbf{u})\nabla \nu \cdot \mathbf{v}+{\int }_{\Omega }\nu \mathbf{\omega }\cdot \mathrm{curl}\mathbf{v}+{\int }_{\Omega }\mathbf{\omega }\cdot (\nabla \nu \times \mathbf{v})-{\int }_{\Omega }p\mathrm{div}\mathbf{v}$$

$${\int }_{\Omega }\nu \mathbf{\theta }\cdot \mathrm{curl}\mathbf{u}-{\int }_{\Omega }\nu \mathbf{\omega }\cdot \mathbf{\theta }$$

$$-{\int }_{\Omega }q\mathrm{div}\mathbf{u}$$

$${\kappa }_1{\nu }_0{\int }_{\Omega }(\mathrm{curl}\mathbf{u}-\mathbf{\omega })\cdot \mathrm{curl}\mathbf{v}$$

$${\kappa }_2{\int }_{\Omega }\mathrm{div}\mathbf{u}\mathrm{div}\mathbf{v}$$

$$\begin{array}{rl}A((\mathbf{u},\mathbf{\omega }),(\mathbf{v},\mathbf{\theta }))+B((\mathbf{v},\mathbf{\theta }),p)=& G(\mathbf{v},\mathbf{\theta })\forall (\mathbf{v},\mathbf{\theta })\in {\mathrm{H}}_0^1(\Omega {)}^3\times {\mathrm{L}}^2(\Omega {)}^3,\\ B((\mathbf{u},\mathbf{\omega }),q)=& 0\forall q\in {\mathrm{L}}_0^2(\Omega ),\end{array}$$

$$A((\mathbf{u},\mathbf{\omega }),(\mathbf{v},\mathbf{\theta }))$$

$$+{\kappa }_2{\int }_{\Omega }\mathrm{div}\mathbf{u}\mathrm{div}\mathbf{v}-{\kappa }_1{\nu }_0{\int }_{\Omega }\mathbf{\omega }\cdot \mathrm{curl}\mathbf{v}-2{\int }_{\Omega }\mathbf{\epsilon }(\mathbf{u})\nabla \nu \cdot \mathbf{v}+{\int }_{\Omega }\mathbf{\omega }\cdot (\nabla \nu \times \mathbf{v}),$$

$$B((\mathbf{v},\mathbf{\theta }),q)$$

$$|{\int }_{\Omega }\nu {\mathbf{K}}^{-1}\mathbf{u}\cdot \mathbf{v}|\le {\sigma }_{\max }{\nu }_1|||\mathbf{u}||{|}_{1,\Omega }|||\mathbf{v}||{|}_{1,\Omega },|{\int }_{\Omega }\nu \mathbf{\omega }\cdot \mathbf{\theta }|\le {\nu }_1\Vert \mathbf{\omega }{\Vert }_{0,\Omega }\Vert \mathbf{\theta }{\Vert }_{0,\Omega },$$

$$|{\int }_{\Omega }\nu \mathbf{\theta }\cdot \mathrm{curl}\mathbf{v}|\le {\nu }_1\Vert \mathbf{\theta }{\Vert }_{0,\Omega }|||\mathbf{v}||{|}_{1,\Omega },|{\int }_{\Omega }\mathbf{\epsilon }(\mathbf{u})\nabla \nu \cdot \mathbf{v}|\le \Vert \nabla \nu {\Vert }_{\infty ,\Omega }\Vert \mathbf{\epsilon }(\mathbf{u}){\Vert }_{0,\Omega }\Vert \mathbf{v}{\Vert }_{0,\Omega },$$

$$|{\int }_{\Omega }\mathbf{\theta }\cdot (\nabla \nu \times \mathbf{v})|\le 2\Vert \nabla \nu {\Vert }_{\infty ,\Omega }\Vert \mathbf{v}{\Vert }_{0,\Omega }\Vert \mathbf{\theta }{\Vert }_{0,\Omega },|G(\mathbf{v},\mathbf{\theta })|\le \Vert \mathbf{f}{\Vert }_{0,\Omega }\Vert \mathbf{v}{\Vert }_{0,\Omega }.$$

$$|A((\mathbf{u},\mathbf{\omega }),(\mathbf{v},\mathbf{\theta }))|\le C_1\Vert (\mathbf{u},\mathbf{\omega }){\Vert }_{{\mathrm{H}}_0^1(\Omega {)}^3\times {\mathrm{L}}^2(\Omega {)}^3}\Vert (\mathbf{v},\mathbf{\theta }){\Vert }_{{\mathrm{H}}_0^1(\Omega {)}^3\times {\mathrm{L}}^2(\Omega {)}^3},$$

$$|B((\mathbf{v},\mathbf{\theta }),q)|\le C_2\Vert (\mathbf{v},\mathbf{\theta }){\Vert }_{{\mathrm{H}}_0^1(\Omega {)}^3\times {\mathrm{L}}^2(\Omega {)}^3}\Vert q{\Vert }_{0,\Omega },|G(\mathbf{v},\mathbf{\theta })|\le C_3\Vert (\mathbf{v},\mathbf{\theta }){\Vert }_{{\mathrm{H}}_0^1(\Omega {)}^3\times {\mathrm{L}}^2(\Omega {)}^3},$$

$$\Vert (\mathbf{v},\mathbf{\theta }){\Vert }_{{\mathrm{H}}_0^1(\Omega {)}^3\times {\mathrm{L}}^2(\Omega {)}^3}^2=|||\mathbf{v}||{|}_{1,\Omega }^2+\Vert \mathbf{\theta }{\Vert }_{0,\Omega }^2.$$

$$\frac{4\Vert \nabla \nu {\Vert }_{\infty ,\Omega }^2}{{\sigma }_{\min }{\nu }_0^2}<1/4.$$

$$A((\mathbf{v},\mathbf{\theta }),(\mathbf{v},\mathbf{\theta }))\ge \alpha \Vert (\mathbf{v},\mathbf{\theta }){\Vert }_{{\mathrm{H}}_0^1(\Omega {)}^3\times {\mathrm{L}}^2(\Omega {)}^3}^2\forall (\mathbf{v},\mathbf{\theta })\in {\mathrm{H}}_0^1(\Omega {)}^3\times {\mathrm{L}}^2(\Omega {)}^3.$$

$$\left|2{\int }_{\Omega }\mathbf{\epsilon }(\mathbf{v})\nabla \nu \cdot \mathbf{v}\right|\le \frac{4\Vert \nabla \nu {\Vert }_{\infty ,\Omega }^2}{{\sigma }_{\min }{\nu }_0}(\Vert \mathrm{curl}\mathbf{v}{\Vert }_{0,\Omega }^2+\Vert \mathrm{div}\mathbf{v}{\Vert }_{0,\Omega }^2)+\frac{{\sigma }_{\min }{\nu }_0}{4}\Vert \mathbf{v}{\Vert }_{0,\Omega }^2,$$

$$\left|{\int }_{\Omega }\mathbf{\theta }\cdot (\nabla \nu \times \mathbf{v})\right|\le$$

$$\left|{\kappa }_1{\nu }_0{\int }_{\Omega }\mathbf{\theta }\cdot \mathrm{curl}\mathbf{v}\right|\le$$

$$A((\mathbf{v},\mathbf{\theta }),(\mathbf{v},\mathbf{\theta }))\ge$$

$$+{\kappa }_2\Vert \mathrm{div}\mathbf{v}{\Vert }_{0,\Omega }^2-\frac{4\Vert \nabla \nu {\Vert }_{\infty ,\Omega }^2}{{\sigma }_{\min }{\nu }_0}(\Vert \mathrm{curl}\mathbf{v}{\Vert }_{0,\Omega }^2+\Vert \mathrm{div}\mathbf{v}{\Vert }_{0,\Omega }^2)-\frac{{\sigma }_{\min }{\nu }_0}{4}\Vert \mathbf{v}{\Vert }_{0,\Omega }^2$$

$$-\frac{4\Vert \nabla \nu {\Vert }_{\infty ,\Omega }^2}{{\sigma }_{\min }{\nu }_0}\Vert \mathbf{\theta }{\Vert }_{0,\Omega }^2-\frac{{\sigma }_{\min }{\nu }_0}{4}\Vert \mathbf{v}{\Vert }_{0,\Omega }^2$$

$$\ge$$

$$+\left({\kappa }_2-\frac{4\Vert \nabla \nu {\Vert }_{\infty ,\Omega }^2}{{\sigma }_{\min }{\nu }_0}\right)\Vert \mathrm{div}\mathbf{v}{\Vert }_{0,\Omega }^2.$$

$$A((\mathbf{v},\mathbf{\theta }),(\mathbf{v},\mathbf{\theta }))\ge$$

$$\ge$$

Full text

Keywords: Brinkman equations; vorticity-based formulation; mixed finite elements; variable viscosity; error analysis.

Mathematics subject classifications (2000): 65N30; 65N12; 76D07; 65N15.

1 Introduction

Formulations for flow equations that use vorticity as an additional unknown enjoy many appealing features [21], and starting from the works [10, 12], they have been employed in many instances (see e.g. [1, 2, 4, 3, 6, 20, 11, 24, 5, 23]). However, a major limitation in all of these contributions, in comparison with competing formulations using solely the primal variables, is that the transformation of the momentum equation introducing vorticity (and subsequently using a convenient structure of the problem to analyse its mathematical properties and devising suitable numerical schemes) is only valid when the viscosity is constant. Plus, a number of applications including Stokes flow and coupled thermal or thermo-haline effects with Brinkman flows (see e.g. [16, 19, 22] and [17, 18, 25], respectively) depend strongly on marked spatial distributions of viscosity.

In this brief note, we provide a way of incorporating variable viscosities while keeping vorticity as field variable. The resulting non-symmetric formulation is augmented via least-squares terms involving the constitutive equation and mass conservation equation and subsequently the problem maintains a saddle-point structure amenable to analysis through classical tools from mixed methods (under the assumption that the viscosity is regular enough). Even if we have decided to provide all steps for the specific case of Brinkman equations, the same ideas in principle carry over to other vorticity-based models such as Oseen, Navier-Stokes, interfacial flows, and coupled Boussinesq or flow-transport problems.

The main advantages of the propose scheme are the direct approximation of vorticity without invoking any postprocessing, and also the simplicity of the analysis and implementation. Indeed, one can use standard inf-sup stable finite elements for the Stokes equations plus any conforming discrete space for vorticity.

Outline. In Section 2, we recall the governing equations and state the least-squares–based augmented formulation. There we also perform the solvability analysis employing standard arguments from the Babuška–Brezzi theory. The finite element discretisation is presented in Section 3, where we also write a stability analysis and derive optimal error estimates. A few numerical tests illustrating the convergence of the proposed method are finally reported in Section 4.

2 Variable viscosity Brinkman equations

Let $\Omega$ be a bounded domain of $\mathbb{R}^{3}$ with Lipschitz boundary $\Gamma=\partial\Omega$, and let us write the following version of the Brinkman equations with variable viscosity where the unknowns are velocity $\boldsymbol{u}$, vorticity $\boldsymbol{\omega}$, and pressure $p$ of the incompressible viscous fluid

[TABLE]

The kinematic viscosity is assumed such that $\nu\in\mathrm{W}^{1,\infty}(\Omega)$ and

[TABLE]

Moreover, $\boldsymbol{f}\in\mathrm{L}^{2}(\Omega)^{3}$ is a force density and $\mathbf{K}\in\mathrm{L}^{\infty}(\Omega)^{3\times 3}$ is the (symmetric and uniformly positive definite) tensor of permeability. In particular, there exist $\sigma_{\min},\sigma_{\max}>0$ such that

[TABLE]

Instead of $\nu\mathbf{K}^{-1}$ some works equivalently use $\hat{\mathbf{K}}^{-1}$ as the drag coefficient in the momentum equation, where $\hat{\mathbf{K}}=\mathbf{K}/\nu$. Note that (2.1) can be derived from the usual momentum equation by invoking the identity

[TABLE]

where $\boldsymbol{\varepsilon}(\boldsymbol{u})$ is the strain rate tensor and where we have also used (2.3) and the additional identity

[TABLE]

2.1 Variational formulation and preliminary results

For any $s\geq 0$, the notation $\left\|\cdot\right\|_{s,\Omega}$ stands for the norm of the Hilbertian Sobolev spaces $\mathrm{H}^{s}(\Omega)$ or $\mathrm{H}^{s}(\Omega)^{3}$, with the usual convention $\mathrm{H}^{0}(\Omega):=\mathrm{L}^{2}(\Omega)$. We also endow the space $\mathrm{H}_{0}^{1}(\Omega)^{3}$ with the following norm:

[TABLE]

We note that in $\mathrm{H}_{0}^{1}(\Omega)^{3}$, the above norm is equivalent to the usual norm. In particular, we have that there exists a positive constant $C_{pf}$ such that:

[TABLE]

the above inequality is a consequence of the identity $\|\nabla\boldsymbol{v}\|_{0,\Omega}^{2}=\|\mathop{\mathbf{curl}}\nolimits\boldsymbol{v}\|_{0,\Omega}^{2}+\|\mathop{\mathrm{div}}\nolimits\boldsymbol{v}\|_{0,\Omega}^{2}$ which follows from (2.7) and the Poincaré inequality.

Testing (2.1)-(2.3) appropriately, using Green’s formula in the following version (see [14, Thm. 2.11])

[TABLE]

and applying the boundary conditions (2.4)-(2.5), we get the following weak formulation

[TABLE]

where $\mathrm{L}_{0}^{2}(\Omega):=\{q\in\mathrm{L}^{2}(\Omega):(q,1)_{0,\Omega}=0\}$. Then, we proceed to augment this formulation with the following residual terms arising from equations (2.2) and (2.3):

[TABLE]

with $\nu_{0}>0$ (cf. (2.6)), and where $\kappa_{1}$ and $\kappa_{2}$ are positive parameters to be specified later. Then, the augmented formulation reads: Find $((\boldsymbol{u},\boldsymbol{\omega}),p)\in(\mathrm{H}_{0}^{1}(\Omega)^{3}\times\mathrm{L}^{2}(\Omega)^{3})\times\mathrm{L}_{0}^{2}(\Omega)$ such that

[TABLE]

where the bilinear forms and the linear functional are defined by

[TABLE]

for all $(\boldsymbol{u},\boldsymbol{\omega}),(\boldsymbol{v},\boldsymbol{\theta})\in\mathrm{H}_{0}^{1}(\Omega)^{3}\times\mathrm{L}^{2}(\Omega)^{3}$, and $q\in\mathrm{L}_{0}^{2}(\Omega)$.

2.2 Unique solvability of the augmented formulation

Problem (2.11) accommodates an analysis directly under the classical Babuška-Brezzi theory [9, 13]. More precisely, the continuity of the bilinear and linear functionals in (2.12) is a direct consequence of Lemma 2.1 below, whose proof is obtained by rather standard arguments. In particular, the penultimate estimate holds owing to the assumption $\nabla\nu\in\mathrm{L}^{\infty}(\Omega)^{3}$ and the fact that $\|\nabla\nu\times\boldsymbol{v}\|_{0,\Omega}\leq 2\|\nabla\nu\|_{\infty,\Omega}\|\boldsymbol{v}\|_{0,\Omega}$. Then, the ellipticity of $A$, stated in Lemma 2.2, follows from adding the redundant terms in (2.9)-(2.10).

Lemma 2.1**.**

The following estimates hold

[TABLE]

Therefore, we have that there exist $C_{1},C_{2},C_{3}>0$ such that

[TABLE]

where

[TABLE]

Lemma 2.2**.**

Assume that

[TABLE]

Suppose that $\kappa_{1}\in(\frac{1}{2},\frac{3}{2})$ and $\kappa_{2}>\frac{\nu_{0}}{4}$. Then, there exists $\alpha>0$ such that

[TABLE]

Proof.

Given $(\boldsymbol{v},\boldsymbol{\theta})\in\mathrm{H}_{0}^{1}(\Omega)^{3}\times\mathrm{L}^{2}(\Omega)^{3}$ first we observe that as a consequence of Lemma 2.1, we have

[TABLE]

where we have used (2.8). Moreover, using that $\|(\nabla\nu\times\boldsymbol{v})\|_{0,\Omega}\leq 2\|\nabla\nu\|_{\infty,\Omega}\|\boldsymbol{v}\|_{0,\Omega}$, we obtain

[TABLE]

and these estimates are put in combination with Cauchy-Schwarz inequality to obtain that

[TABLE]

Now, using (2.13), we have that

[TABLE]

where $\alpha$ depends on $\kappa_{1},\kappa_{2},\nu_{0}$ and $\sigma_{\min}$. ∎

Finally, recall the inf-sup condition (cf. [13]): there exists $C>0$ only depends on $\Omega$ such that

[TABLE]

Lemma 2.3**.**

There exists $\beta>0$, independent of $\nu$, such that

[TABLE]

Proof.

The result is a consequence of (2.14) and the fact that

[TABLE]

where the term in the righ-hand side has the usual norm in $\mathrm{H}_{0}^{1}(\Omega)^{3}$. ∎

All these steps lead to the unique solvability of the problem.

Theorem 2.1**.**

There exists a unique solution $((\boldsymbol{u},\boldsymbol{\omega}),p)\in(\mathrm{H}_{0}^{1}(\Omega)^{3}\times\mathrm{L}^{2}(\Omega)^{3})\times\mathrm{L}_{0}^{2}(\Omega)$ to (2.11) and there exists a constant $C>0$ such that the following continuous dependence result holds:

[TABLE]

Proof.

By virtue of Lemmas 2.2 and 2.3, the proof is a straightforward application of [9, Thm. II.1.1]. ∎

3 Finite element discretisations

Taking generic subspaces for the approximation of velocity, vorticity, and pressure, a Galerkin scheme associated with (2.11) reads: Find $((\boldsymbol{u}_{h},\boldsymbol{\omega}_{h}),p_{h})\in(\mathrm{H}_{h}\times\mathrm{Z}_{h})\times\mathrm{Q}_{h}$ such that

[TABLE]

We can adopt in particular

[TABLE]

where $k\geq 1,\ell\geq 0$. Here $\{\mathcal{T}_{h}(\Omega)\}_{h>0}$ is a shape-regular family of partitions of $\bar{\Omega}$ by tetrahedra $T$ of diameter $h_{T}$. The meshsize is $h:=\max\{h_{T}:\;T\in\mathcal{T}_{h}(\Omega)\}$, and $\mathbb{P}_{m}(S)$ denotes the space of polynomials with total degree up to $m$, defined on a generic set $S$.

We recall that $\mathrm{H}_{h}\times\mathrm{Q}_{h}$ in the generalised Hood-Taylor finite element pair for the Stokes equations [15]. As we will see, the schemes coming from (3.1)-(3.4) are well-posed for any approximation order of the discrete vorticity $\ell$ (and being continuous or discontinuous polynomials); however, an appropriate choice is to take $\ell=k$ and discontinuous elements, which deliver a consistent overall rate of convergence for all unknowns.

Next, we proceed to show that the proposed method is stable and convergent.

Lemma 3.1**.**

Assuming (2.13), and choosing $\kappa_{1}\in\left(\frac{1}{2},\frac{3}{2}\right)$ and $\kappa_{2}>\frac{\nu_{0}}{4}$, there exists $\alpha>0$, such that

[TABLE]

Remark 3.1**.**

The values for the augmentation parameters $\kappa_{1}$ and $\kappa_{2}$ are chosen such that the largest ellipticity constant in Lemma 3.1 is achieved. This means that we take $\kappa_{1}=1$ (the middle point of the relevant interval, see e.g. [4, Sect. 3]) and $\kappa_{2}=\frac{\nu_{0}}{2}$.

Moreover, since for the pair of spaces (3.2),(3.4) one has an inf-sup condition of the form

[TABLE]

where $\tilde{\beta}_{2}$ is independent of $h$ (see [7, 8]), then it is straightforward to prove the following result.

Lemma 3.2**.**

There exists $\tilde{\beta}>0$, such that

[TABLE]

Recall now that the Lagrange interpolant $\Pi:\mathrm{H}^{1+s}(\Omega)^{3}\to\mathrm{H}_{h}$ satisfies the following error estimate: There exists $C>0$, independent of $h$, such that for all $s\in(1/2,k+1]$:

[TABLE]

Likewise, denoting by ${\mathcal{P}}$ the orthogonal projection from $\mathrm{L}^{2}(\Omega)$ (or from $\mathrm{L}^{2}(\Omega)^{3}$) onto the subspace $\mathrm{Q}_{h}$ (or onto the subspace $\mathrm{Z}_{h}$), we have an estimate valid for all $s>0$:

[TABLE]

Thanks to Lemmas 3.1 and 3.2, we can state the stability and Céa estimate of the method as follows.

Theorem 3.1**.**

Let $\mathrm{H}_{h}$, $\mathrm{Z}_{h}$ and $\mathrm{Q}_{h}$ be specified as in (3.2), (3.3) and (3.4), respectively. Then, there exists a unique $((\boldsymbol{u}_{h},\boldsymbol{\omega}_{h}),p_{h})\in(\mathrm{H}_{h}\times\mathrm{Z}_{h})\times\mathrm{Q}_{h}$ solution of the Galerkin scheme (3.1). Furthermore, there exist positive constants $\hat{C}_{1},\,\hat{C}_{2}>0$, independent of $h$, such that

[TABLE]

and

[TABLE]

where $((\boldsymbol{u},\boldsymbol{\omega}),p)\in(\mathrm{H}_{0}^{1}(\Omega)^{3}\times\mathrm{L}^{2}(\Omega)^{3})\times\mathrm{L}_{0}^{2}(\Omega)$ is the unique solution to variational problem (2.11).

And finally the convergence of the augmented scheme can be formulated as follows.

Theorem 3.2**.**

Let $\mathrm{H}_{h},\mathrm{Z}_{h}$ and $\mathrm{Q}_{h}$ be given by (3.2), (3.3), and (3.4), respectively, setting $\ell=k$ with $k\geq 1$. Let $(\boldsymbol{u},\boldsymbol{\omega},p)\in\mathrm{H}_{0}^{1}(\Omega)^{3}\times\mathrm{L}^{2}(\Omega)^{3}\times\mathrm{L}_{0}^{2}(\Omega)$ and $(\boldsymbol{u}_{h},\boldsymbol{\omega}_{h},p_{h})\in\mathrm{H}_{h}\times\mathrm{Z}_{h}\times\mathrm{Q}_{h}$ be the unique solutions to the continuous and discrete problems (2.11) and (3.1), respectively. Assume that $\boldsymbol{u}\in\mathrm{H}^{1+s}(\Omega)^{3}$, $\boldsymbol{\omega}\in\mathrm{H}^{s}(\Omega)^{3}$ and $p\in\mathrm{H}^{s}(\Omega)$, for some $s\in(1/2,k+1]$. Then, there exists $\hat{C}>0$, independent of $h$, such that

[TABLE]

Proof.

It follows from (3.8)-(3.9) and (3.6)-(3.7). ∎

Remark 3.2**.**

Instead of Hood–Taylor finite elements (3.2),(3.4), we can also consider any other Stokes inf-sup stable pairs. For instance, using the MINI-element for velocity and pressure (piecewise linear velocities enriched with quartic bubbles, or cubic bubbles in 2D, and piecewise linear and continuous pressures, see e.g. [9]) and piecewise constant elements for vorticity, we can easily adapt the analysis to obtain the error estimate

[TABLE]

4 Numerical results

We proceed to verify numerically the convergence rates predicted by (3.10). Following [16], on $\Omega=(0,1)^{2}$ we take $\phi(x,y)=1000x^{2}(1-x)^{4}y^{3}(1-y)^{2}$ and define exact velocity, vorticity, and pressure as

[TABLE]

which satisfy the incompressibility constraint as well as the homogeneous boundary and compatibility conditions. Two specifications for viscosity are considered, with a mild and with a higher gradient

[TABLE]

and we use $\nu_{0}=10^{-4}$, $\nu_{1}=1$. A current restriction in our analysis is (2.13) that only permits sufficiently small permeability such that the lower bound for its inverse, $\sigma_{\min}$ is large enough (in any case, for most relevant applications in porous media flow these values are reasonable). We use $\mathbf{K}=10^{-6}\mathbf{I}$. Sample solutions are shown in Figure 1 and the convergence history (produced on a sequence of successively refined meshes and computing errors for all fields and rates as usual) is presented in Table 1. At least for these two cases, we observe a higher convergence of the pressure and that a steeper viscosity does not affect the accuracy.

We close with a 3D example simulating the cavity flow in the presence of a viscosity boundary layer. The domain $\Omega=(0,1)^{3}$ is discretised with a structured tetrahedral mesh and we employ the scheme from Remark 3.2 (the MINI-element for the velocity-pressure pair together with piecewise constant vorticity approximation) resulting in a system with 560165 DoF. We use $\boldsymbol{f}=\boldsymbol{0}$ and the velocity $\boldsymbol{u}=(1,0,0)^{\tt t}$ is prescribed on the top lid (at $z=1$) while no-slip velocities are set on the other sides of the boundary. We set $\mathbf{K}=10^{-4}\mathbf{I}$, and choosing now $\nu_{0}=10^{-5}$, $\nu_{1}=10$, the variable viscosity field is

[TABLE]

The approximate solutions are depicted in Figure 2 where we observe how the velocity and pressure lose the usual symmetry expected in lid-driven cavity flows, and it separates due to the viscosity boundary layer.