A cholera mathematical model with vaccination and the biggest outbreak of world’s history

1.   Introduction

Scope and related work. Brinkman equations constitute one of the simplest homogenized models for viscous fluid flow in highly heterogeneous porous media. The analysis of the underlying PDE system exhibits challenges due to the presence of two parameters: viscosity and permeability. Likewise, the design and analysis of numerical methods that maintain robustness with respect to the very relevant cases of low effective viscosity (Darcy-like regime) vs very large permeability (Stokes-like regime), is still an issue of high interest.

For classical velocity–pressure formulations, a number of contributions are available to deal with the inf-sup compatibility of velocity and pressure spaces irrespective of the regime (Darcy-like or Stokes-like). See, e.g., [17,18,20,21,23]. On the other hand, mixed finite element formulations for Brinkman equations (that is, using other fields apart from the classical velocity–pressure pair) include pseudostress–based methods [14,19] (see also VEM counterparts in [8,15], and [26] for the DG case with strong imposition of pseudostress symmetry, closer to the present contribution), vorticity–velocity–pressure schemes in augmented and non-augmented form [3-6,28]. There, the analysis of continuous and discrete problems depends either on the Babuška–Brezzi theory, or on a generalized inf-sup argument.

The method advanced in this paper is based on a pure–stress formulation, obtained from taking the deviatoric part of the stress and eliminating pressure, and using the momentum equation to write velocity as a function of the external and internal forces (forcing plus the divergence of the Cauchy stress). We consider the case of mixed boundary conditions, and use a symmetric interior penalty DG discretization. The space for discrete stresses incorporates symmetry strongly, as in [26] (see also [29] for elasticity and [25] for viscoelasticity). At the discrete level, the velocity and the pressure are easily recovered in terms of the discrete Cauchy stress through post-processing. While the pressure post-processing is derivative-free, the usual post-process of velocity through the momentum balance entails taking the discrete divergence of the approximate stress. We also propose a more involved reconstruction, starting from the initial post-processed velocity and solving an additional problem in mixed form using $H( \mathop{\mathrm{div}}\nolimits)$-conforming approximations. The newly computed velocity is exactly divergence-free and we prove that it enjoys the same convergence order as the stress approximation.

Other advantages of the present formulation include the physically accurate imposition of outflow boundary conditions fixing the normal traces of the full Cauchy stress (which is not straightforward to incorporate with vorticity- or pseudostress-based methods), having symmetric and positive definite linear systems after discretization, and straightforwardly handling heterogeneous and anisotropic permeability distributions. In this regard, the analysis of the formulation uses a DG-energy norm chosen in such a way that the error estimates are independent of the jumps of the permeability. This implies a similar property concerning the estimates for the velocity approximation. However, for the pressure we can only achieve stability and optimal rates of convergence in terms of the mesh parameter.

Outline. The contents of the paper have been laid out in the following manner. The remainder of this Section lists useful notation regarding tensors, and it recalls the definition of Sobolev spaces, inner products, and integration by parts. In Section 2 we state the Brinkman problem, specify assumptions on the model coefficients, and derive a weak formulation written only in terms of stress. Section 3 gathers the preliminary concepts needed for the construction and analysis of the discrete problem. Here we include mesh properties, recall useful trace and inverse inequalities, define suitable projectors, and establish approximation properties in conveniently chosen norms. The definition of the stress-based DG method and the proof its consistency and optimal convergence in the stress variable is presented in Section 4. The velocity and pressure post-processing, together with the corresponding error estimates, are given in Section 5. Next, in Section 6 we derive optimal error bounds for the stress in the $L^2$-norm, and Section 7 contains a collection of numerical examples that confirm experimentally the convergence of the method, and also showcase its application into typical flow problems in 2D and 3D.

Notational preliminaries. We denote the space of real matrices of order $d\times d$ by $\mathbb{M}$ and let $\mathbb{S}: = \left\{ {{ \boldsymbol{\tau}\in \mathbb{M};\ \boldsymbol{\tau} = \boldsymbol{\tau}^{ \mathtt{t}} }} \right\}$ be the subspace of symmetric matrices, where $\boldsymbol{\tau}^{ \mathtt{t}}: = (\tau_{ji})$ stands for the transpose of $\boldsymbol{\tau} = (\tau_{ij})$. The component-wise inner product of two matrices $\boldsymbol{\sigma}, \, \boldsymbol{\tau} \in \mathbb{M}$ is defined by $\boldsymbol{\sigma}: \boldsymbol{\tau}: = \sum_{i,j}\sigma_{ij}\tau_{ij}$. We also introduce the deviatoric part $\boldsymbol{\tau}^{\mathtt{D}}: = \boldsymbol{\tau}-\frac{1}{d}\left( \mathop{\mathrm{tr}}\nolimits \boldsymbol{\tau}\right) \boldsymbol{I}$ of a tensor $\boldsymbol{\tau}$, where $\mathop{\mathrm{tr}}\nolimits \boldsymbol{\tau}: = \sum_{i = 1}^d\tau_{ii}$ and $\boldsymbol{I}$ stands here for the identity in $\mathbb{M}$.

Let $\Omega$ be a polyhedral Lipschitz bounded domain of $\mathbb{R}^d$ $(d = 2,3)$, with boundary $\partial \Omega$. Along this paper we convene to apply all differential operators row-wise. Hence, given a tensorial function $\boldsymbol{\sigma}:\Omega\to \mathbb{M}$ and a vector field $\boldsymbol{u}:\Omega\to \mathbb{R}^d$, we set the divergence $\mathop{\bf{div}}\nolimits \boldsymbol{\sigma}:\Omega \to \mathbb{R}^d$, the gradient $\boldsymbol{\nabla} \boldsymbol{u}:\Omega \to \mathbb{M}$, and the linearized strain tensor $\boldsymbol{\varepsilon}( \boldsymbol{u}) : \Omega \to \mathbb{S}$ as

For $s\in \mathbb{R}$, $H^s(\Omega,E)$ stands for the usual Hilbertian Sobolev space of functions with domain $\Omega$ and values in $E$, where $E$ is either $\mathbb{R}$, $\mathbb{R}^d$ or $\mathbb{M}$. In the case $E = \mathbb{R}$ we simply write $H^s(\Omega)$. The norm of $H^s(\Omega,E)$ is denoted $\left\| {{\cdot}} \right\|_{s,\Omega}$ and the corresponding semi-norm $|\cdot|_{s,\Omega}$, indistinctly for $E = \mathbb{R}, \mathbb{R}^d, \mathbb{M}$. We use the convention $H^0(\Omega, E): = L^2(\Omega,E)$ and let $(\cdot, \cdot)$ be the inner product in $L^2(\Omega, E)$, for $E = \mathbb{R}, \mathbb{R}^d, \mathbb{M}$, namely,

The space $H( \mathop{\mathrm{div}}\nolimits, \Omega)$ stands for the vector fields $\boldsymbol{v}\in L^2(\Omega, \mathbb{R}^d)$ satisfying $\mathop{\mathrm{div}}\nolimits \boldsymbol{v}\in L^2(\Omega)$. We denote the corresponding norm $\left\| {{ \boldsymbol{v}}} \right\|^2_{H( \mathop{\mathrm{div}}\nolimits,\Omega)}: = \left\| {{ \boldsymbol{v}}} \right\|_{0,\Omega}^2+\left\| {{ \mathop{\mathrm{div}}\nolimits \boldsymbol{v}}} \right\|^2_{0,\Omega}$. Similarly, the space of tensors in $L^2(\Omega, \mathbb{S})$ with divergence in $L^2(\Omega, \mathbb{R}^d)$ is denoted $H( \mathop{\bf{div}}\nolimits, \Omega, \mathbb{S})$. We maintain the same notation $\left\| {{ \boldsymbol{\tau}}} \right\|^2_{H( \mathop{\bf{div}}\nolimits,\Omega)}: = \left\| {{ \boldsymbol{\tau}}} \right\|_{0,\Omega}^2+\left\| {{ \mathop{\bf{div}}\nolimits \boldsymbol{\tau}}} \right\|^2_{0,\Omega}$ for the corresponding norm. Let $\boldsymbol{n}$ be the outward unit normal vector to $\partial \Omega$. The Green formula

valid for $\boldsymbol{\tau} \in H( \mathop{\bf{div}}\nolimits, \Omega, \mathbb{S})$, can be used to extend the normal trace operator $\mathcal{C}^\infty(\overline \Omega, \mathbb{S})\ni \boldsymbol{\tau} \to ( \boldsymbol{\tau}|_{\partial \Omega}) \boldsymbol{n}$ to a linear continuous mapping $(\cdot|_{\partial \Omega}) \boldsymbol{n}:\, H( \mathop{\bf{div}}\nolimits, \Omega, \mathbb{S}) \to H^{-\frac{1}{2}}(\partial \Omega, \mathbb{R}^d)$, where $H^{-\frac{1}{2}}(\partial \Omega, \mathbb{R}^d)$ is the dual of $H^{\frac{1}{2}}(\partial \Omega, \mathbb{R}^d)$.

Throughout this paper, we shall use the letter $C$ to denote a generic positive constant independent of the mesh size $h$ and the physical parameters $\kappa$ and $\mu$, that may stand for different values at its different occurrences. Moreover, given any positive expressions $X$ and $Y$ depending on $h$, $\kappa$, and $\mu$, the notation $X \,\lesssim\, Y$ means that $X \,\le\, C\, Y$.

2.   The pure-stress formulation of the Brinkman problem

Let $\Omega\subset \mathbb{R}^d$ be a bounded and connected Lipschitz domain with boundary $\Gamma: = \partial \Omega$. Our purpose is to solve the Brinkman model

that describes the flow of a fluid with dynamic viscosity $\mu>0$, with velocity field $\boldsymbol{u}:\Omega \to \mathbb{R}^d$ and pressure $p:\Omega \to \mathbb{R}$, in a porous medium characterized by a permeability coefficient $\kappa\in L^\infty(\Omega)$. The volume force $\boldsymbol{f}\in L^2(\Omega, \mathbb{R}^d)$ is given and we assume that

We assume a no-slip boundary condition $\boldsymbol{u} = \boldsymbol{0}$ on a subset $\Gamma_D\subset\Gamma : = \partial \Omega$ of positive surface measure and impose the normal stress boundary condition $\boldsymbol{\sigma} \boldsymbol{n} = \boldsymbol{0}$ on its complement $\Gamma_N: = \Gamma \setminus \Gamma_D$, where $\boldsymbol{n}$ represents the exterior unit normal vector on $\Gamma$. In the case $\Gamma_D = \Gamma$ we impose the zero mean value restriction $\left( {{ p, 1}} \right) = 0$ on the pressure to enforce uniqueness.

We want to impose the stress tensor $\boldsymbol{\sigma}$ as a primary variable. To this end, we write the deviatoric part of (1a) and use equation (1b) to eliminate $p$ and $\boldsymbol{u}$, respectively, and end up with the boundary value problem

The no-slip boundary condition (2b) should be understood as integrating by parts in (2a) and then using again (1b). In the case of a subset $\Gamma_N$ with positive surface measure, the essential boundary condition (2c) requires the introduction of the closed subspace of $H( \mathop{\bf{div}}\nolimits, \Omega, \mathbb{S})$ given by

where $\left\langle {{\cdot, \cdot}} \right\rangle_\Gamma$ holds for the duality pairing between $H^{\frac{1}{2}}(\Gamma, \mathbb{R}^d)$ and $H^{-\frac{1}{2}}(\Gamma, \mathbb{R}^d)$. Hence, the energy space is given by

Testing (2a) with $\boldsymbol{\tau}\in X$ and integrating by parts yields

This leads us to propose the following variational formulation of the problem: Find $\boldsymbol{\sigma}\in X$ such that

where

with $\theta = 1$ if a no-slip boundary condition is imposed everywhere on $\Gamma$ (namely, if $\Gamma_D = \Gamma$) and $\theta = 0$ otherwise. We point out that, in the case $\theta = 1$, testing problem (3) with $\boldsymbol{\tau} = \boldsymbol{I}$ gives $\left( {{ \mathop{\mathrm{tr}}\nolimits( \boldsymbol{\sigma}),1}} \right) = 0$, which corresponds to the zero mean value restriction on the pressure. We are then opting for a variational insertion of this condition in order to free the energy space $X = H( \mathop{\bf{div}}\nolimits, \Omega, \mathbb{S})$ from this constraint.

Theorem 2.1. The variational problem $(3)$ admits a unique solution and there exists a constant $C>0$, depending $\Omega$ and $\kappa$, such that

Proof. Let us first notice that, by virtue of

the bilinear form defining the variational problem (3) is bounded:

for all $\boldsymbol{\sigma}, \boldsymbol{\tau} \in X$. Moreover, as a consequence of the Poincaré–Friedrichs inequalities (see, e.g., [7,Proposition 9.1.1])

and (see [14,Lemma 2.4])

the bilinear form is also coercive on $X$ and the well-posedness of (3) is a consequence of Lax–Milgram Lemma. Indeed, if $\Gamma_N$ has positive surface measure (which corresponds to $\theta = 0$) the coercivity of the bilinear form follows directly from (5). In the case $\Gamma_D = \Gamma$ (and $\theta = 1$), we can take advantage of the $L^2(\Omega, \mathbb{M})$-orthogonal decomposition $\boldsymbol{\tau} = \boldsymbol{\tau}_0 + \frac{1}{d|\Omega|} \left( {{ \mathop{\mathrm{tr}}\nolimits \boldsymbol{\tau}, 1}} \right) \boldsymbol{I}$ and the properties $\boldsymbol{\tau}^\mathtt{D} = ( \boldsymbol{\tau}_0)^\mathtt{D}$, $\mathop{\bf{div}}\nolimits \boldsymbol{\tau}_0 = \mathop{\bf{div}}\nolimits \boldsymbol{\tau}$ and $\left( {{ \mathop{\mathrm{tr}}\nolimits \boldsymbol{\tau}_0 , 1}} \right) = 0$ to deduce the coercivity from (4) as follows:

for all $\boldsymbol{\tau} \in H( \mathop{\bf{div}}\nolimits, \Omega, \mathbb{M})$, with $\beta = \max\left\{(1 + \frac{1}{\alpha})\frac{1}{\kappa_-}, \frac{2}{\alpha}, \frac{1}{d|\Omega|} \right\}$.

Remark 2.2. Once the stress tensor $\boldsymbol{\sigma}$ is known, the remaining variables can be recovered from

Moreover, we point out that testing (3) with $\boldsymbol{\tau} = \phi \boldsymbol{I}$, with $\phi:\Omega \to \mathbb{R}$ smooth and compactly supported in $\Omega$, we readily deduce the incompressibility condition $\mathop{\mathrm{div}}\nolimits \boldsymbol{u} = 0$ in $\Omega$.

3.   Auxiliary results concerning discretization

From now on, we assume that there exists a polygonal/polyhedral disjoint partition $\big\{\Omega_j,\ j = 1,\ldots,J\big\}$ of $\bar \Omega$ such that $\kappa|_{\Omega_j}: = \kappa_j$ constant, for all $j = 1,\ldots,J$.

We consider a sequence $\{\mathcal{T}_h\}_h$ of shape regular meshes that subdivide the domain $\bar \Omega$ into triangles/tetrahedra $K$ of diameter $h_K$. This means that there exists $\vartheta>0$ such that

where $\varrho_K$ stands for the the diameter of the largest ball that can be inscribed in $K$. The parameter $h: = \max_{K\in \mathcal{T}_h} \{h_K\}$ represents the mesh size of $\mathcal{T}_h$. We assume that $\mathcal{T}_h$ is aligned with the partition $\bar\Omega = \cup_{j = 1}^J \bar{\Omega}_j$ and that $\mathcal{T}_h(\Omega_j) : = \left\{ {{K\in \mathcal{T}_h;\ K\subset \Omega_j }} \right\}$ is a shape regular mesh of $\bar\Omega_j$ for all $j = 1,\cdots, J$ and all $h$.

For all $s\geq 0$, we consider the broken Sobolev space

corresponding to the partition $\bar\Omega = \cup_{j = 1}^J \bar{\Omega}_j$. Its vectorial and tensorial versions are denoted $H^s(\cup_j\Omega_j, \mathbb{R}^d)$ and $H^s(\cup_j\Omega_j, \mathbb{M})$, respectively. Likewise, the broken Sobolev space with respect to the subdivision of $\bar \Omega$ into $\mathcal{T}_h$ is

For each $\boldsymbol{v}: = \left\{ {{ \boldsymbol{v}_K}} \right\}\in H^s( \mathcal{T}_h, \mathbb{R}^d)$ and $\boldsymbol{\tau}: = \left\{ {{ \boldsymbol{\tau}_K}} \right\}\in H^s( \mathcal{T}_h, \mathbb{M})$ the components $\boldsymbol{v}_K$ and $\boldsymbol{\tau}_K$ represent the restrictions $\boldsymbol{v}|_K$ and $\boldsymbol{\tau}|_K$. When no confusion arises, the restrictions of these functions will be written without any subscript.

Hereafter, given an integer $m\geq 0$ and a domain $D\subset \mathbb{R}^d$, $\mathcal{P}_m(D)$ denotes the space of polynomials of degree at most $m$ on $D$. We introduce the space

of piecewise polynomial functions relative to $\mathcal{T}_h$. Moreover, $\mathcal{P}_m( \mathcal{T}_h,E)$ is the space of functions with values in $E$ and entries in $\mathcal{P}_m( \mathcal{T}_h)$, where $E$ is either $\mathbb{R}^d$, $\mathbb{M}$ or $\mathbb{S}$.

Let us introduce now notations related to DG approximations of $H(\text{div})$-type spaces. We say that a closed subset $F\subset \overline{\Omega}$ is an interior edge/face if $F$ has a positive $(d-1)$-dimensional measure and if there are distinct elements $K$ and $K'$ such that $F = \bar K\cap \bar K'$. A closed subset $F\subset \overline{\Omega}$ is a boundary edge/face if there exists $K\in \mathcal{T}_h$ such that $F$ is an edge/face of $K$ and $F = \bar K\cap \Gamma$. We consider the set $\mathcal{F}_h^0$ of interior edges/faces, the set $\mathcal{F}_h^\partial$ of boundary edges/faces and let $\mathcal{F}(K): = \left\{ {{F\in \mathcal{F}_h:\ F\subset \partial K}} \right\}$ be the set of edges/faces composing the boundary of $K$.

We assume that the boundary mesh $\mathcal{F}_h^\partial$ is compatible with the partition $\partial \Omega = \Gamma_D \cup \Gamma_N$ in the sense that, if

then $\Gamma_D = \cup_{F\in \mathcal{F}_h^D} F$ and $\Gamma_N = \cup_{F\in \mathcal{F}_h^N} F$. We denote

Obviously, in the case $\Gamma_D = \Gamma$ we have that $\mathcal{F}^*_h = \mathcal{F}^0_h$.

We will need the space given on the skeletons of the triangulations $\mathcal{T}_h$ by $L^2( \mathcal{F}^*_h): = \bigoplus_{F\in \mathcal{F}^*_h} L^2(F)$. Its vector valued version is denoted $L^2( \mathcal{F}^*_h, \mathbb{R}^d)$. Here again, the components $\boldsymbol{v}_F$ of $\boldsymbol{v} : = \left\{ {{ \boldsymbol{v}_F}} \right\}\in L^2( \mathcal{F}^*_h, \mathbb{R}^d)$ coincide with the restrictions $\boldsymbol{v}|_F$. We endow $L^2( \mathcal{F}^*_h, \mathbb{R}^d)$ with the inner product

and denote the corresponding norm $\left\| {{ \boldsymbol{v}}} \right\|^2_{0, \mathcal{F}^*_h}: = ( \boldsymbol{v}, \boldsymbol{v})_{ \mathcal{F}^*_h}$. From now on, $h_ \mathcal{F}\in L^2( \mathcal{F}^*_h)$ is the piecewise constant function defined by $h_ \mathcal{F}|_F : = h_F$ for all $F \in \mathcal{F}^*_h$ with $h_F$ denoting the diameter of edge/face $F$. By virtue of our hypotheses on $\kappa$ and on the triangulation $\mathcal{T}_h$, we may consider that $\kappa$ is an element of $\mathcal{P}_0( \mathcal{T}_h)$ and denote $\kappa_K: = \kappa|_K$ for all $K\in \mathcal{T}_h$. Moreover, we introduce $\gamma_ \mathcal{F}\in L^2( \mathcal{F}^*_h)$ defined by $\gamma_F : = \min\{ \kappa_K^{-1}, \kappa_{K'}^{-1} \}$ if $K\cap K' = F$ and $\gamma_F : = \kappa_K^{-1}$ if $F\cap K \in \mathcal{F}^N_h$.

Given $\boldsymbol{v}\in H^s( \mathcal{T}_h, \mathbb{R}^d)$ and $\boldsymbol{\tau}\in H^s( \mathcal{T}_h, \mathbb{M})$, with $s>\frac{1}{2}$, we define averages $\left\{\kern-1.ex\left\{ { \boldsymbol{v}} \right\}\kern-1.ex\right\}\in L^2( \mathcal{F}^*_h, \mathbb{R}^d)$ and jumps $\left[\!\left[ { { \boldsymbol{\tau}} } \right]\!\right]\in L^2( \mathcal{F}^*_h, \mathbb{R}^d)$ by

with the conventions

where $\boldsymbol{n}_K$ is the outward unit normal vector to $\partial K$.

For any $k\geq 1$, we let $\mathcal{X}^k(h) : = X + \mathcal{P}_{k}( \mathcal{T}_h, \mathbb{S})$. Given $\boldsymbol{\tau} \in \mathcal{P}_{k}( \mathcal{T}_h, \mathbb{S})$, we define $\mathop{\bf{div}}\nolimits_h \boldsymbol{\tau} \in L^2(\Omega, \mathbb{R}^d)$ by $\mathop{\bf{div}}\nolimits_h \boldsymbol{\tau}|_{K} : = \mathop{\bf{div}}\nolimits \boldsymbol{\tau}_K$ for all $K\in \mathcal{T}_h$ and endow $\mathcal{X}^k(h)$ with the norm

Under the condition $\mathop{\bf{div}}\nolimits_h \boldsymbol{\tau} \in H^s( \mathcal{T}_h, \mathbb{R}^d)$ ($s>\frac{1}{2}$), we also introduce

We end this section by recalling technical results needed for the convergence analysis of problem (3). We begin with the following well-known trace inequality, see for example [10].

Lemma 3.1. There exists a constant $C>0$ independent of $h$ such that

for all $v\in H^1(K)$ and all $K\in \mathcal{T}_h$.

The following discrete trace inequality will also be useful in our analysis.

Lemma 3.2. There exists a constant $C_{ \text{tr}}>0$ depending only on $d$ and on the shape regularity of $\{ \mathcal{T}_h\}_h$ such that

Proof. It is shown in [30] that, for any $K\in \mathcal{T}_h$ and $F\in \mathcal{F}(K)$, it holds

As a consequence of (12), there exists $C_{0}>0$ depending only on the shape-regularity constant $\vartheta$ and $d$ such that

By definition of $\gamma_ \mathcal{F}$, for any $\boldsymbol{v} \in \mathcal{P}_k( \mathcal{T}_h, \mathbb{R}^d)$,

and the result follows from (13) with $C_{ \text{tr}} = \sqrt{d+1}C_0$.

The Scott–Zhang like quasi-interpolation operator $\pi_h:\, L^2(\Omega) \to \mathcal{P}_{k}( \mathcal{T}_h)\cap H^1(\Omega)$, obtained in [12] by applying an $L^2$-orthogonal projection onto $\mathcal{P}_{k}( \mathcal{T}_h)$ followed by an averaging procedure with range in the space of continuous and piecewise $\mathcal{P}_{k}$ functions, will be especially useful in the forthcoming analysis. We recall in the next lemma the local approximation properties given by [12,Theorem 5.2]. Let us first introduce some notations. For any $K\in \mathcal{T}_h$, we introduce the subset of $\mathcal{T}_h$ defined by $\mathcal{T}_h^K : = \left\{ {{K'\in \mathcal{T}_h:\, K\cap K' \neq \emptyset}} \right\}$ and let $D_K = \text{interior}\left(\cup_{K'\in \mathcal{T}_h^K} K' \right)$.

Lemma 3.3. The quasi-interpolation operator $\pi_h$ is invariant in the space $\mathcal{P}_{k}( \mathcal{T}_h)\cap H^1(\Omega)$ and there exists a constant $C>0$ independent of $h$ such that

for all real numbers $0\leq r \leq k+1$, all natural numbers $0\leq m \leq \lfloor r\rfloor$, all $v\in H^r(D_K)$, and all $K\in \mathcal{T}_h$. Here $\lfloor r\rfloor$ stands for the largest integer less than or equal to $r$.

We point out that, as a consequence of (14) and the triangle inequality, it holds

for all natural number $0\leq m \leq k+1$, all $v\in H^m(D_K)$ and all $K\in \mathcal{T}_h$. Moreover, it is straightforward to deduce from (14) and the trace inequality (10) that

for all $2 \leq r \leq k+1$ ($k\geq 1$), all $v\in H^r(D_K)$ and all $K\in \mathcal{T}_h$.

We can deduce from (15) a global stability property for $\pi_h$ on $H^m(\cup_j\Omega_j)$, $0\leq m \leq k+1$, by taking advantage of the fact that the cardinal $\#( \mathcal{T}_h^K)$ of $\mathcal{T}_h^K$ is uniformly bounded for all $K\in \mathcal{T}_h$ and all $h$, as a consequence of the shape-regularity of the mesh sequence $\{ \mathcal{T}_h\}$. Indeed, given $v\in H^m(\cup_j\Omega_j)$, we let $\mathcal{T}_h^K(\Omega_j) : = \left\{ {{K'\in \mathcal{T}_h(\Omega_j):\, K\cap K' \neq \emptyset}} \right\}$ be the subset of elements in $\mathcal{T}_h^K$ that are contained in $\bar\Omega_j$ and denote $D_K^j : = \text{interior}\left(\cup_{K'\in \mathcal{T}_h^K(\Omega_j)} K' \right)$. It follows from (15) that

Summing over $K\in \mathcal{T}_h^K(\Omega_j)$ and using that $\#( \mathcal{T}_h^{K}(\Omega_j)) \leq \#( \mathcal{T}_h^{K}) \leq c$ for all $1\leq j \leq J$ and all $h$, we deduce that

In what follows, we use the same notation for the tensorial version $\pi_h:\, L^{2}(\Omega, \mathbb{S}) \to \mathcal{P}_{k}( \mathcal{T}_h, \mathbb{S})\cap H^1(\Omega, \mathbb{S})$ of the quasi-interpolation operator, which is obtained by applying the scalar operator componentwise. Namely, for any $\boldsymbol{\tau} \in L^2(\Omega, \mathbb{S})$, we let $(\pi_h \boldsymbol{\tau})_{ij} : = \pi_h(\tau_{ij})$, $1 \leq i, j\leq d$, which ensures that $\pi_h \boldsymbol{\tau}$ is a symmetric tensor when $\boldsymbol{\tau}$ is symmetric. As a consequence of (14) and (16) we have the following result.

Lemma 3.4. There exists a constant $C>0$ independent of $h$ and $\kappa$ such that

for all $\boldsymbol{\tau} \in H( \mathop{\bf{div}}\nolimits, \Omega, \mathbb{S})\cap H^{r+1}(\cup_j\Omega_j, \mathbb{M})$, $r\geq 1$, provided $h\leq \sqrt \kappa_+$.

Proof. Given $K\in \mathcal{T}_h(\Omega_j)$, $1\leq j \leq J$, we obtain from (14) that

as well as

For the last term in the left-hand side of (18), we notice that

and (16) yields

Summing (19), (20) and (21) over $K\in \mathcal{T}_h(\Omega_j)$ and then over $j = 1,\ldots, J$ and invoking the shape-regularity of the mesh sequence give the result.

Given $s>1/2$ and $m\geq 1$, the tensorial version of the canonical interpolation operator $\Pi^{\texttt{BDM}}_h: H^s(\cup_j\Omega_j, \mathbb{R}^d) \cap H( \mathop{\mathrm{div}}\nolimits,\Omega) \to \mathcal{P}_m( \mathcal{T}_h, \mathbb{R}^d)\cap H( \mathop{\mathrm{div}}\nolimits,\Omega)$ associated with the Brezzi–Douglas–Marini (BDM) mixed finite element satisfies the following classical error estimate, see [7, Proposition 2.5.4],

for all $\boldsymbol{v} \in H^s(\cup_j\Omega_j, \mathbb{R}^d) \cap H( \mathop{\mathrm{div}}\nolimits,\Omega)$, $s>1/2$. Moreover, we have the well-known commutativity property,

where $Q^{m-1}_h$ stands for the $L^2(\Omega)$-orthogonal projection onto $\mathcal{P}_{m-1}( \mathcal{T}_h)$. Applying $\Pi^{\texttt{BDM}}_h$ row-wise to matrices we obtain $\boldsymbol{\Pi}^{\texttt{BDM}}_h:\, H^s(\cup_j\Omega_j, \mathbb{M}) \cap H( \mathop{\mathrm{div}}\nolimits,\Omega, \mathbb{M}) \to \mathcal{P}_m( \mathcal{T}_h, \mathbb{M}) \cap H( \mathop{\mathrm{div}}\nolimits,\Omega, \mathbb{M})$. Obviously, this tensorial version of the Brezzi-Douglas-Marini interpolation operator also satisfies

and

for all $\boldsymbol{\tau} \in H^s(\cup_j\Omega_j, \mathbb{M}) \cap H( \mathop{\bf{div}}\nolimits,\Omega, \mathbb{M})$, $s>1/2$, where $\boldsymbol{Q}^{m-1}_h$ is the $L^2(\Omega, \mathbb{R}^d)$-orthogonal projection onto $\mathcal{P}_{m-1}( \mathcal{T}_h, \mathbb{R}^d)$.

4.   The mixed-DG method and its convergence analysis

We assume that $\boldsymbol{f}\in H^1(\cup_j \Omega_j, \mathbb{R}^d)$ and consider the following DG discretization of (3): Find $\boldsymbol{\sigma}_h\in \mathcal{P}_{k}( \mathcal{T}_h, \mathbb{S})$ such that

for all $\boldsymbol{\tau} \in \mathcal{P}_{k}( \mathcal{T}_h, \mathbb{S})$, where $\mathtt{a}>0$ is a large enough given parameter. In the case $\Gamma_D = \Gamma$, we notice that taking $\boldsymbol{\tau} = \boldsymbol{I}$ in (25) implies that $\left( {{ \mathop{\mathrm{tr}}\nolimits \boldsymbol{\sigma}_h,1}} \right) = 0$.

Remark 4.1. Should the boundary conditions be modified to be non-homogeneous

for sufficiently regular Dirichlet velocity $\boldsymbol{g}_D$ and sufficiently regular normal Cauchy stress $\boldsymbol{g}_N$, then (25) needs to be rewritten as follows: Find $\boldsymbol{\sigma}_h\in \mathcal{P}_{k}( \mathcal{T}_h, \mathbb{S})$ such that

for all $\boldsymbol{\tau} \in \mathcal{P}_{k}( \mathcal{T}_h, \mathbb{S})$.

Proposition 4.2. The linear systems of equations corresponding to $(25)$ are symmetric and positive definite, provided $\mathtt{a} \geq 2 C_{ \text{tr}}^2 + \frac{1}{2}$.

Proof. We first point out that the mapping $\boldsymbol{\tau}_h \mapsto {\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert { \boldsymbol{\tau}_h} \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert}$ defined in (9) is actually a norm on $\mathcal{P}_k( \mathcal{T}_h, \mathbb{S})$. Indeed, if $\boldsymbol{\tau}_h \in \mathcal{P}_k( \mathcal{T}_h, \mathbb{S})$ satisfies ${\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert { \boldsymbol{\tau}_h} \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert} = 0$, then $\boldsymbol{\tau}_h$ is H(div)-conforming since the jumps of its normal components vanish across all the internal faces $F\in \mathcal{F}_h^0$. Moreover, it holds $\boldsymbol{\tau}_h \boldsymbol{n} = {\bf 0}$ on $\Gamma_N$. Hence, $\boldsymbol{\tau}_h = {\bf 0}$ as a consequence of (5).

By virtue of the Cauchy–Schwarz inequality, Young's inequality together with the discrete trace inequality (11) it holds,

for all $\boldsymbol{\tau} \in \mathcal{P}_{k}( \mathcal{T}_h, \mathbb{S})$. It follows from (28) that

for all ${\bf 0} \neq \boldsymbol{\tau}_h\in \mathcal{P}_k( \mathcal{T}_h, \mathbb{S})$, if $\mathtt{a}\geq 2 C_{ \text{tr}}^2 + \frac{1}{2}$, and the result follows.

Proposition 4.3. The solution $\boldsymbol{\sigma}$ of $(3)$ satisfies the following consistency property

Proof. Taking into account that $\left( {{ \mathop{\mathrm{tr}}\nolimits \boldsymbol{\sigma}, 1}} \right) = 0$ in the case $\Gamma_D = \Gamma$, and testing (3) with a tensor $\boldsymbol{\tau}:\Omega \to \mathbb{S}$ whose entries are indefinitely differentiable and compactly supported in $\Omega$, we deduce that $\boldsymbol{u} : = \frac{\kappa}{{\mu}} ( \mathop{\bf{div}}\nolimits \boldsymbol{\sigma} + \boldsymbol{f})$ satisfies $\boldsymbol{\varepsilon}( \boldsymbol{u}) = \frac{1}{2{\mu}} \boldsymbol{\sigma}^\mathtt{D} \in L^2(\Omega, \mathbb{S})$. Moreover, applying a Green's formula in (3) yields $\boldsymbol{u} = {\bf 0}$ on $\Gamma_D$. Hence, by virtue of Korn's inequality, $\boldsymbol{u}\in H^1_D(\Omega, \mathbb{R}^d): = \left\{ {{ \boldsymbol{v} \in H^1(\Omega, \mathbb{R}^d):\ \boldsymbol{v}|_{\Gamma_D} = {\bf 0}}} \right\}$, which ensures that

Furthermore, an integration by parts on each element $K\in \mathcal{T}_h$ gives

and the result follows by substituting $\boldsymbol{u} : = {\frac{\kappa}{\mu}} ( \mathop{\bf{div}}\nolimits \boldsymbol{\sigma} + \boldsymbol{f})$ in the last expression.

Lemma 4.4. Let $\boldsymbol{\sigma}$ and $\boldsymbol{\sigma}_h$ be the solutions of $(3)$ and $(25)$, respectively. There exists $\mathtt{a}_0>0$ such that

whenever the stability parameter $\mathtt{a}$ of problem $(25)$ satisfies $\mathtt{a} \geq \mathtt{a}_0$. Here, the constants $C$ and $\mathtt{a}_0$ are independent of $h$, $\kappa$, and $\mu$.

Proof. We introduce $\boldsymbol{e}_h = \pi_h \boldsymbol{\sigma} - \boldsymbol{\sigma}_h \in \mathcal{P}_k( \mathcal{T}_h, \mathbb{S})$ and notice that, thanks to (30),

for all $\boldsymbol{\tau}\in \mathcal{P}_{k}( \mathcal{T}_h, \mathbb{S})$, where

Taking $\boldsymbol{\tau} = \boldsymbol{e}_h$ in (32) yields

Let us bound now each of terms of the right-hand side of (33) by means of the Cauchy–Schwarz and Young inequalities. For the first term, proceeding as in (28) gives

Next, we estimate $G( \boldsymbol{e}_h)$ in two steps: from the one hand,

and from the other hand,

Substituting (34)-(36) in (33) and rearranging terms we deduce that (31) is satisfied if the stability parameter $\mathtt{a}$ of problem (25) fulfils the condition

Remark 4.5. The unknown constant $C_{ \text{tr}}>0$ only depends on the shape-regularity of the mesh sequence (see Lemma 3.2). Hence, the stability condition (37) is independent of the mesh size $h$, the polynomial degree $k$, and the data $\boldsymbol{f}$, $\kappa$ and $\mu$ of the problem.

Theorem 4.6. Let $\boldsymbol{\sigma}$ and $\boldsymbol{\sigma}_h$ be the solutions of problems $(3)$ and $(25)$, respectively. If $\boldsymbol{\sigma}\in H^{r+1}(\cup_j\Omega_j, \mathbb{M})$, with $r\geq 1$, then

for all $\mathtt{a} \geq \mathtt{a}_0$, with $C>0$ independent of $h$, $\kappa$, and $\mu$.

Proof. The result is a direct combination of (31) and the interpolation error estimates provided by Lemma 3.4.

5.   Discrete post-processing of velocity and pressure

We recall that the velocity field and the pressure can be recovered from the stress tensor at the continuous level by

These same expressions permit us to reconstruct, with local and independent calculations on each element, a discrete velocity $\boldsymbol{u}_h$ and a discrete pressure $p_h$ that converge to their continuous counterpart at an optimal rate of convergence, as shown in the following results.

Corollary 5.1. Let $\boldsymbol{\sigma}$ and $\boldsymbol{\sigma}_h$ be the solutions of problems $(3)$ and $(25)$, respectively. We introduce $\boldsymbol{u}_h : = {\frac{\kappa}{\mu}} ( \mathop{\bf{div}}\nolimits_h \boldsymbol{\sigma}_h + Q^{k-1}_h \boldsymbol{f} )$. If $\boldsymbol{\sigma}\in H^{r+1}(\cup_j\Omega_j, \mathbb{M})$ and $\boldsymbol{f}\in H^{r}(\cup_j\Omega_j, \mathbb{R}^d)$ with $r\geq 1$, there exists a constant $C>0$ independent of $h$, $\kappa$, and $\mu$ such that,

if the stability parameter $\mathtt{a}$ of problem $(25)$ satisfies $\mathtt{a} \geq \mathtt{a}_0$.

Proof. The result follows immediately from (38) and from the fact that

Corollary 5.2. Let $\boldsymbol{\sigma}$ and $\boldsymbol{\sigma}_h$ be the solutions of problems $(3)$ and $(25)$, respectively. We consider $p_h : = -\tfrac{1}{d} \mathop{\mathrm{tr}}\nolimits \boldsymbol{\sigma}_h$. If $\boldsymbol{\sigma}\in H^{r+1}(\cup_j\Omega_j, \mathbb{M})$ with $r\geq 1$, there exists a constant $C>0$ independent of $h$ and $\mu$, but depending on $\kappa_+/\kappa_-$ such that,

if the stability parameter $\mathtt{a}$ of problem $(25)$ satisfies $\mathtt{a} \geq \mathtt{a}_0$.

Proof. We notice that $p - p_h = \tfrac{1}{d} \mathop{\mathrm{tr}}\nolimits( \boldsymbol{\sigma}_h - \boldsymbol{\sigma}) \in Q$ where

It follows from the well-known inf-sup condition (see for example [11,Lemma 53.9])

that

Using an elementwise integration by parts formula gives

for all $\boldsymbol{v} \in H^1_D(\Omega, \mathbb{R}^d)$. Hence, using the Cauchy–Schwarz inequality

and the trace inequality (10) to estimate the term $\left\| {{h_ \mathcal{F}^{\frac{1}{2}} \boldsymbol{v}}} \right\|_{0, \mathcal{F}^*_h}$, we deduce that

Plugging the forgoing estimate in (39) yields

and it follows from (38) that there exists a constant $C>0$ independent of $h$, but depending on $\kappa_+/\kappa_-$ such that

which gives the result.

We aim now to obtain a velocity field $\boldsymbol{u}_h^*\simeq \boldsymbol{u}$ that preserves exactly the incompressibility condition $\mathop{\mathrm{div}}\nolimits \boldsymbol{u}_h^* = 0$ in $\Omega$. The computational cost for such an enhancement is more demanding. For $k\geq 2$, it is achieved by solving the following elliptic problem in mixed form: Find $\boldsymbol{u}_h^*\in H( \mathop{\mathrm{div}}\nolimits,\Omega)\cap \mathcal{P}_{k-1}( \mathcal{T}_h, \mathbb{R}^d)$ and $\lambda_h\in \mathcal{P}_{k-2}( \mathcal{T}_h)$ such that

Clearly, $\mathop{\mathrm{div}}\nolimits \boldsymbol{u}_h^* = 0$ in $\Omega$. We need now to estimate the error $\boldsymbol{u} - \boldsymbol{u}_h^*$ in the $L^2(\Omega, \mathbb{R}^d)$-norm.

Lemma 5.3. Let $\boldsymbol{\sigma}$ and $\boldsymbol{\sigma}_h$ be the solutions of problems $(3)$ and $(25)$, respectively. If $\boldsymbol{\sigma}\in H^{r+1}(\cup_j\Omega_j, \mathbb{M})$ and $\boldsymbol{f}\in H^{r}(\cup_j\Omega_j, \mathbb{R}^d)$ with $r\geq 1$, there exists a constant $C>0$ independent of $h$, $\kappa$, and $\mu$ such that, for all $\mathtt{a} \geq \mathtt{a}_0$,

where $\boldsymbol{u}_h^*$ is obtained by solving the auxiliary problem $(40)$.

Proof. We begin by considering the following auxiliary problem: Find $\boldsymbol{w}\in H( \mathop{\mathrm{div}}\nolimits, \Omega)$ and $\lambda\in L^2(\Omega)$ such that

We denote by $V$ the kernel of the bilinear form $H( \mathop{\mathrm{div}}\nolimits,\Omega) \times L^2(\Omega) \ni ( \boldsymbol{v}, \eta)\mapsto \left( {{\eta, \mathop{\mathrm{div}}\nolimits \boldsymbol{v}}} \right)$, that is

The fact that $( \boldsymbol{v}, \boldsymbol{v}) = \left\| {{ \boldsymbol{v}}} \right\|_{0,\Omega}^2 = \left\| {{ \boldsymbol{v}}} \right\|^2_{H( \mathop{\mathrm{div}}\nolimits, \Omega)}$ for all $\boldsymbol{v} \in V$ and the well-known inf-sup condition (cf. [7])

permit us to apply the Babuška–Brezzi theory to deduce that problem (41) is well-posed. Moreover, it can easily be seen that its unique solution is given by $\boldsymbol{w} = \boldsymbol{u}$ and $\lambda = 0$.

We recall that the BDM-mixed finite element pair

satisfies the discrete inf-sup condition

Moreover, we notice that discrete kernel $V_h$ of the bilinear form

is a subspace of $V$ since

Hence, the Galerkin method based on the BDM-element and defined by: Find $\boldsymbol{w}_h\in H( \mathop{\mathrm{div}}\nolimits, \Omega)\cap \mathcal{P}_{k-1}( \mathcal{T}_h, \mathbb{R}^d)$ and $\lambda_h\in \mathcal{P}_{k-2}( \mathcal{T}_h)$ such that

is stable, convergent and we have the Céa estimate (recall that $\boldsymbol{w} = \boldsymbol{u}$ and $\lambda = 0$)

We notice now that problem (40) is none other than a discretization of problem (41) obtained from the Galerkin method (42) after replacing the right-hand side $\left( {{ \boldsymbol{u}, \boldsymbol{v}}} \right)$ by $\left( {{ \boldsymbol{u}_h, \boldsymbol{v}}} \right)$ for all $\boldsymbol{v}\in H( \mathop{\mathrm{div}}\nolimits, \Omega)\cap \mathcal{P}_{k-1}( \mathcal{T}_h, \mathbb{R}^d)$. Hence, by virtue of Strang's Lemma, it immediately holds

Finally, we deduce from (22) and Corollary 5.1 that there exists a constant $C>0$ independent of $h$, $\kappa$, and $\mu$ such that

and the result follows.

6.   $L^2$-error estimates for the stress

The analysis of this section is restricted to the case $\Gamma_D = \Gamma$, so that $\theta = 1$, $X = H( \mathop{\bf{div}}\nolimits,\Omega, \mathbb{S})$ and $\mathcal{F}_h^* = \mathcal{F}_h^0$. The deduction of error estimates for the stress in the $L^2$-norm relies on the construction of adequate approximations of the exact solution $\boldsymbol{\sigma}\in H( \mathop{\bf{div}}\nolimits, \Omega, \mathbb{S})$ and of the discrete solution $\boldsymbol{\sigma}_h\in \mathcal{P}_k( \mathcal{T}_h, \mathbb{S})$ in the space $H( \mathop{\bf{div}}\nolimits, \Omega, \mathbb{S})\cap \mathcal{P}_k( \mathcal{T}_h, \mathbb{S})$.

One of the tools that are needed to achieve this goal is the averaging operator $\mathcal{A}_h:\, \mathcal{P}_k( \mathcal{T}_h, \mathbb{M}) \to \mathcal{P}_k( \mathcal{T}_h, \mathbb{M})\cap H( \mathop{\bf{div}}\nolimits,\Omega, \mathbb{M})$ defined on each $K\in \mathcal{T}_h$ and for any $\boldsymbol{\tau}\in \mathcal{P}_k( \mathcal{T}_h, \mathbb{M})$, by the conditions

Lemma 6.1. The projector $\mathcal{A}_h:\, \mathcal{P}_k( \mathcal{T}_h, \mathbb{M}) \to \mathcal{P}_k( \mathcal{T}_h, \mathbb{M})\cap H( \mathop{\bf{div}}\nolimits,\Omega, \mathbb{M})$ is uniquely characterized by the conditions $\rm (43a)-(43c)$ and it satisfies

with $C>0$ independent of $h$.

Proof. See [24,Proposition 5.2].

A combination of (44) and (38) shows that, under a sufficient regularity assumption on the exact solution $\boldsymbol{\sigma}\in H( \mathop{\bf{div}}\nolimits, \Omega, \mathbb{S})$, the term $\mathcal{A}_h \boldsymbol{\sigma}_h- \boldsymbol{\sigma}_h$ converges to zero in the $L^2$-norm. However, it is clear that the operator $\mathcal{A}_h$ does not preserve symmetry. To remedy this drawback, we follow [13,16,26,29] and use a symmetrization procedure that requires the stability the Scott–Vogelius element [27] for the Stokes problem. We refer to [11,Section 55.3] for a detailed account on the conditions (on the mesh $\mathcal{T}_h$ and on $k$) under which this stability property is guaranteed in the 2D and 3D cases.

Lemma 6.2. Assume that the pair $\left\{ {{ \mathcal{P}_{k+1}( \mathcal{T}_h, \mathbb{R}^d)\cap H^1(\Omega, \mathbb{R}^d), \mathcal{P}_{k}( \mathcal{T}_h, \mathbb{R}^d)}} \right\}$ is stable for the Stokes problem on the mesh $\mathcal{T}_h$. Then, there exists a linear operator

such that, for all $\boldsymbol{\tau}_h \in \mathcal{P}_k( \mathcal{T}_h, \mathbb{M})\cap H( \mathop{\bf{div}}\nolimits,\Omega, \mathbb{M})$,

i) $\mathop{\bf{div}}\nolimits ( \boldsymbol{\tau}_h - \mathcal S_h \boldsymbol{\tau}_h) = {\bf 0}$ in $\Omega$,

ii) and $\left\| {{ \boldsymbol{\tau}_h - \mathcal S_h \boldsymbol{\tau}_h}} \right\|_{0,\Omega}\lesssim \left\| {{ \boldsymbol{\tau}_h - \boldsymbol{\tau}_h^{ \mathtt{t}}}} \right\|_{0,\Omega}.$

Proof. See for example [29,Lemma 5.2].

Remark 6.3. The reason we are limiting the analysis of this section to the case $\Gamma_D = \Gamma$ is due to the fact that, to the authors knowledge, the symmetrization procedure provided by operator $\mathcal S_h$ has only been addressed in the literature for homogeneous Dirichlet or Neumann boundary conditions.

The next result shows that, under sufficient regularity assumptions on the exact solution $\boldsymbol{\sigma}$ of problem (3), $\boldsymbol{\sigma}_h - \mathcal S_h(\mathcal{A}_h \boldsymbol{\sigma}_h)$ and $\boldsymbol{\sigma} - \mathcal S_h( \boldsymbol{\Pi}^{\texttt{BDM}}_h \boldsymbol{\sigma} )$ converge to zero at optimal order in the $L^2$-norm.

Corollary 6.4. Let $\boldsymbol{\sigma}$ and $\boldsymbol{\sigma}_h$ be the solutions of problems $(3)$ and $(25)$, respectively. Assume that the pair $\left\{ {{ \mathcal{P}_{k+1}(\Omega, \mathbb{R}^d)\cap H^1(\Omega, \mathbb{R}^d), \mathcal{P}_{k}(\Omega, \mathbb{R}^d)}} \right\}$ is Stokes-stable on the mesh $\mathcal{T}_h$. Then, if $\boldsymbol{\sigma}\in H^{r+1}(\cup_j\Omega_j, \mathbb{M})$, with $r\geq 1$, there exists $C>0$ independent of $h$ and $\mu$, but depending on $\kappa_+/\kappa_-$, such that

for all $\mathtt{a} \geq \mathtt{a}_0$.

Proof. We point out that, as a consequence of property ii) in Lemma 6.2 and the symmetry of $\boldsymbol{\sigma}$, we have that

Using this time property ii) of Lemma 6.2 in combination with the symmetry of $\boldsymbol{\sigma}_h$ and (44) give

The result follows now by using (23) in (46) and (38) in (47).

Theorem 6.5. Let $\boldsymbol{\sigma}$ and $\boldsymbol{\sigma}_h$ be the solutions of problems $(3)$ and $(25)$, respectively. Assume that the pair $\left\{ {{ \mathcal{P}_{k+1}(\Omega, \mathbb{R}^d)\cap H^1(\Omega, \mathbb{R}^d), \mathcal{P}_{k}(\Omega, \mathbb{R}^d)}} \right\}$ is Stokes-stable on the mesh $\mathcal{T}_h$. Then, if $\boldsymbol{\sigma}\in H^{r+1}(\cup_j\Omega_j, \mathbb{M})$, with $r\geq 1$, there exists $C>0$ independent of $h$ and $\mu$, but depending on $\kappa_+/\kappa_-$, such that

for all $\mathtt{a} \geq \mathtt{a}_0$.

Proof. It follows from the consistency property (30) that

for all $\boldsymbol{\tau}_h \in \mathcal{P}_{k}( \mathcal{T}_h, \mathbb{S})\cap H( \mathop{\bf{div}}\nolimits,\Omega, \mathbb{S})$. Using an integration by parts in each $K\in \mathcal{T}_h$ gives

and taking advantage of (43b) yields

for all $\boldsymbol{\tau}_h \in \mathcal{P}_{k}( \mathcal{T}_h, \mathbb{S})\cap H( \mathop{\bf{div}}\nolimits,\Omega, \mathbb{S})$. Next, by construction of $\mathcal{A}_h$ it holds

which means that (49) can be equivalently written

Substituting back the last identity in (48) yields

for all $\boldsymbol{\tau}_h \in \mathcal{P}_{k}( \mathcal{T}_h, \mathbb{S})\cap H( \mathop{\bf{div}}\nolimits,\Omega, \mathbb{S})$. On the other hand, we notice that by virtue of (24) and property i) of Theorem 6.2,

and hence, taking $\boldsymbol{\tau}_h = \mathcal S_h( \boldsymbol{\Pi}^{\texttt{BDM}}_h \boldsymbol{\sigma} - \mathcal{A}_h \boldsymbol{\sigma}_h) \in \mathcal{P}_{k}( \mathcal{T}_h, \mathbb{S})\cap H( \mathop{\bf{div}}\nolimits,\Omega, \mathbb{S})$ in (50) gives

Consequently,

and the Cauchy–Schwarz inequality implies that

The result follows from (45) and from the fact that $\left( {{ \mathop{\mathrm{tr}}\nolimits( \boldsymbol{\sigma} - \boldsymbol{\sigma}_h), 1}} \right) = 0$.

It only remains now to show that the convergence of $p_h: = -\tfrac{1}{d} \mathop{\mathrm{tr}}\nolimits \boldsymbol{\sigma}_h$ to $p: = -\tfrac{1}{d} \mathop{\mathrm{tr}}\nolimits \boldsymbol{\sigma}$ is also enhanced by one order in the $L^2$-norm, where $\boldsymbol{\sigma}$ and $\boldsymbol{\sigma}_h$ are the solutions of problems (3) and (25), respectively.

Corollary 6.6. Assume that $\left\{ {{ \mathcal{P}_{k+1}(\Omega, \mathbb{R}^d)\cap H^1(\Omega, \mathbb{R}^d), \mathcal{P}_{k}(\Omega, \mathbb{R}^d)}} \right\}$ is a stable Stokes pair on the mesh $\mathcal{T}_h$. If $\boldsymbol{\sigma}\in H^{r+1}(\cup_j\Omega_j, \mathbb{M})$, with $r\geq 1$, there exists $C>0$ independent of $h$ and $\mu$, but depending on $\kappa_+/\kappa_-$, such that

for all $\mathtt{a} \geq \mathtt{a}_0$.

Proof. Using the triangle inequality and (52) yields

which permits us to deduce from (51) that

Consequently, it follows from (53), (54) and the coercivity property (6) that

Finally, by virtue of the triangle inequality,

and the result follows from (45).

7.   Numerical results

We now include a set of numerical examples that illustrate the convergence properties of the proposed DG method, and the usability of the formulation in simulating typical viscous flow in porous media. The stabilisation constant $\mathtt{a}$ is experimentally adjusted for each example so that it satisfies (37). All computational tests were conducted using the open source finite element library FEniCS [2].

Convergence tests. First we compare approximate and closed-form exact solutions for various levels of uniform mesh refinement. Let us consider the unit square domain $\Omega = (0,1)^2$, the parameters $\mu = 10^{-3}$, $\kappa = 1$, $\mathtt{a} = 10$, and the following manufactured smooth velocity and pressure

from which the exact Cauchy stress $\boldsymbol{\sigma}_{\mathrm{manuf}}$, forcing term, imposed non-homogeneous velocity $\boldsymbol{g}_D$ on $\Gamma_D = \{ 0 \}\times (0,1) \cup (0,1)\times\{1\}$ (left and top sides of the square), and imposed non-homogeneous normal stress $\boldsymbol{g}_N$ on $\Gamma_N = \partial\Omega\setminus\Gamma_D$ are constructed. The boundary partition indicates that we use the parameter $\theta = 0$.

Errors between exact and approximate solutions computed using (27) are measured in the following norms

and the experimental rates of convergence are computed as

where ${\tt e},\tilde{{\tt e}}$ denote errors generated on two consecutive meshes of sizes $h$ and $\tilde{h}$, respectively. Such an error history is displayed, for uniform meshes composed of squares split along one diagonal, in Table 1.

As anticipated by Theorem 4.6, a $k$-order of convergence is observed for the Cauchy stress in the norm defined in (9). An agreement with the error estimates from Corollaries 5.1 and 5.2 is also observed for the post-processed velocity and pressure in the $L^2$-norms. The error bounds from Lemma 5.3 are also confirmed for the discrete velocity field $\boldsymbol{u}_h^*$, which is obtained by solving problem (40). We show the decay of the error in the $H( \mathop{\mathrm{div}}\nolimits)$-norm, but we recall that the divergence of $\boldsymbol{u}_h^*$ is zero to machine precision. Note that, for $k = 1$, this second velocity post-processing requires the solution of (40) with the lowest order BDM-mixed pair $\left\{ {{H( \mathop{\mathrm{div}}\nolimits, \Omega)\cap \mathcal{P}_{1}( \mathcal{T}_h, \mathbb{R}^2), \mathcal{P}_{0}( \mathcal{T}_h)}} \right\}$.

We point out that, for the type of meshes considered in Table 1, the additional order of convergence for the deviatoric stress and for the pressure in the $L^2$-norm are not achieved. However, if we consider instead special partitions using, for example, simplicial barycentric trisected – Hsieh–Clough–Tocher meshes – or twice quadrisected crisscrossed meshes, the additional order shown in Theorem 6.5 and Corollary 6.6 is clearly obtained, as seen in Table 2, because in this case the Scott–Vogelius pair $\left\{ {{ \mathcal{P}_{k+1}( \mathcal{T}_h, \mathbb{R}^d)\cap H^1(\Omega, \mathbb{R}^d), \mathcal{P}_{k}( \mathcal{T}_h, \mathbb{R}^d)}} \right\}$ is Stokes inf-sup stable [11,Section 55.3.1].

We remark that the limit case of $\kappa = 0$ cannot be studied with the present formulation. Moreover, while for the other limit case of $\mu = 0$ the method converges optimally in all the stress norms as well as in the pressure post-processing, we cannot use (7) to recover velocity (and the second post-process depends on the first one). Nevertheless, we point out that the method still converges optimally for stress and all post-processed fields, for very low values of viscosity and/or permeability, and for very high values of permeability (confirmed through tests over several orders of magnitude).

2D maze and channel flow with heterogeneous permeability. Now we focus on two different geometries. For the first 2D example we consider a maze-shaped geometry of length 2.2 and height 1 (in adimensional units). We use an unstructured simplicial barycentric trisected grid (which guaranties the stability the Scott–Vogelius element $\left\{ {{ \mathcal{P}_{k+1}( \mathcal{T}_h, \mathbb{R}^d)\cap H^1(\Omega, \mathbb{R}^d), \mathcal{P}_{k}( \mathcal{T}_h, \mathbb{R}^d)}} \right\}$ for the Stokes problem [11,Section 55.3.1] and ensures the results of Section 6) with 66006 elements, representing, for $k = 1$, a total of 594054 degrees of freedom. The rightmost segment of the boundary is the outlet ($\Gamma_N$), where the outflow condition $\boldsymbol{\sigma} \boldsymbol{n} = \boldsymbol{0}$ is imposed. The remainder of the boundary is $\Gamma_D$, split between the leftmost segment of the boundary (the inlet, where we impose the parabolic profile $\boldsymbol{g}_{\mathrm{in}} = (100(y-0.45)(0.55-y),0)^{\tt t}$) and the walls where $\boldsymbol{g}_{\mathrm{wall}} = \boldsymbol{0}$. The external force is zero, the viscosity is $\mu = 10^{-6}$, and the permeability $\kappa( \boldsymbol{x})$ {is characterized by a non-homogeneous field taking the value $\kappa_{+} = 10^{-5}$ everywhere on the domain, except on 60 small disks distributed randomly, where the permeability smoothly goes down to the much smaller value $\kappa_{-} = 10^{-10}$:

where $(q_x(i),q_y(i))$ denote the coordinates of the randomly located points}. The stabilization parameter is $\mathtt{a} = 15$. Figure 1 displays the permeability field, the Cauchy stress magnitude, the line integral convolutions of post-processed velocity, and the distribution of post-processed pressure, where we observe that the expected symmetry of the flow is disrupted by the heterogeneous permeability. These results indicate well-resolved approximations.

In the second 2D example we follow [20] and compute channel flow solutions and take the permeability distribution from the Cartesian SPE10 benchmark dataset for reservoir simulations / model 2 (we choose layer 45, corresponding to a fluvial fan pattern with channeling. See also, e.g., [1,9]). The permeability data has a very large contrast: the minimum and maximum values are $\kappa_{-} = 1.3\cdot10^{-18}\,[\text{m}^2]$ and $\kappa_{+} = 2.0\cdot10^{-11}\,[\text{m}^2]$, and it is projected onto a twice quadrisected crisscrossed mesh discretizing the rectangular domain $\Omega = (0,6.096)\,[\text{m}]\times(0,3.048)$ [m]. The remaining parameters are $\mathtt{a} = 25$ and $\mu = 10^{-6}\,[\text{KPa}\cdot\text{s}]$. On the inlet (the left segment of the boundary) we impose the inflow velocity $\boldsymbol{u} = (10,0)^{\tt t}\,$[m/s], on the outlet (the right vertical segment) we set zero normal stress, and on the horizontal walls we impose no-slip conditions.

The flow patterns are shown in Figure 2, and the qualitative behaviour coincides with the expected filtration mechanisms observed elsewhere (see, e.g., [22]). The magnitude of velocity in its second post-process $| \boldsymbol{u}^*_h|$, and the stress magnitude are plotted in logarithmic scale.

3D flow through porous media. To close this section, in this test we use a cylindrical geometry of radius 0.2 [m] and height 0.7 [m], discretized into an unstructured simplicial barycentric quadrisected mesh of 48756 tetrahedra. We set the viscosity to the value of water $\mu = 8.9\cdot 10^{-4}\,[\mathrm{Pa}\cdot\mathrm{s}]$ and again use inlet (the face $z = 0$), wall (the surface $r = 0.2$ [m]), and outlet (the face $z = 0.7$ [m]) type of boundary conditions mimicking a channel flow with homogeneous right-hand side in the momentum balance $\boldsymbol{f} = \boldsymbol{0}$. The inlet velocity is $\boldsymbol{g}_{\mathrm{in}} = (0,0,1)^{\tt t}\,[\mathrm{m/s}]$. A synthetic permeability field is generated based on a porosity of approximately 0.4. For this we generate a normalized uniform random distribution, apply a Gaussian filter with standard deviation of 2, and then use a quantile of 40%. The resulting field is rescaled between $\kappa_{-} = 10^{-12}\,[\text{m}^2]$ and $\kappa_{+} = 10^{-4}\,[\text{m}^2]$. Figure 3 presents the permeability and contour iso-surfaces of velocity magnitude, which indicate the formation of wormhole-like structures with higher velocity.

Acknowledgments

This research was supported by Spain's Ministry of Economy Project PID2020-116287GB-I00, by the Monash Mathematics Research Fund S05802-3951284, by the Ministry of Science and Higher Education of the Russian Federation within the framework of state support for the creation and development of World-Class Research Centres "Digital biodesign and personalized healthcare" No. 075-15-2022-304, and by the Australian Research Council through the Discovery Project grant DP220103160.