An H^1-conforming Virtual Element Method for Darcy equations and Brinkman equations

TL;DR

This paper introduces an H^1-conforming Virtual Element Method for Darcy and Brinkman equations, achieving optimal convergence and stability across different flow regimes with rigorous analysis and numerical validation.

Contribution

It develops a new Virtual Element space with computable projections and divergence-free properties, applicable to both Darcy and Brinkman equations, ensuring stability and optimal convergence.

Findings

Optimal order convergence for Darcy equations

Stable and robust method for Brinkman equations

Numerical tests confirm theoretical error estimates

Abstract

The focus of the present paper is on developing a Virtual Element Method for Darcy and Brinkman equations. In [15] we presented a family of Virtual Elements for Stokes equations and we defined a new Virtual Element space of velocities such that the associated discrete kernel is pointwise divergence-free. We use a slightly different Virtual Element space having two fundamental properties: the L^2-projection onto P_k is exactly computable on the basis of the degrees of freedom, and the associated discrete kernel is still pointwise divergence-free. The resulting numerical scheme for the Darcy equation has optimal order of convergence and H^1 conforming velocity solution. We can apply the same approach to develop a robust virtual element method for the Brinkman equation that is stable for both the Stokes and Darcy limit case. We provide a rigorous error analysis of the method and several numerical tests.

Figures

Tables

Table 1. Table 1: Test 6.2 : Error for the velocities and the pressures for the “div-free” method.

	DoFs	$\mathrm{error}(𝐮,H^1)$	$\mathrm{error}(𝐮,L^2)$	$\mathrm{error}(p,L^2)$
$\mu =1e-01$	$182$	$2.049871825e-01$	$9.414645391e-03$	$1.531569296e-02$
	$702$	$4.616835760e-02$	$8.379142208e-04$	$2.796060666e-03$
	$2794$	$1.102679000e-02$	$9.416547836e-05$	$5.322283997e-04$
	$11210$	$2.654465229e-03$	$1.104272204e-05$	$1.261317758e-04$
$\mu =1e-04$	$182$	$2.563406238e-01$	$1.296095515e-02$	$6.431351247e-03$
	$702$	$5.462263791e-02$	$1.090600946e-03$	$1.887150783e-03$
	$2794$	$1.246452741e-02$	$1.179870900e-04$	$4.203480846e-04$
	$11210$	$2.790603844e-03$	$1.238676976e-05$	$1.026912579e-04$
$\mu =1e-14$	$182$	$2.572957705e-01$	$1.301886694e-02$	$6.431351247e-03$
	$702$	$5.539413175e-02$	$1.111681710e-03$	$1.887150783e-03$
	$2794$	$1.299961549e-02$	$1.253090639e-04$	$4.203480846e-04$
	$11210$	$3.003059376e-03$	$1.394861018e-05$	$1.026912579e-04$

Table 2. Table 2: Test 6.2 : Error for the velocities and the pressures for the “non div-free” method.

	DoFs	$\mathrm{error}(𝐮,H^1)$	$\mathrm{error}(𝐮,L^2)$	$\mathrm{error}(p,L^2)$
$\mu =1e-01$	$246$	$2.074920846e-01$	$9.923688300e-03$	$1.597079423e-02$
	$958$	$4.660732867e-02$	$8.775839313e-04$	$2.901160116e-03$
	$3818$	$1.113343018e-02$	$1.007780759e-04$	$5.396527876e-04$
	$15306$	$2.677875898e-03$	$1.189832113e-05$	$1.265284965e-04$
$\mu =1e-04$	$246$	$2.477687988e-01$	$1.163768487e-02$	$6.584759134e-03$
	$958$	$8.730792836e-02$	$1.749981712e-03$	$1.922719665e-03$
	$3818$	$5.060500911e-02$	$5.351452940e-04$	$4.258968745e-04$
	$15306$	$3.007229784e-02$	$1.597800852e-04$	$1.036140919e-04$
$\mu =1e-14$	$246$	$2.485968435e-01$	$1.168014371e-02$	$6.581632957e-03$
	$958$	$9.149517231e-02$	$1.837032059e-03$	$1.922719665e-03$
	$3818$	$6.422656370e-02$	$6.779428128e-04$	$4.260009082e-04$
	$15306$	$6.289228510e-02$	$3.325972217e-04$	$1.036819751e-04$

Equations

$$\{\begin{array}{rlrl}& \text{find (u,p) such that}& & \\ & {\mathbb{K}}^{-1}\mathbf{u}+\nabla p=0& & \text{in Ω,}\\ & \mathrm{div}\mathbf{u}=f& & \text{in Ω,}\\ & \mathbf{u}\cdot \mathbf{n}=0& & \text{on ∂Ω,}\end{array}$$

$$\mathbf{V}:=\left\{\mathbf{u}\in H(\mathrm{div},\Omega ),\text{s.t}\mathbf{u}\cdot \mathbf{n}=0\text{on ∂Ω}\right\},Q:=L_0^2(\Omega )=\left\{q\in L^2(\Omega )\text{s.t.}{\int }_{\Omega }q\mathrm{d}\Omega =0\right\}$$

$$\Vert \mathbf{v}{\Vert }_{\mathbf{V}}^2:=\Vert \mathbf{v}{\Vert }_{{\left[L^2(\Omega )\right]}^2}^2+\Vert \mathrm{div}\mathbf{v}{\Vert }_{L^2(\Omega )}^2,\Vert q{\Vert }_Q:=\Vert q{\Vert }_{L^2(\Omega )},$$

$$a(\mathbf{u},\mathbf{v}):={\int }_{\Omega }{\mathbb{K}}^{-1}\mathbf{u}\cdot \mathbf{v}\mathrm{d}\Omega ,\text{for all u,v∈V}$$

$$b(\mathbf{v},q):={\int }_{\Omega }\mathrm{div}\mathbf{v}q\mathrm{d}\Omega \text{for all v∈V, q∈Q.}$$

$$\{\begin{array}{rlr}& \text{find (u,p)∈V×Q, such that}& \\ & a(\mathbf{u},\mathbf{v})+b(\mathbf{v},p)=0& \text{for all v∈V,}\\ & b(\mathbf{u},q)=(f,q)& \text{for all q∈Q,}\end{array}$$

$$(f,q):={\int }_{\Omega }fq\mathrm{d}\Omega \text{for all q∈Q.}$$

$$\mathbf{Z}:=\{\mathbf{v}\in \mathbf{V}\text{s.t.}b(\mathbf{v},q)=0\text{for all q∈Q}\};$$

$$\Vert \mathbf{v}{\Vert }_{\mathbf{V}}:=\Vert \mathbf{v}{\Vert }_{{\left[L^2(\Omega )\right]}^2}\text{for all v∈Z.}$$

$$|a(\mathbf{u},\mathbf{v})|\le \Vert a\Vert \Vert \mathbf{u}{\Vert }_{\mathbf{V}}\Vert \mathbf{v}{\Vert }_{\mathbf{V}}\text{for all u,v∈V,}$$

$$|b(\mathbf{v},q)|\le \Vert b\Vert \Vert \mathbf{v}{\Vert }_{\mathbf{V}}\Vert q{\Vert }_Q\text{for all v∈V and q∈Q;}$$

$$a(\mathbf{v},\mathbf{v})\ge \alpha \Vert \mathbf{v}{\Vert }_{\mathbf{V}}^2\text{for all v∈Z;}$$

$$\exists \beta >0\text{such that}\underset{\mathbf{v}\in \mathbf{V}\mathbf{v}\ne 0}{\sup }\frac{b(\mathbf{u},q)}{\Vert \mathbf{v}{\Vert }_{\mathbf{V}}}\ge \beta \Vert q{\Vert }_Q\text{for all q∈Q.}$$

$$\Vert \mathbf{u}{\Vert }_{\mathbf{V}}+\Vert p{\Vert }_Q\le C\Vert f{\Vert }_{L^2(\Omega )}$$

$$h_K:=\mathrm{diameter}(K),h:=\underset{K\in {\mathcal{T}}_h}{\sup }{h}_K.$$

$$\mathbf{W}_{h}^{K}:=\biggl{\{}\mathbf{v}\in[H^{1}(K)]^{2}\quad\text{s.t}\quad\mathbf{v}_{|{\partial K}}\in[\mathbb{B}_{k}(\partial K)]^{2}\,,\biggr{.}\\ \left.\biggl{\{}\begin{aligned} &-\boldsymbol{\Delta}\mathbf{v}-\nabla s\in\mathcal{G}_{k-2}(K)^{\perp},\\ &{\rm div}\,\mathbf{v}\in\mathbb{P}_{k-1}(K),\end{aligned}\biggr{.}\qquad\text{ for some $s\in L^{2}(K)$}\quad\right\}$$

$$\begin{array}{rl}\dim \left({\mathbf{W}}_h^K\right)& =\dim \left([{\mathbb{B}}_k(\partial K){]}^2\right)+\dim \left({\mathcal{G}}_{k-2}(K{)}^{\perp }\right)+\left(\dim ({\mathbb{P}}_{k-1}(K))-1\right)\\ & =2n_{K}k+\frac{(k-1)(k-2)}{2}+\frac{(k+1)k}{2}-1.\end{array}$$

$${\int }_K\mathbf{v}\cdot {\mathbf{g}}_{k-2}^{\perp }\mathrm{d}K\text{for all gk−2⊥∈Gk−2(K)⊥,}$$

$${\int }_K(\mathrm{div}\mathbf{v})q_{k-1}\mathrm{d}K\text{for all qk−1∈Pk−1(K)/R.}$$

$$\{\begin{array}{rl}& {\int }_K\nabla {\mathbf{q}}_k:\nabla ({\mathbf{v}}_h-{\Pi }_k^{\nabla ,K}{\mathbf{v}}_h)\mathrm{d}K=0\text{for all qk∈[Pk(K)]2,}\\ & {\Pi }_0^{0,K}({\mathbf{v}}_h-{\Pi }_k^{\nabla ,K}{\mathbf{v}}_h)=0,\end{array}$$

$${\int }_K{\mathbf{q}}_n\cdot ({\mathbf{v}}_h-{\Pi }_n^{0,K}{\mathbf{v}}_h)\mathrm{d}K=0\text{for all qn∈[Pn(K)]2.}$$

$$\mathbf{U}_{h}^{K}:=\biggl{\{}\mathbf{v}\in[H^{1}(K)]^{2}\quad\text{s.t}\quad\mathbf{v}_{|{\partial K}}\in[\mathbb{B}_{k}(\partial K)]^{2}\,,\biggr{.}\\ \left.\biggl{\{}\begin{aligned} &-\boldsymbol{\Delta}\mathbf{v}-\nabla s\in\mathcal{G}_{k}(K)^{\perp},\\ &{\rm div}\,\mathbf{v}\in\mathbb{P}_{k-1}(K),\end{aligned}\biggr{.}\qquad\text{ for some $s\in L^{2}(K)$}\quad\right\}$$

$${\mathbf{V}}_h^K:=\left\{\mathbf{v}\in {\mathbf{U}}_h^K\text{s.t.}{\left(\mathbf{v}-{\Pi }_k^{\nabla ,K}\mathbf{v},{\mathbf{g}}_k^{\perp }\right)}_{[L^2(K){]}^2}=0\text{for all gk⊥∈Gk(K)⊥/Gk−2(K)⊥}\right\},$$

$$\dim \left([{\mathbb{B}}_k(\partial K){]}^2\right)=2n_{K}k,\dim \left({\mathbb{P}}_{k-1}(K)\right)=\frac{k(k+1)}{2},\dim \left({\mathcal{G}}_k(K{)}^{\perp }\right)=\frac{k(k+1)}{2}$$

$$\begin{array}{rl}\dim \left({\mathbf{U}}_h^K\right)& =2n_{K}k+\frac{k(k+1)}{2}+\frac{(k+1)k}{2}-1.\end{array}$$

$${\mathbf{D}}_{\mathbf{U}}:{\int }_K\mathbf{v}\cdot {\mathbf{g}}_k^{\perp }\mathrm{d}K\text{for all gk⊥∈Gk(K)⊥/Gk−2(K)⊥.}$$

$${\int }_{K}g\mathrm{d}\Omega ={\int }_{\partial K}{\mathbf{g}}_b\cdot \mathbf{n}\mathrm{d}s,$$

$${\mathbf{v}}_{|\partial K}={\mathbf{g}}_b,\mathrm{div}\mathbf{v}=g,-\mathbf{\Delta }\mathbf{v}-\nabla s=\mathbf{h}.$$

$$\dim \left({\mathbf{U}}_h^K\right)=\dim \left([{\mathbb{B}}_k(\partial K){]}^2\right)+\dim \left({\mathcal{G}}_k(K{)}^{\perp }\right)+\left(\dim ({\mathbb{P}}_{k-1}(K))-1\right)$$

$$\begin{array}{rl}\dim \left({\mathbf{V}}_h^K\right)& =2n_{K}k+\frac{(k-1)(k-2)}{2}+\frac{(k+1)k}{2}-1.\end{array}$$

Full text

1 Introduction

The Virtual Element Methods (in short, VEM or VEMs) is a recent technique for solving PDEs. VEMs were recently introduced in [5] as a generalization of the finite element method on polyhedral or polygonal meshes. In the numerical analysis and engineering literature there has been a recent growth of interest in developing numerical methods that can make use of general polygonal and polyhedral meshes, as opposed to more standard triangular/quadrilateral (tetrahedral/hexahedral) grids. Indeed, making use of polygonal meshes brings forth a range of advantages, including for instance automatic hanging node treatment, more efficient approximation of geometric data features, better domain meshing capabilities, more efficient and easier adaptivity, more robustness to mesh deformation, and others. This interest in the literature is also reflected in commercial codes, such as CD-Adapco, that have recently included polytopal meshes.

We refer to the recent papers and monographs [25, 11, 20, 21, 37, 42, 44, 43, 45, 50, 51, 29, 30, 41, 28] as a brief representative sample of the increasing list of technologies that make use of polygonal/polyhedral meshes. We mention here in particular the polygonal finite elements, that generalize finite elements to polygons/polyhedrons by making use of generalized non-polynomial shape functions, and the mimetic discretisation schemes [38, 12], that combine ideas from the finite difference and finite element methods.

The principal idea behind VEM is to use approximated discrete bilinear forms that require only integration of polynomials on the (polytopal) element in order to be computed. The resulting discrete solution is conforming and the accuracy granted by such discrete bilinear forms turns out to be sufficient to achieve the correct order of convergence. Following this approach, VEM is able to make use of very general polygonal/polyhedral meshes without the need to integrate complex non-polynomial functions on the elements and without loss of accuracy. Moreover, VEM is not restricted to low order converge and can be easily applied to three dimensions and use non convex (even non simply connected) elements. The Virtual Element Method has been developed successfully for a large range of problems, see for instance [5, 26, 6, 1, 17, 23, 3, 32, 19, 18, 40, 10, 13, 27, 52, 54, 48, 47, 4, 33, 31]. A helpful paper for the implementation of the method is [7].

The focus of this paper is on developing a new Virtual Element Method for the Darcy equation that is suitable for a robust extension to the (more complex) Brinkman problem. For such a problem, other VEM numerical schemes have been proposed, see for example [23, 8].

In [15] the authors developed a new Virtual Element Method for Stokes problems by exploiting the flexibility of the Virtual Element construction in a new way. In particular, they define a new Virtual Element space of velocities carefully designed to solve the Stokes problem. In connection with a suitable pressure space, the new Virtual Element space leads to an exactly divergence-free discrete velocity, a favorable property when more complex problems, such as the Navier-Stokes problem, are considered. We highlight that this feature is not shared by the method defined in [6] or by most of the standard mixed Finite Element methods, where the divergence-free constraint is imposed only in a weak (relaxed) sense.

In the present contribution we develop the Virtual Element Method for Darcy equations by introducing a slightly different virtual space for the velocities such that the local $L^{2}$ orthogonal projection onto the space of polynomials of degree less or equal than $k$ (where $k$ is the polynomial degree of accuracy of the method) can be computed using the local degrees of freedom. The resulting Virtual Elements family inherits the advantages on the scheme proposed in [15], in particular it yields an exactly divergence-free discrete kernel. Thus we obtain a stable Darcy element that is also uniformly stable for the Stokes problem. A sample of uniformly stable methods for Darcy-Stokes model is for instance [39, 53, 36, 49].

The last part of the paper deals with the analysis of a new mixed finite element method for Brinkman equations that stems from the above scheme for the Darcy problem. Mathematically, the Brinkman problem resembles both the Stokes problem for fluid flow and the Darcy problem for flow in porous media (see [35, 2, 34]). Constructing finite element methods to solve the Brinkman equation that are robust for both (Stokes and Darcy) limits is challenging. We will see how the above Virtual Element approach offers a natural and straightforward framework for constructing stable numerical algorithms for the Brinkman equations.

We remark that the proposed scheme belongs to the class of the pressure-robust method, i.e. delivers a velocity error independent of the continuous pressure.

The paper is organized as follows. In Section 2 we introduce the model continuous Darcy problem. In Section 3 we present its VEM discretisation. In Section 4 we detail the theoretical features and the convergence analysis of the problem. In Section 5 we develop a stable numerical methods for Brinkman equations. In Section 6 we show the numerical tests. Finally in the Appendix we present the theoretical analysis of the extension to the Darcy equation of the scheme of [6]. Even though this latter method is not recommended for the Darcy problem, the numerical experiments showed an unexpected optimal convergence rate for the pressure. We theoretically prove this behaviour, developing an inverse inequality for the VEM space, which is interesting on its own.

2 The continuous problem

We consider the classical Darcy equation that describes the flow of a fluid through a porous medium. Let $\Omega\subseteq\mathbb{R}^{2}$ be a bounded polygon then the Darcy equation in mixed form is

[TABLE]

where $\mathbf{u}$ and $p$ are respectively the velocity and the pressure fields, $f\in L^{2}(\Omega)$ is the source term and $\mathbb{K}$ is a uniformly symmetric, positive definite tensor that represents the permeability of the medium. From (1), since we have assumed no flux boundary conditions all over $\partial\Omega$, the external force $f$ has zero mean value on $\Omega$. We consider the spaces

[TABLE]

equipped with the natural norms

[TABLE]

and the bilinear forms $a(\cdot,\cdot)\colon\mathbf{V}\times\mathbf{V}\to\mathbb{R}$ and $b(\cdot,\cdot)\colon\mathbf{V}\times Q\to\mathbb{R}$ defined by:

[TABLE]

[TABLE]

Then the variational formulation of Problem (1) is

[TABLE]

where

[TABLE]

Let us introduce the kernel

[TABLE]

then it is straightforward to see that

[TABLE]

It is well known that (see for instance [24]):

• •

$a(\cdot,\cdot)$ and $b(\cdot,\cdot)$ are continuous, i.e.

[TABLE]

[TABLE]

• •

$a(\cdot,\cdot)$ is coercive on the kernel $\mathbf{Z}$, i.e. there exists a positive constant $\alpha$ depending on $\mathbb{K}$ such that

[TABLE]

• •

$b(\cdot,\cdot)$ satisfies the inf-sup condition, i.e.

[TABLE]

Therefore, Problem (4) has a unique solution $(\mathbf{u},p)\in\mathbf{V}\times Q$ such that

[TABLE]

with the constant $C$ depending only on $\Omega$ and $\mathbb{K}$.

3 Virtual formulation for Darcy equations

3.1 Decomposition and the original virtual element spaces

We outline the Virtual Element discretization of Problem (4). Here and in the rest of the paper the symbol $C$ will indicate a generic positive constant independent of the mesh size that may change at each occurrence. Moreover, given any subset $\omega$ in ${\mathbb{R}}^{2}$ and $k\in{\mathbb{N}}$, we will denote by $\mathbb{P}_{k}(\omega)$ the polynomials of total degree at most $k$ defined on $\omega$, with the extended notation $\mathbb{P}_{-1}(\omega)=\emptyset$. Let $\set{\mathcal{T}_{h}}_{h}$ be a sequence of decompositions of $\Omega$ into general polygonal elements $K$ with

[TABLE]

We suppose that for all $h$, each element $K$ in $\mathcal{T}_{h}$ fulfils the following assumptions:

• •

$\mathbf{(A1)}$ $K$ is star-shaped with respect to a ball of radius $\geq\,\gamma\,h_{K}$,

• •

$\mathbf{(A2)}$ the distance between any two vertexes of $K$ is $\geq c\,h_{K}$,

where $\gamma$ and $c$ are positive constants. We remark that the hypotheses above, though not too restrictive in many practical cases, can be further relaxed, as noted in  [5, 14]. From now on we assume that $\mathbb{K}$ is piecewise constant with respect to $\mathcal{T}_{h}$ on $\Omega$.

Using standard VEM notation, for $k\in\mathbb{N}$, let us define the spaces

• •

$\mathbb{P}_{k}(K)$ the set of polynomials on $K$ of degree $\leq k$,

• •

$\mathbb{B}_{k}(K):=\{v\in C^{0}(\partial K)\quad\text{s.t}\quad v_{|e}\in\mathbb{P}_{k}(e)\quad\forall\mbox{ edge }e\subset\partial K\}$,

• •

$\mathcal{G}_{k}(K):=\nabla(\mathbb{P}_{k+1}(K))\subseteq[\mathbb{P}_{k}(K)]^{2}$,

• •

$\mathcal{G}_{k}(K)^{\perp}:=\mathbf{x}^{\perp}[\mathbb{P}_{k-1}(K)]\subseteq[\mathbb{P}_{k}(K)]^{2}$ with $\mathbf{x}^{\perp}:=(x_{2},-x_{1})$.

In [15] the authors have introduced a new family of Virtual Elements for the Stokes problem on polygonal meshes. In particular, by a proper choice of the Virtual space of velocities, the virtual local spaces are associated to a Stokes-like variational problem on each element. The main ideas of the method are

• •

the Virtual space contains the space of all the polynomials of the prescribed order plus suitable non polynomial functions,

• •

the degrees of freedom are carefully chosen so that the $H^{1}$ semi-norm projection onto the space of polynomials can be exactly computed,

• •

the choice of the Virtual space of velocities and the associated degrees of freedom guarantee that the final discrete velocity is pointwise divergence-free and more generally the discrete kernel is contained in the continuous one.

In this section we briefly recall from [15] the notations, the main properties of the Virtual spaces and some details of the construction of the $H^{1}$ semi-norm projection. Let $k\geq 2$ the polynomial degree of accuracy of the method, then we define on each element $K\in\mathcal{T}_{h}$ the finite dimensional local virtual space

[TABLE]

where all the operators and equations above are to be interpreted in the distributional sense. It is easy to check that $[\mathbb{P}_{k}(K)]^{2}\subseteq\mathbf{W}_{h}^{K}$, and that (see [15] for the proof) the dimension of $\mathbf{W}_{h}^{K}$ is

[TABLE]

The corresponding degrees of freedom are chosen prescribing, given a function $\mathbf{v}\in\mathbf{W}_{h}^{K}$, the following linear operators $\mathbf{D_{V}}$, split into four subsets (see Figure 1):

• •

$\mathbf{D_{V}1}$: the values of $\mathbf{v}$ at the vertices of the polygon $K$,

• •

$\mathbf{D_{V}2}$: the values of $\mathbf{v}$ at $k-1$ distinct points of every edge $e\in\partial K$ (for example we can take the $k-1$ internal points of the $(k+1)$-Gauss-Lobatto quadrature rule in $e$, as suggested in [7]),

• •

$\mathbf{D_{V}3}$: the moments of $\mathbf{v}$

[TABLE]

• •

$\mathbf{D_{V}4}$: the moments up to order $k-1$ and greater than zero of ${\rm div}\,\mathbf{v}$ in $K$, i.e.

[TABLE]

For all $K\in\mathcal{T}_{h}$, we introduce the $H^{1}$ semi-norm projection ${\Pi}_{k}^{\nabla,K}\colon\mathbf{W}_{h}^{K}\to[\mathbb{P}_{k}(K)]^{2}$, defined by

[TABLE]

where $\Pi_{0}^{0,K}$ is the $L^{2}$-projection operator onto the constant functions defined on $K$. It is immediate to check that the energy projection is well defined and it clearly holds ${\Pi}_{k}^{\nabla,K}\mathbf{q}_{k}=\mathbf{q}_{k}$ for all $\mathbf{q}_{k}\in\mathbb{P}_{k}(K)$. Moreover the operator ${\Pi}_{k}^{\nabla,K}$ is computable in terms of the degrees of freedom $\mathbf{D_{V}}$ (see equations $(27)-(29)$ in [15] and the subsequent discussion).

3.2 The modified virtual space and the projection $\Pi^{0,K}_{k}$

Let $n$ a positive integer, then for all $K\in\mathcal{T}_{h}$, the $L^{2}$-projection $\Pi^{0,K}_{n}\colon\mathbf{W}^{K}_{h}\to[\mathbb{P}_{n}(K)]^{2}$ is defined by

[TABLE]

It is possible to check (see Section 3.3 of [15] for the proof) that the degrees of freedom $\mathbf{D_{V}}$ allow us to compute exactly the $L^{2}$-projection $\Pi^{0,K}_{k-2}$. On the other hand we can observe that we can not compute exactly from the DoFs the $L^{2}$-projection onto the space of polynomials of degree $\leq k$. The goal of the present section is to introduce, taking the inspiration from [1], a new virtual space $\mathbf{V}_{h}^{K}$ to be used in place of $\mathbf{W}_{h}^{K}$ in such a way that

• •

the DoFs $\mathbf{D_{V}}$ can still be used for $\mathbf{V}_{h}^{K}$,

• •

$[\mathbb{P}_{k}(K)]^{2}\subseteq\mathbf{V}_{h}^{K}$,

• •

the projection $\Pi_{k}^{0,K}\colon\mathbf{V}_{h}^{K}\to[\mathbb{P}_{k}(K)]^{2}$ can be exactly computable by the DoFs $\mathbf{D_{V}}$.

To construct $\mathbf{V}_{h}^{K}$ we proceed as follows: first of all we define an augmented virtual local space $\mathbf{U}_{h}^{K}$ by taking

[TABLE]

Now we define the enhanced Virtual Element space $\mathbf{V}_{h}^{K}$ as the restriction of $\mathbf{U}_{h}^{K}$ given by

[TABLE]

where the symbol $\mathcal{G}_{k}(K)^{\perp}/\mathcal{G}_{k-2}(K)^{\perp}$ denotes the polynomials in $\mathcal{G}_{k}(K)^{\perp}$ that are $L^{2}-$orthogonal to all polynomials of $\mathcal{G}_{k-2}(K)^{\perp}$. We proceed by investigating the dimension and by choosing suitable DoFs of the virtual space $\mathbf{V}_{h}$. First of all we recall from [9] the following facts

[TABLE]

where $n_{K}$ is the number of edges of the polygon $K$.

Lemma 3.1**.**

The dimension of $\mathbf{U}_{h}^{K}$ is

[TABLE]

Moreover as DoFs for $\mathbf{U}_{h}^{K}$ we can take the linear operators $\mathbf{D_{V}}$ and plus the moments

[TABLE]

Proof.

The proof is virtually identical to that given in [15] for $\mathbf{W}_{h}^{K}$ and it is based (see for instance [24]) on the fact that given

• •

a polynomial function $\mathbf{g}_{b}\in[\mathbb{B}_{k}({\partial K})]^{2}$,

• •

a polynomial function $\mathbf{h}\in\mathcal{G}_{k}(K)^{\perp}$,

• •

a polynomial function $g\in\mathbb{P}_{k-1}(K)$ satisfying the compatibility condition

[TABLE]

there exists a unique pair $(\mathbf{v},s)\in\mathbf{U}_{h}^{K}\times L^{2}(K)/\mathbb{R}$ such that

[TABLE]

Moreover, since from [9], ${\rm rot}\colon\mathcal{G}_{k}(K)^{\perp}\to\mathbb{P}_{k-1}(K)$ is an isomorphism, we can conclude that the map that associates a given compatible data set $(\mathbf{g}_{b},\,\mathbf{h},\,g)$ to the velocity field $\mathbf{v}$ that solves (12) is an injective map. Then

[TABLE]

and the thesis follows from (11). ∎

Proposition 3.1**.**

The dimension of $\mathbf{V}_{h}^{K}$ is equal to that of $\mathbf{W}_{h}^{K}$ that is, as in (8)

[TABLE]

As DoFs in $\mathbf{V}_{h}^{K}$ we can take $\mathbf{D_{V}}$.

Proof.

From (11) it is straightforward to check that

[TABLE]

Hence, neglecting the independence of the additional $2k-1$ conditions in (10), it holds that

[TABLE]

We now observe that a function $\mathbf{v}\in\mathbf{V}_{h}^{K}$ such that $\mathbf{D_{V}}(\mathbf{v})=0$ is identically zero. Indeed, from (9), it is immediate to check that in this case the $\Pi_{k}^{\nabla,K}\,\mathbf{v}$ would be zero, implying that all its moment are zero, in particular, since $\mathbf{v}\in\mathbf{V}_{h}^{K}$, all the moments $\mathbf{D_{U}}$ of $\mathbf{v}$ are also zero. Now, from Lemma 3.1, we have that $\mathbf{v}$ is zero. Therefore, from (14), we obtain that the dimension of $\mathbf{V}_{h}^{K}$ is actually the same of $\mathbf{W}_{h}^{K}$, and that the DoFs $\mathbf{D_{V}}$ are unisolvent for $\mathbf{V}_{h}^{K}$. ∎

Proposition 3.2**.**

The degrees of freedom $\mathbf{D_{V}}$ allow us to compute exactly the $L^{2}$-projection $\Pi^{0,K}_{k}\colon\mathbf{V}_{h}\to[\mathbb{P}_{k}(K)]^{2}$, i.e. the moments

[TABLE]

for all $\mathbf{v}\in\mathbf{V}_{h}$ and for all $\mathbf{q}_{k}\in[\mathbb{P}_{k}(K)]^{2}$.

Proof.

Let us set

[TABLE]

with $q_{k+1}\in\mathbb{P}_{k+1}(K)/\mathbb{R}$, $\mathbf{g}_{k-2}^{\perp}\in\mathcal{G}_{k-2}^{\perp}(K)$ and $\mathbf{g}_{k}^{\perp}\in\mathcal{G}_{k}^{\perp}(K)/\mathcal{G}_{k-2}^{\perp}(K)$. Therefore using the Green formula and since $\mathbf{v}\in\mathbf{V}_{h}$, we get

[TABLE]

Now, since ${\rm div}\,\mathbf{v}$ is a polynomial of degree less or equal than $k-1$ we can reconstruct its value from $\mathbf{D_{V}4}$ and compute exactly the first term. The second term is computable from $\mathbf{D_{V}3}$. The third term is computable from all the $\mathbf{D_{V}}$ using the projection $\Pi_{k}^{\nabla,K}\,\mathbf{v}$. Finally from $\mathbf{D_{V}1}$ and $\mathbf{D_{V}2}$ we can reconstruct $\mathbf{v}$ on the boundary and so compute exactly the boundary term. ∎

For what concerns the pressures we take the standard finite dimensional space

[TABLE]

having dimension

[TABLE]

The corresponding degrees of freedom are chosen defining for each $q\in Q_{h}^{K}$ the following linear operators $\mathbf{D_{Q}}$:

• •

$\mathbf{D_{Q}}$: the moments up to order $k-1$ of $q$, i.e.

[TABLE]

Finally we define the global virtual element spaces as

[TABLE]

and

[TABLE]

with the obvious associated sets of global degrees of freedom. A simple computation shows that:

[TABLE]

and

[TABLE]

where $n_{P}$ is the number of elements, $n_{E}$, $n_{V}$ (resp., $n_{E,B}$, $n_{V,B}$) is the number of internal edges and vertexes (resp., boundary edges and vertexes) in $\mathcal{T}_{h}$. As observed in [15], we remark that

[TABLE]

Remark 3.1*.*

By definition (16) it is clear that our discrete velocities field is $H^{1}$-conforming, in particular we obtain continuous velocities, whereas the natural discretization is only $H({\rm div})$- conforming. This property, in combination with (18), will make our method suitable for a (robust) extension to the Brinkman problem.

3.3 The discrete bilinear forms

The next step in the construction of our method is to define on the virtual spaces $\mathbf{V}_{h}$ and $Q_{h}$ a discrete version of the bilinear forms $a(\cdot,\cdot)$ and $b(\cdot,\cdot)$ given in (2) and (3). For simplicity we assume that the tensor $\mathbb{K}$ is piecewise constant with respect to the decomposition $\mathcal{T}_{h}$, i.e. $\mathbb{K}$ is constant on each polygon $K\in\mathcal{T}_{h}$. First of all we decompose into local contributions the bilinear forms $a(\cdot,\cdot)$ and $b(\cdot,\cdot)$, the norms $\|\cdot\|_{\mathbf{V}}$ and $\|\cdot\|_{Q}$ by defining

[TABLE]

[TABLE]

and

[TABLE]

We now define discrete versions of the bilinear form $a(\cdot,\cdot)$ (cf. (2)), and of the bilinear form $b(\cdot,\cdot)$ (cf. (3)). For what concerns $b(\cdot,\cdot)$, we simply set

[TABLE]

i.e. as noticed in [15] we do not introduce any approximation of the bilinear form. We notice that (19) is computable from the degrees of freedom $\mathbf{D_{V}1}$, $\mathbf{D_{V}2}$ and $\mathbf{D_{V}4}$, since $q$ is polynomial in each element $K\in\mathcal{T}_{h}$. On the other hand, the bilinear form $a(\cdot,\cdot)$ needs to be dealt with in a more careful way. First of all, by Proposition 3.2, we observe that for all $\mathbf{q}_{k}\in[\mathbb{P}_{k}(K)]^{2}$ and for all $\mathbf{v}\in\mathbf{V}_{h}^{K}$, the quantity

[TABLE]

is exactly computable by the DoFs. However, for an arbitrary pair $(\mathbf{u},\mathbf{v})\in\mathbf{V}_{h}^{K}\times\mathbf{V}_{h}^{K}$, the quantity $a_{h}^{K}(\mathbf{w},\mathbf{v})$ is clearly not computable. In the standard procedure of VEM framework, we define a computable discrete local bilinear form

[TABLE]

approximating the continuous form $a^{K}(\cdot,\cdot)$ and satisfying the following properties:

• •

$\mathbf{k}$-consistency: for all $\mathbf{q}_{k}\in[\mathbb{P}_{k}(K)]^{2}$ and $\mathbf{v}_{h}\in\mathbf{V}_{h}^{K}$

[TABLE]

• •

stability: there exist two positive constants $\alpha_{*}$ and $\alpha^{*}$, independent of $h$ and $K$, such that, for all $\mathbf{v}_{h}\in\mathbf{V}_{h}^{K}$, it holds

[TABLE]

Let $\mathcal{R}^{K}\colon\mathbf{V}_{h}^{K}\times\mathbf{V}_{h}^{K}\to\mathbb{R}$ be a (symmetric) stabilizing bilinear form, satisfying

[TABLE]

with $c_{*}$ and $c^{*}$ positive constants independent of $h$ and $K$. Then, we can set

[TABLE]

for all $\mathbf{u}_{h},\mathbf{v}_{h}\in\mathbf{V}_{h}^{K}$.

It is straightforward to check that Definition (9) and properties (23) imply the consistency and the stability of the bilinear form $a_{h}^{K}(\cdot,\cdot)$.

Remark 3.2*.*

In the construction of the stabilizing form $\mathcal{R}^{K}$ with condition (23) we essentially require that the stabilizing term $\mathcal{R}^{K}(\mathbf{v}_{h},\mathbf{v}_{h})$ scales as $a^{K}(\mathbf{v}_{h},\mathbf{v}_{h})$. Following the standard VEM technique (cf. [5, 7] for more details), denoting with $\bar{\mathbf{u}}_{h}$, $\bar{\mathbf{v}}_{h}\in\mathbb{R}^{N_{K}}$ the vectors containing the values of the $N_{K}$ local degrees of freedom associated to $\mathbf{u}_{h},\mathbf{v}_{h}\in\mathbf{V}_{h}^{K}$, we set

[TABLE]

where $\alpha^{K}$ is a suitable positive constant that scales as $|K|$. For example, in the numerical tests presented in Section 6, we have chosen $\alpha^{K}$ as the mean value of the eigenvalues of the matrix stemming from the term $a^{K}\left({\Pi}_{k}^{0,K}\mathbf{u}_{h},\,{\Pi}_{k}^{0,K}\mathbf{v}_{h}\right)$ in (24).

Finally we define the global approximated bilinear form $a_{h}(\cdot,\cdot)\colon\mathbf{V}_{h}\times\mathbf{V}_{h}\to\mathbb{R}$ by simply summing the local contributions:

[TABLE]

3.4 The discrete problem

We are now ready to state the proposed discrete problem. Referring to (16), (17),  (19), and  (25) we consider the virtual element problem:

[TABLE]

We point out that the symmetry of $a_{h}(\cdot,\cdot)$ together with (22) easily implies that $a_{h}(\cdot,\cdot)$ is (uniformly) continuous with respect to the $L^{2}$ norm. Moreover, as observed in [15], introducing the discrete kernel:

[TABLE]

it is immediate to check that

[TABLE]

Then the bilinear form $a_{h}(\cdot,\cdot)$ is also uniformly coercive on the discrete kernel $\mathbf{Z}_{h}$ with respect to the $\mathbf{V}$ norm. Moreover as a direct consequence of Proposition 4.3 in [15], we have the following stability result.

Proposition 3.3**.**

Given the discrete spaces $\mathbf{V}_{h}$ and $Q_{h}$ defined in (16) and (17), there exists a positive $\tilde{\beta}$, independent of $h$, such that:

[TABLE]

In particular, the the inf-sup condition of Proposition 3.3, along with property (18), implies that:

[TABLE]

Finally we can state the well-posedness of virtual problem (26).

Theorem 3.1**.**

Problem (26) has a unique solution $(\mathbf{u}_{h},p_{h})\in\mathbf{V}_{h}\times Q_{h}$, verifying the estimate

[TABLE]

4 Theoretical results

We begin by proving an approximation result for the virtual local space $\mathbf{V}_{h}$. First of all, let us recall a classical result by Brenner-Scott (see [22]).

Lemma 4.1**.**

Let $K\in\mathcal{T}_{h}$, then for all $\mathbf{u}\in[H^{s+1}(K)]^{2}$ with $0\leq s\leq k$, there exists a polynomial function $\mathbf{u}_{\pi}\in[{\mathbb{P}}_{k}(K)]^{2}$, such that

[TABLE]

We have the following approximation results (for the proof see [16]).

Proposition 4.1**.**

Let $\mathbf{u}\in\mathbf{V}\cap[H^{s+1}(\Omega)]^{2}$ with $0\leq s\leq k$. Under the assumption $\mathbf{(A1)}$ and $\mathbf{(A2)}$ on the decomposition $\mathcal{T}_{h}$, there exists $\mathbf{u}_{int}\in\mathbf{W}_{h}$ such that

[TABLE]

where $C$ is a constant independent of $h$.

For what concerns the pressures, from classic polynomial approximation theory [22], for $q\in H^{k}(\Omega)$ it holds

[TABLE]

We are ready to state the following convergence theorem.

Theorem 4.1**.**

Let $(\mathbf{u},p)\in\mathbf{V}\times Q$ be the solution of problem (4) and $(\mathbf{u}_{h},p_{h})\in\mathbf{V}_{h}\times Q_{h}$ be the solution of problem (26). Then it holds

[TABLE]

Proof.

We begin by remarking that as a consequence of the inf-sup condition with classical arguments (see for instance Proposition 2.5 in [24]), there exists $\mathbf{u}_{I}\in\mathbf{V}_{h}$ such that

[TABLE]

Let us set $\boldsymbol{\delta}_{h}=\mathbf{u}_{I}-\mathbf{u}_{h}$. From (30) and (26), we have that ${\rm div}\,\boldsymbol{\delta}_{h}=0$ and thus $\boldsymbol{\delta}_{h}\in\mathbf{Z}_{h}$. Now, using (5), (22), (26) and introducing the piecewise polynomial approximation (28) together with (21), we have

[TABLE]

then

[TABLE]

The $L^{2}$-estimate follows easily by the triangle inequality. It is also straightforward to see from (4) and (26) that

[TABLE]

than we get ${\rm div}\,\mathbf{u}_{h}=\Pi_{k-1}^{0,K}({\rm div}\,\mathbf{u})$ for all $K\in\mathcal{T}_{h}$ and therefore

[TABLE]

from which the estimate in the ${\mathbf{V}}$ norm. We proceed by analysing the error on the pressure field. Let $q_{h}\in Q_{h}$, then from the discrete inf-sup condition (27), we infer:

[TABLE]

Since $(\mathbf{u},p)$ and $(\mathbf{u}_{h},p_{h})$ are respectively the solution of (4) and (26), it follows that

[TABLE]

Therefore, we get

[TABLE]

Using (21), the continuity of $a_{h}(\cdot,\cdot)$ and the triangle inequality, we get:

[TABLE]

where $\mathbf{u}_{\pi}$ is the piecewise polynomial of degree $k$ defined in Lemma 4.1. Then, from estimate (28) and the previous estimate on the velocity error, we obtain

[TABLE]

Moreover, we have

[TABLE]

Then, using (33) and (34) in (32), we infer

[TABLE]

Finally, using (35) and the triangular inequality, we get

[TABLE]

Passing to the infimum with respect to $q_{h}\in Q_{h}$, and using estimate (29), we get the thesis. ∎

Remark 4.1*.*

We observe that the estimates on the velocity errors in Theorem 4.1 do not depend on the continuous pressure, whereas the velocity errors of the classical methods have a pressure contribution. Therefore the proposed scheme belongs to the class of the pressure-robust methods.

5 A Stable VEM for Brinkman Equations

5.1 The continuous problem

The Brinkman equation describes fluid flow in complex porous media with a viscosity coefficient highly varying so that the flow is dominated by the Darcy equations in some regions of the domain and by the Stokes equation in others. We consider the Brinkman equation on a polygon $\Omega\subseteq\mathbb{R}^{2}$ with homogeneous Dirichlet boundary conditions:

[TABLE]

where $\mathbf{u}$ and $p$ are the unknown velocity and pressure fields, $\mu$ is the fluid viscosity, $\mathbb{K}$ denotes the permeability tensor of the porous media and $\mathbf{f}\in[L^{2}(\Omega)]^{2}$ is the external source term. We assume that $\mathbb{K}$ is a symmetric positive definite tensor and that there exist two positive (uniform) constants $\lambda_{1}$, $\lambda_{2}>0$ such that

[TABLE]

For what concerns the fluid viscosity we consider $0<\mu\leq C$, this include the case where $\mu$ approaches zero and equation (36) becomes a singular perturbation of the classic Darcy equations. Let us consider the spaces

[TABLE]

with the usual norms, and let $A(\cdot,\cdot)\colon\mathbf{V}\times\mathbf{V}\to\mathbb{R}$ be the bilinear form defined by:

[TABLE]

where

[TABLE]

and $a(\cdot,\cdot)$ is the bilinear form defined in (2). Then the variational formulation of Problem (36) is:

[TABLE]

where and $b(\cdot,\cdot)\colon\mathbf{V}\times Q\to\mathbb{R}$ is the bilinear form in (3) and using standard notation

[TABLE]

The natural energy norm for the velocities is induced by the symmetric an positive bilinear form $A(\cdot,\cdot)$ and is defined by (e.g. [35])

[TABLE]

We can observe that the equivalence with the $\mathbf{V}$ norm is not uniform, i.e.

[TABLE]

where $c_{1}$, $c_{2}$ here and in the follows denote two positive constant independent of $h$ and $\mu$. For what concerns the pressures, we consider the norm (see for instance [35])

[TABLE]

Using the inf-sup condition in the usual norm it is possible to check the equivalence between the norms for the pressure but again the equivalence is not uniform, i.e.

[TABLE]

Since, considering the modified norm, the bilinear form $A(\cdot,\cdot)$ is uniformly continuous and coercive, and inf-sup condition is clearly fulfilled, Problem (37) has a unique solution $(\mathbf{u},p)\in\mathbf{V}\times Q$ such that

[TABLE]

where the constant $C$ depends only on $\Omega$.

5.2 Virtual formulation for Brinkman equations

Mathematically, Brinkman equations can be viewed as a combination of the Stokes and the Darcy equation, that can change from place to place in the computational domain. Therefore, numerical schemes for Brinkman equations have to be carefully designed to accommodate both Stokes and Darcy simultaneously. In this section we propose a Virtual Element scheme that is accurate for both Darcy and Stokes flows. For this goal we combine the ideas developed in the previous sections with the argument in [15].

Let us consider the virtual spaces $\mathbf{V}_{h}$ and $Q_{h}$ (cfr. (16) and (17)). As usual in the VEM framework we need to define a computable approximation of the continuous bilinear forms. Using obvious notations we split the bilinear form $A(\cdot,\cdot)$ as

[TABLE]

We begin by observing that, from [15] (in particular c.f. $(27)-(29)$) and from Section 3.2, $A^{K}(\mathbf{q}_{k},\mathbf{v})$ is computable on the basis of the DoFs $\mathbf{D_{V}}$ for all $\mathbf{q}_{k}\in[\mathbb{P}_{k}(K)]^{2}$ and for all $\mathbf{v}\in\mathbf{V}_{h}$. Starting from this observation we can approximate the continuous form $A^{K}(\cdot,\cdot)$ with the bilinear form $A_{h}^{K}(\cdot,\cdot)\colon\mathbf{V}_{h}^{K}\times\mathbf{V}_{h}^{K}\to\mathbb{R}$, given by

[TABLE]

where $a_{h}^{\nabla,K}(\cdot,\cdot)$ is the bilinear form defined in equation $(35)$ in [15] and $a_{h}^{K}(\cdot,\cdot)$ is defined in (20). It is clear that the bilinear form $A_{h}^{K}(\cdot,\cdot)$ satisfies the $k$-consistency and the stability properties. As usual we build the global approximated bilinear form $A_{h}(\cdot,\cdot)\colon\mathbf{V}_{h}\times\mathbf{V}_{h}\to\mathbb{R}$ by simply summing the local contributions. For what concerns the bilinear form $b(\cdot,\cdot)$, as observed in Section 3.3, it can be computed exactly. The last step consists in constructing a computable approximation of the right-hand side $(\mathbf{f},\,\mathbf{v})$ in (37). We define the approximated load term $\mathbf{f}_{h}$ as

[TABLE]

and consider:

[TABLE]

We observe that (40) can be exactly computed from $\mathbf{D_{V}}$ for all $\mathbf{v}_{h}\in\mathbf{V}_{h}$ (see Proposition 3.2). Furthermore, the following result concerning a $L^{2}$ and $H^{1}$-type norm, can be proved using standard arguments [5].

Lemma 5.1**.**

Let $\mathbf{f}_{h}$ be defined as in (39), and let us assume $\mathbf{f}\in H^{k+1}(\Omega)$. Then, for all $\mathbf{v}_{h}\in\mathbf{V}_{h}$, it holds

[TABLE]

In the light of the previous definitions, we consider the virtual element approximation of the Brinkman problem:

[TABLE]

Equation (41) is well posed since the discrete bilinear form $A_{h}(\cdot,\cdot)$ is (uniformly) stable with respect to the norm $\|\cdot\|_{\mathbf{V},\mu}$ by construction and the inf-sup condition is fulfilled (the proof follows the guidelines of Proposition 4.2 in [15] and the linearity of the Fortin operator). Then we have the following result.

Theorem 5.1**.**

Problem (41) has a unique solution $(\mathbf{u}_{h},p_{h})\in\mathbf{V}_{h}\times Q_{h}$, verifying the estimate

[TABLE]

We now notice that, if $\mathbf{u}\in\mathbf{V}$ is the velocity solution to Problem (37), then it is the solution to Problem:

[TABLE]

Analogously, if $\mathbf{u}_{h}\in\mathbf{V}_{h}$ is the velocity solution to Problem (41), then it is the solution to Problem:

[TABLE]

For what concerns the convergence results we state the following theorem. The proof can be derived by extending the techniques of the previous section and is therefore omitted.

Theorem 5.2**.**

Let $\mathbf{u}\in\mathbf{Z}$ be the solution of problem (42) and $\mathbf{u}_{h}\in\mathbf{Z}_{h}$ be the solution of problem (43). Then

[TABLE]

Let $(\mathbf{u},p)\in\mathbf{V}\times Q$ be the solution of Problem (37) and $(\mathbf{u}_{h},p_{h})\in\mathbf{W}_{h}\in Q_{h}$ be the solution of Problem (41). Then it holds:

[TABLE]

The constants $C$ above are independent of $h$ and $\mu$.

In the last part of this section we present a brief discussion about the construction of a reduced virtual element method for Brinkman equations equivalent to Problem (41) but involving significantly fewer degrees of freedom, especially for large $k$. This construction essentially follows the guidelines of Section 5 in [15] (where we refer the reader for a deeper presentation). Let us define the original reduced local virtual spaces, for $k\geq 2$:

[TABLE]

As before we enlarge the virtual space $\widehat{\mathbf{W}}_{h}^{K}$ and we consider

[TABLE]

Finally we define the enhanced Virtual Element space, the restriction $\widehat{\mathbf{V}}_{h}^{K}$ of $\widehat{\mathbf{U}}_{h}^{K}$ given by

[TABLE]

where as before the symbol $\mathcal{G}_{k}(K)^{\perp}/\mathcal{G}_{k-2}(K)^{\perp}$ denotes the polynomials i $\mathcal{G}_{k}(K)^{\perp}$ that are $L^{2}-$orthogonal to all polynomials of $\mathcal{G}_{k-2}(K)^{\perp}$. For the pressures we consider the reduced space

[TABLE]

As sets of degrees of freedom for the reduced spaces, combining the argument in Section 3 and [15] we may consider the following. For every function $\mathbf{v}\in\widehat{\mathbf{V}}_{h}^{K}$ we take the following linear operators $\mathbf{\widehat{D}_{V}}$, split into three subsets (see Figure 2):

• •

$\mathbf{\widehat{D}_{V}1}$: the values of $\mathbf{v}$ at each vertex of the polygon $K$,

• •

$\mathbf{\widehat{D}_{V}2}$: the values of $\mathbf{v}$ at $k-1$ distinct points of every edge $e\in\partial K$,

• •

$\mathbf{\widehat{D}_{V}3}$: the moments of $\mathbf{v}$

[TABLE]

For every $q\in\widehat{Q}_{h}$ we consider

• •

$\mathbf{\widehat{D}_{Q}}$: the moment

[TABLE]

Therefore we have that:

[TABLE]

and

[TABLE]

where $n_{K}$ is the number of vertexes in $K$.

We define the global reduced virtual element spaces in the standard fashion. The reduced virtual element discretization of the Brinkman problem (37) is then:

[TABLE]

Above, the bilinear forms $A_{h}(\cdot,\cdot)$ and $b(\cdot,\cdot)$, and the loading term $\mathbf{f}_{h}$ are the same as before. The following proposition states the relation between Problem (41) and the reduced Problem (44) (the proof is equivalent to that of Proposition 5.1 in [15]).

Proposition 5.1**.**

Let $(\mathbf{u}_{h},p_{h})\in\mathbf{V}_{h}\times Q_{h}$ be the solution of problem (41) and $(\widehat{\mathbf{u}}_{h},\widehat{p}_{h})\in\widehat{\mathbf{V}}_{h}\times\widehat{Q}_{h}$ be the solution of problem (44). Then

[TABLE]

6 Numerical tests

In this section we present two numerical experiments to test the practical performance of the method. The first experiment is focused on the method introduced in Section 3 for the Darcy Problem, whereas in the second experiment we test the method in Section 5 for the Brinkman equations.

Since the VEM velocity solution $\mathbf{u}_{h}$ is not explicitly known point-wise inside the elements, we compute the method error comparing $\mathbf{u}$ with a suitable polynomial projection of the approximated $\mathbf{u}_{h}$. In particular we consider the computable error quantities:

[TABLE]

Regarding the computational domain, in our tests we always take the square domain $\Omega=[0,1]^{2}$, which is partitioned using the following sequences of polygonal meshes:

• •

$\{\mathcal{V}_{h}\}_{h}$: sequence of Voronoi meshes with $h=1/4,1/8,1/16,1/32$,

• •

$\{\mathcal{T}_{h}\}_{h}$: sequence of triangular meshes with $h=1/2,1/4,1/8,1/16$,

• •

$\{\mathcal{Q}_{h}\}_{h}$: sequence of square meshes with $h=1/4,1/8,1/16,1/32$.

• •

$\{\mathcal{W}_{b}\}_{h}$: sequence of WEB-like meshes with $h=4/10,2/10,1/10,1/20$.

An example of the adopted meshes is shown in Figure 3.

For the generation of the Voronoi meshes we use the code Polymesher [46]. The non convex WEB-like meshes are composed by hexagons, generated starting from the triangular meshes $\{\mathcal{T}_{h}\}_{h}$ and randomly displacing the midpoint of each (non boundary) edge.

Test 6.1*.*

In this example we consider the Darcy problem (4) where we set $\mathbb{K}=I$, and we choose the load term $\mathbf{f}$ in such a way that the analytical solution is

[TABLE]

We analyse the practical performance of the virtual method by studying the errors versus the diameter $h$ of the meshes. In addition we compare the results obtained with the scheme of Section 3, labeled as “div-free”, with those obtained with the method in the Appendix, labeled as “non div-free” (in both cases we consider polynomial degrees $k=2$).

Remark 6.1*.*

We notice that the "non div-free" method is a naive extension to the Darcy equation of the inf-sup stable scheme proposed in [6]. Since the scheme lacks a uniform ellipticity-on-the-kernel condition, it is not recommended for the problem under consideration. The purpose of the comparison is thus to underline the importance of the property $Z_{h}\subseteq Z$ (cf. Section 3.4) in the present context.

In Figure 4 and 5, we display the results for the sequence of Voronoi meshes $\mathcal{V}_{h}$. In Figure 6 and 7, we show the results for the sequence of meshes $\mathcal{T}_{h}$, while in Figure 8 and 9 we plot the results for the sequence of meshes $\mathcal{Q}_{h}$, finally in 10 and 11 we exhibit the results for the sequence of meshes $\mathcal{W}_{h}$.

We notice that the theoretical predictions of Section 4 and the Appendix are confirmed for both the $L^{2}$ norm and the $H({\rm div})$ norm. Note that for the $H({\rm div})$ norm we plot only the error for the “div-free” method since such scheme guarantees by construction, a better approximation of the divergence. Indeed let $u_{h}$ (resp. $\widetilde{u}_{h}$) be the solution obtained with the “div-free” method (“non div-free” method) then $u_{h}$ satisfies

[TABLE]

whereas $\widetilde{u}_{h}$ satisfies the same equation only in a projected sense, i.e.

[TABLE]

We can observe that the convergence rate of the $L^{2}$ norm for the pressure is optimal also for the “non div-free” method as proved in the Appendix. Finally, we can observe that using a square mesh decomposition holds a convergence rate that is slightly better than what predicted by the theory.

Test 6.2*.*

In this example we test the Brinkman equation (37) with different values of the fluid viscosity $\mu$ and fixed permeability tensor $\mathbb{K}=I$. We choose the load term $\mathbf{f}$ and the Dirichlet boundary conditions in such a way that the analytical solution is

[TABLE]

The aim of this test is to check the practical performance of the method introduced in Section 5 in the reduced formulation (c.f. (44). In Table 1 and Table 2 we display the total amount of DoFs and the errors for the family of meshes $\mathcal{V}_{h}$ choosing $k=2$ respectively for the “div-free” method (cf. Section 5) and the “non div-free” method (cf. Reamrk 6.1 and the Appendix). We observe that also in the limit case, when the equation becomes a singular perturbation of the classic Darcy equations (e.g. for “small” $\mu$), the proposed “div-free” method preserves the optimal order of accuracy.

7 Acknowledgements

The author wishes to thank L. Beirão da Veiga and C. Lovadina for several interesting discussions and suggestions on the paper. The author was partially supported by the European Research Council through the H2020 Consolidator Grant (grant no. 681162) CAVE, Challenges and Advancements in Virtual Elements. This support is gratefully acknowledged.

Appendix: Non divergence-free virtual space

We have built a new $H^{1}$-conforming (vector valued) virtual space for the velocity vector field different from the more standard one presented in [6] for the elasticity problem. The topic of the present section is to analyse the extension to the Darcy equation of the scheme of [6]. Even though the method should not be used for the Darcy problem (cf. Remark 6.1), the numerical experiments have shown an optimal error convergence rate for the pressure variable. In this Section, we theoretically explain such a behaviour, under a convexity assumption on $\Omega$ (essentially, a regularity assumption on the problem). To this end, we develop an inverse estimate for the VEM spaces which is interesting on its own, and can be used in other contexts. We briefly describe the method by making use of various tools from the Virtual Element technology, and we refer the interested reader to the papers [5, 1, 7, 6]) for a deeper presentation. We consider the local virtual space

[TABLE]

with local degrees of freedom $\mathbf{\widetilde{D}_{V}}$:

• •

$\mathbf{\widetilde{D}_{V}1}$: the values of $\mathbf{v}$ at each vertex of the polygon $K$,

• •

$\mathbf{\widetilde{D}_{V}2}$: the values of $\mathbf{v}$ at $k-1$ distinct points of every edge $e\in\partial K$,

• •

$\mathbf{\widetilde{D}_{V}3}$: the moments of $\mathbf{v}$ up to order $k-2$, i.e.

[TABLE]

As observed in [1], the DoFs $\mathbf{\widetilde{D}_{V}}$ allow us to compute the operator $\widetilde{{\Pi}}_{k}^{\nabla,K}\colon\widetilde{\mathbf{W}}_{h}^{K}\to[\mathbb{P}_{k}(K)]^{2}$ defined as the analogous of the $H^{1}$ semi-norm projection (c.f. (9)). For all $K\in\mathcal{T}_{h}$, the augmented virtual local space $\widetilde{\mathbf{U}}_{h}^{K}$ is defined by

[TABLE]

Now we define the enhanced Virtual Element space, the restriction $\mathbf{\widetilde{V}}_{h}^{K}$ of $\widetilde{\mathbf{U}}_{h}^{K}$ given by

[TABLE]

where the symbol $\mathbb{P}_{k}(K)/\mathbb{P}_{k-2}(K)$ denotes the polynomials of degree $k$ living on $K$ that are $L^{2}-$orthogonal to all polynomials of degree $k-2$ on $K$. The enhanced space $\widetilde{\mathbf{V}}_{h}^{K}$ has three fundamental properties (see [1] for a proof):

• •

$[\mathbb{P}_{k}(K)]^{2}\subseteq\widetilde{\mathbf{V}}_{h}^{K}$,

• •

the set of linear operators $\mathbf{\widetilde{D}_{V}}$ constitutes a set of DoFs for the space $\widetilde{\mathbf{V}}_{h}^{K}$,

• •

the $L^{2}$-projection operator $\widetilde{\Pi}^{0,K}_{k}\colon\widetilde{\mathbf{V}}_{h}^{K}\to[\mathbb{P}_{k}(K)]^{2}$ is exactly computable by the DoFs.

Recalling (13) and from [6] it holds that ${\rm dim}\left(\widetilde{\mathbf{V}}_{h}^{K}\right)={\rm dim}\left(\mathbf{V}_{h}^{K}\right)$. For the pressures we use the space of the piecewise polynomials $Q_{h}$ (c.f. (17)). For what concerns the construction of the approximated bilinear forms, it is straightforward to see that

[TABLE]

is computable from the DoFs for all $\mathbf{v}\in\widetilde{\mathbf{V}}_{h}$, $q\in Q_{h}$. Moreover using standard arguments [1, 48] we can define a computable bilinear form

[TABLE]

approximating the continuous form $a^{K}(\cdot,\cdot)$, and satisfying the $k$-consistency (c.f. (21)) and the stability properties (c.f. (22)). Finally we define the global approximated bilinear form $\widetilde{a}_{h}(\cdot,\cdot)\colon\widetilde{\mathbf{V}}_{h}\times\widetilde{\mathbf{V}}_{h}\to\mathbb{R}$ by simply summing the local contributions. By construction (see for instance [1]) the discrete bilinear form $\widetilde{a}_{h}(\cdot,\cdot)$ is (uniformly) stable with respect to the $L^{2}$ norm. We are now ready to state the proposed discrete virtual element problem:

[TABLE]

We shall first prove an inverse inequality for the virtual element functions in $\widetilde{\mathbf{V}}_{h}$.

Lemma 7.1**.**

Under the assumption $\mathbf{(A1)}$, $\mathbf{(A2)}$, let $K\in\mathcal{T}_{h}$ and let $\mathbf{v}_{h}\in\widetilde{\mathbf{V}}_{h}^{K}$. Then the following inverse estimate holds

[TABLE]

where the constant $c_{inv}$ is independent of $\mathbf{v}_{h}$, $h_{K}$ and $K$.

Proof.

We only sketch the proof, since we follow the guidelines of Lemma 3.1 and 3.3 in [14]. Let $\mathbf{v}_{h}\in\widetilde{\mathbf{V}}_{h}^{K}$, then

[TABLE]

Under the assumption $\mathbf{(A1)}$, $\mathbf{(A2)}$ and by Lemma 3.3 in [14] we get

[TABLE]

where the constant $C_{1}$ is independent of $\mathbf{v}_{h}$, $h_{K}$ and $K$. For what concerns the second addend in the right side of (47), under the assumptions $\mathbf{(A1)}$, $\mathbf{(A2)}$, and using Lemma 3.1 in [14], for all $\mathbf{w}\in[H^{1/2}(\partial K)]^{2}$ the following holds: there exists an extension $\widetilde{\mathbf{w}}\in[H^{1}(K)]^{2}$ of $\mathbf{w}$ such that

[TABLE]

where we consider the scaled norm

[TABLE]

By definition it holds

[TABLE]

Now, using the definition (50), an inverse estimate ($\mathbf{v}_{h}$ is polynomial on $\partial K$) and the trace theorem [22], it holds that

[TABLE]

for any real $\varepsilon>0$. For the last term, using (50), (49) and (48) we get

[TABLE]

From (51) and (52) we can conclude that

[TABLE]

Finally, choosing $\varepsilon=\frac{1}{2C_{3}}$ and collecting (48) and (53) in (47) we have

[TABLE]

from which follows the thesis. ∎

Let us analyse the theoretical properties of the method. We consider the discrete kernel:

[TABLE]

therefore the divergence-free property is satisfied only in a relaxed (projected) sense. As a consequence the bilinear form $\widetilde{a}_{h}(\cdot,\cdot)$ is not uniformly coercive on the discrete kernel $\widetilde{\mathbf{Z}}_{h}$; nevertheless it holds the following $h$-dependent coercivity property

[TABLE]

that can be derived by using inverse estimate (46). Recalling that $\widetilde{a}_{h}(\cdot,\cdot)$ is continuous with respect the ${\mathbf{V}}$ norm and that the discrete inf-sup condition is fulfilled [6]

[TABLE]

problem (45) has a unique solution but we expect a worse order of accuracy since the bilinear form $\widetilde{a}_{h}(\cdot,\cdot)$ is not uniformly stable. In fact we have the following convergence results that are, perhaps surprisingly, still optimal in the pressure variable.

Theorem 7.1**.**

Let $(\mathbf{u},p)\in\mathbf{V}\times Q$ be the solution of problem (4) and $(\widetilde{\mathbf{u}}_{h},\widetilde{p}_{h})\in\widetilde{\mathbf{V}}_{h}\times Q_{h}$ be the solution of problem (45). Then

[TABLE]

Assuming further that $\Omega$ is convex, the following estimate holds:

[TABLE]

Proof.

As observed in the proof of Theorem 4.1 the inf-sup condition (55) implies the existence of a function $\widetilde{\mathbf{u}}_{I}\in\widetilde{\mathbf{V}}_{h}$ such that

[TABLE]

Now let us set $\boldsymbol{\delta}_{h}=\widetilde{\mathbf{u}}_{h}-\widetilde{\mathbf{u}}_{I}$. For what concerns the $L^{2}$ norm, using the stability of the bilinear form $\widetilde{a}_{h}(\cdot,\cdot)$ and (45) together with (4)

[TABLE]

By (45) and property (56), it is straightforward to see that

[TABLE]

so that $\boldsymbol{\delta}_{h}\in\widetilde{\mathbf{Z}}_{h}$. Therefore

[TABLE]

for all $q_{h}\in Q_{h}$. Using the inverse estimate (46) and standard approximation theory we get

[TABLE]

By standard technique in VEM convergence theory, it holds that

[TABLE]

Collecting (59) and (60) in (58) we get the $L^{2}$ estimate. Whereas the $\mathbf{V}$ norm estimate follows from an inverse estimate (46).

For what concerns the estimate on the pressure, let $p_{\pi}$ the piecewise polynomial with respect to $\mathcal{T}_{h}$ defined by $p_{\pi}=\Pi_{k-1}^{0,K}\,p$ for all $K\in\mathcal{T}_{h}$. Let us set

[TABLE]

From (4) and (45), it is straightforward to see that the couple $(\boldsymbol{\chi}_{h},\varrho_{h})$ solves the Darcy problem

[TABLE]

To prove the estimate for the pressure we employ the usual duality argument. Let therefore $\varphi$ be the solution of the auxiliary problem

[TABLE]

that, due to the convexity assumption, satisfies

[TABLE]

where the constant $C$ depends only on $\Omega$. For all $\mathbf{v}\in\widetilde{\mathbf{V}}_{h}$ let us denote with $\mathbf{v}_{I}$ its interpolant defined in (56) and (57). Therefore Green formula together with (61) yields

[TABLE]

We analyse separately the three terms. For the first one, using the consistency property of $\widetilde{a}_{h}(\cdot,\cdot)$, the polynomial approximation of $\mathbf{u}$ and $\nabla\varphi$, the estimate on the velocity error and (62) we get

[TABLE]

For what concerns the second term we have

[TABLE]

Finally, for the third term we begin by observing that from (61)

[TABLE]

for all $\varphi_{\pi}\in Q_{h}$, and by the Green formula

[TABLE]

Therefore, by collecting (66), (67), and using the previous error estimate, it holds that

[TABLE]

Finally by collecting (64), (65) and (68) in (63) we get the thesis. ∎