An H^{1}-conforming Virtual Element Methods for Darcy equations and Brinkman equations

# An H1-conforming Virtual Element Methods for Darcy equations and Brinkman equations

Giuseppe Vacca Dipartimento di Matematica e Applicazioni, Università degli Studi di Milano Bicocca, Via Roberto Cozzi 55 - 20125 Milano, Italy; E-mail: giuseppe.vacca@unimib.it.
###### Abstract

The focus of the present paper is on developing a Virtual Element Method for Darcy and Brinkman equations. In [stokes] 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 -projection onto 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 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.

## 1 Introduction

The Virtual Element Methods (in short, VEM or VEMs) is a recent technique for solving PDEs. VEMs were recently introduced in [VEM-volley] 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 [BLS05bis, BLM11book, Bishop13, JA12, LMSXX, RW12, POLY37, ST04, TPPM10, VW13, Wachspress11, DiPietro-Ern-1, di2016discontinuous, Gillette-1, PolyDG-1] 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 [lopezvacca, da2015symplectic], 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 [VEM-volley, Brezzi:Marini:plates, VEM-elasticity, VEM-enhanced, BM13, BFMXX, VEM-stream, GTP14, BBPSXX, benedetto2016hybrid, mora2015virtual, hpvem, lovadina, caceres2015mixed, wriggers, plates-zhao, vaccabis, vaccahyper, giani, gardini, frittelli]. A helpful paper for the implementation of the method is [VEM-hitchhikers].

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 [BFMXX, da2016mixed].

In [stokes] 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 [VEM-elasticity] 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 orthogonal projection onto the space of polynomials of degree less or equal than (where 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 [stokes], 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 [mardal2002robust, xie2008uniformly, karper2009unified, vassilev].

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 [stenberg, anaya2015priori, neilan]). 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 [VEM-elasticity]. 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 be a bounded polygon then the Darcy equation in mixed form is

 ⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩find (u,p) such thatK−1u+∇p=0in Ω,divu=fin Ω,u⋅n=0on ∂Ω, (1)

where and are respectively the velocity and the pressure fields, is the source term and 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 , the external force has zero mean value on . We consider the spaces

 V:={u∈H(div,Ω),s.tu⋅n=0on ∂Ω},Q:=L20(Ω)={q∈L2(Ω)s.t.∫ΩqdΩ=0}

equipped with the natural norms

 ∥v∥2V:=∥v∥2[L2(Ω)]2+∥divv∥2L2(Ω),∥q∥Q:=∥q∥L2(Ω),

and the bilinear forms and defined by:

 (2)
 b(v,q):=∫ΩdivvqdΩfor all v∈V, q∈Q. (3)

Then the variational formulation of Problem (1) is

 (4)

where

 (f,q):=∫ΩfqdΩfor all q∈Q.

Let us introduce the kernel

 Z:={v∈Vs.t.b(v,q)=0for all q∈Q};

then it is straightforward to see that

 ∥v∥V:=∥v∥[L2(Ω)]2for all v∈Z.

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

• and are continuous, i.e.

 |a(u,v)|≤∥a∥∥u∥V∥v∥Vfor all u,v∈V,
 |b(v,q)|≤∥b∥∥v∥V∥q∥Qfor % all v∈V and q∈Q;
• is coercive on the kernel , i.e. there exists a positive constant depending on such that

 a(v,v)≥α∥v∥2V% for all v∈Z; (5)
• satisfies the inf-sup condition, i.e.

 ∃β>0such thatsupv∈Vv≠0b(u,q)∥v∥V≥β∥q∥Qfor all q∈Q. (6)

Therefore, Problem (4) has a unique solution such that

 ∥u∥V+∥p∥Q≤C∥f∥L2(Ω)

with the constant depending only on and .

## 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 will indicate a generic positive constant independent of the mesh size that may change at each occurrence. Moreover, given any subset in and , we will denote by the polynomials of total degree at most defined on , with the extended notation . Let be a sequence of decompositions of into general polygonal elements with

 hK:=diameter(K),h:=supK∈ThhK.

We suppose that for all , each element in fulfils the following assumptions:

• is star-shaped with respect to a ball of radius ,

• the distance between any two vertexes of is ,

where and are positive constants. We remark that the hypotheses above, though not too restrictive in many practical cases, can be further relaxed, as noted in  [VEM-volley, 2016stability]. From now on we assume that is piecewise constant with respect to on .

Using standard VEM notation, for , let us define the spaces

• the set of polynomials on of degree ,

• ,

• ,

• with .

In [stokes] 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 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 [stokes] the notations, the main properties of the Virtual spaces and some details of the construction of the semi-norm projection. Let the polynomial degree of accuracy of the method, then we define on each element the finite dimensional local virtual space

 (7)

where all the operators and equations above are to be interpreted in the distributional sense. It is easy to check that , and that (see [stokes] for the proof) the dimension of is

 dim(WKh)=dim([Bk(∂K)]2)+dim(Gk−2(K)⊥)+(dim(Pk−1(K))−1)=2nKk+(k−1)(k−2)2+(k+1)k2−1. (8)

The corresponding degrees of freedom are chosen prescribing, given a function , the following linear operators , split into four subsets (see Figure 1):

• : the values of at the vertices of the polygon ,

• : the values of at distinct points of every edge (for example we can take the internal points of the -Gauss-Lobatto quadrature rule in , as suggested in [VEM-hitchhikers]),

• : the moments of

 ∫Kv⋅g⊥k−2dKfor all g⊥k−2∈Gk−2(K)⊥,
• : the moments up to order and greater than zero of in , i.e.

 ∫K(divv)qk−1dKfor all qk−1∈Pk−1(K)/R.

For all , we introduce the semi-norm projection , defined by

 ⎧⎪⎨⎪⎩∫K∇qk:∇(vh−Π∇,Kkvh)dK=0for all qk∈[Pk(K)]2,Π0,K0(vh−Π∇,Kkvh)=0, (9)

where is the -projection operator onto the constant functions defined on . It is immediate to check that the energy projection is well defined and it clearly holds for all . Moreover the operator is computable in terms of the degrees of freedom (see equations in [stokes] and the subsequent discussion).

### 3.2 The modified virtual space and the projection Π0,Kk

Let a positive integer, then for all , the -projection is defined by

 ∫Kqn⋅(vh−Π0,Knvh)dK=0for all qn∈[Pn(K)]2.

It is possible to check (see Section 3.3 of [stokes] for the proof) that the degrees of freedom allow us to compute exactly the -projection . On the other hand we can observe that we can not compute exactly from the DoFs the -projection onto the space of polynomials of degree . The goal of the present section is to introduce, taking the inspiration from [VEM-enhanced], a new virtual space to be used in place of in such a way that

• the DoFs can still be used for ,

• ,

• the projection can be exactly computable by the DoFs .

To construct we proceed as follows: first of all we define an augmented virtual local space by taking

Now we define the enhanced Virtual Element space as the restriction of given by

 VKh:={v∈UKhs.t.(v−Π∇,Kkv,g⊥k)[L2(K)]2=0for all g⊥k∈Gk(K)⊥/Gk−2(K)⊥}, (10)

where the symbol denotes the polynomials in that are orthogonal to all polynomials of . We proceed by investigating the dimension and by choosing suitable DoFs of the virtual space . First of all we recall from [conforming] the following facts

 dim([Bk(∂K)]2)=2nKk,dim(Pk−1(K))=k(k+1)2,dim(Gk(K)⊥)=k(k+1)2 (11)

where is the number of edges of the polygon .

###### Lemma 3.1.

The dimension of is

 dim(UKh)=2nKk+k(k+1)2+(k+1)k2−1.

Moreover as DoFs for we can take the linear operators and plus the moments

 DU:∫Kv⋅g⊥kdKfor all g⊥k∈Gk(K)⊥/Gk−2(K)⊥.
###### Proof.

The proof is virtually identical to that given in [stokes] for and it is based (see for instance [fortin1991mixed]) on the fact that given

• a polynomial function ,

• a polynomial function ,

• a polynomial function satisfying the compatibility condition

 ∫KgdΩ=∫∂Kgb⋅nds,

there exists a unique pair such that

 v|∂K=gb,divv=g,−Δv−∇s=h. (12)

Moreover, since from [conforming], is an isomorphism, we can conclude that the map that associates a given compatible data set to the velocity field that solves (12) is an injective map. Then

 dim(UKh)=dim([Bk(∂K)]2)+dim(Gk(K)⊥)+(dim(Pk−1(K))−1)

and the thesis follows from (11). ∎

###### Proposition 3.1.

The dimension of is equal to that of that is, as in (8)

 dim(VKh)=2nKk+(k−1)(k−2)2+(k+1)k2−1. (13)

As DoFs in we can take .

###### Proof.

From (11) it is straightforward to check that

 dim(Gk(K)⊥/Gk−2(K)⊥)=dim(Gk(K)⊥)−dim(Gk−2(K)⊥)=2k−1.

Hence, neglecting the independence of the additional conditions in (10), it holds that

 dim(VKh)≥dim(UKh)−(2k−1)=2nKk+(k−1)(k−2)2+(k+1)k2−1=dim(WKh). (14)

We now observe that a function such that is identically zero. Indeed, from (9), it is immediate to check that in this case the would be zero, implying that all its moment are zero, in particular, since , all the moments of are also zero. Now, from Lemma 3.1, we have that is zero. Therefore, from (14), we obtain that the dimension of is actually the same of , and that the DoFs are unisolvent for . ∎

###### Proposition 3.2.

The degrees of freedom allow us to compute exactly the -projection , i.e. the moments

 ∫Kv⋅qkdK

for all and for all .

###### Proof.

Let us set

 qk=∇qk+1+g⊥k−2+g⊥k.

with , and . Therefore using the Green formula and since , we get

 ∫Kv⋅qkdK=∫Kv⋅(∇qk+1+g⊥k−2+g⊥k)dK=−∫Kdivvqk+1dK+∫Kv⋅g⊥k−2dK+∫KΠ∇,Kkv⋅g⊥kdK+∫∂Kqk+1v⋅nds.

Now, since is a polynomial of degree less or equal than we can reconstruct its value from and compute exactly the first term. The second term is computable from . The third term is computable from all the using the projection . Finally from and we can reconstruct on the boundary and so compute exactly the boundary term. ∎

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

 QKh:=Pk−1(K) (15)

having dimension

 dim(QKh)=dim(Pk−1(K))=(k+1)k2.

The corresponding degrees of freedom are chosen defining for each the following linear operators :

• : the moments up to order of , i.e.

 ∫Kqpk−1dKfor all pk−1∈Pk−1(K).

Finally we define the global virtual element spaces as

 Vh:={v∈[H1(Ω)]2s.tv⋅n=0on ∂Ωandv|K∈VKhfor all K∈Th} (16)

and

 Qh:={q∈L20(Ω)s.t.q|K∈QKhfor all K∈Th}, (17)

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

 dim(Vh)=nP((k+1)k2−1+(k−1)(k−2)2)+2(nV+(k−1)nE)+(nV,B+(k−1)nE,B)

and

 dim(Qh)=nP(k+1)k2−1,

where is the number of elements, , (resp., , ) is the number of internal edges and vertexes (resp., boundary edges and vertexes) in . As observed in [stokes], we remark that

 divVh⊆Qh. (18)
###### Remark 3.1.

By definition (16) it is clear that our discrete velocities field is -conforming, in particular we obtain continuous velocities, whereas the natural discretization is only - 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 and a discrete version of the bilinear forms and given in (2) and (3). For simplicity we assume that the tensor is piecewise constant with respect to the decomposition , i.e. is constant on each polygon . First of all we decompose into local contributions the bilinear forms and , the norms and by defining

 a(u,v)=:∑K∈ThaK(u,v)for all u,v∈V
 b(v,q)=:∑K∈ThbK(v,q)for % all v∈V and q∈Q,

and

 ∥v∥V=:⎛⎝∑K∈Th∥v∥2V,K⎞⎠1/2for all v∈V,∥q∥Q=:⎛⎝∑K∈Th∥q∥2Q,K⎞⎠1/2for all q∈Q.

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

 b(v,q)=∑K∈ThbK(v,q)=∑K∈Th∫KdivvqdKfor all % v∈Vh, q∈Qh, (19)

i.e. as noticed in [stokes] we do not introduce any approximation of the bilinear form. We notice that (19) is computable from the degrees of freedom , and , since is polynomial in each element . On the other hand, the bilinear form needs to be dealt with in a more careful way. First of all, by Proposition 3.2, we observe that for all and for all , the quantity

 aK(qk,v)=∫KK−1qk⋅vdK.

is exactly computable by the DoFs. However, for an arbitrary pair , the quantity is clearly not computable. In the standard procedure of VEM framework, we define a computable discrete local bilinear form

 aKh(⋅,⋅):VKh×VKh→R (20)

approximating the continuous form and satisfying the following properties:

• -consistency: for all and

 aKh(qk,vh)=aK(qk,vh); (21)
• stability: there exist two positive constants and , independent of and , such that, for all , it holds

 α∗aK(vh,vh)≤aKh(vh,vh)≤α∗aK(vh,vh). (22)

Let be a (symmetric) stabilizing bilinear form, satisfying

 c∗aK(vh,vh)≤RK(vh,vh)≤c∗aK(vh,vh)for % all vh∈Vh such that Π0,Kkvh=0 (23)

with and positive constants independent of and . Then, we can set

 aKh(uh,vh):=aK(Π0,Kkuh,Π0,Kkvh)+RK((I−Π0,Kk)uh,(I−Π0,Kk)vh) (24)

for all .

It is straightforward to check that Definition (9) and properties (23) imply the consistency and the stability of the bilinear form .

###### Remark 3.2.

In the construction of the stabilizing form with condition (23) we essentially require that the stabilizing term scales as . Following the standard VEM technique (cf. [VEM-volley, VEM-hitchhikers] for more details), denoting with , the vectors containing the values of the local degrees of freedom associated to , we set

 RK(uh,vh)=αK¯uTh¯vh,

where is a suitable positive constant that scales as . For example, in the numerical tests presented in Section 6, we have chosen as the mean value of the eigenvalues of the matrix stemming from the term in (24).

Finally we define the global approximated bilinear form by simply summing the local contributions:

 ah(uh,vh):=∑K∈ThaKh(uh,vh)for all uh,vh∈Vh. (25)

### 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:

 ⎧⎪⎨⎪⎩find (uh,ph)∈Vh×Qh, such thatah(uh,vh)+b(vh,ph)=0for all vh∈Vh,b(uh,qh)=(f,qh)for all qh∈Qh. (26)

We point out that the symmetry of together with (22) easily implies that is (uniformly) continuous with respect to the norm. Moreover, as observed in [stokes], introducing the discrete kernel:

 Zh:={vh∈Vhs.t.b(vh,qh)=0for all qh∈Qh},

it is immediate to check that

 Zh⊆Z.

Then the bilinear form is also uniformly coercive on the discrete kernel with respect to the norm. Moreover as a direct consequence of Proposition 4.3 in [stokes], we have the following stability result.

###### Proposition 3.3.

Given the discrete spaces and defined in (16) and (17), there exists a positive , independent of , such that:

 supvh∈Vhvh≠0b(vh,qh)∥vh∥V≥~β∥qh∥Qfor all qh∈Qh. (27)

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

 divVh=Qh.

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

###### Theorem 3.1.

Problem (26) has a unique solution , verifying the estimate

 ∥uh∥V+∥ph∥Q≤C∥f∥0.

## 4 Theoretical results

We begin by proving an approximation result for the virtual local space . First of all, let us recall a classical result by Brenner-Scott (see [MR2373954]).

###### Lemma 4.1.

Let , then for all with , there exists a polynomial function , such that

 ∥u−uπ∥0,K+hK|u−uπ|1,K≤Chs+1K|u|s+1,K. (28)

We have the following approximation results (for the proof see [navier-stokes]).

###### Proposition 4.1.

Let with . Under the assumption and on the decomposition , there exists such that

 ∥u−uint∥0+hK|u−uint|1,K≤Chs+1K|u|s+1,K.

where is a constant independent of .

For what concerns the pressures, from classic polynomial approximation theory [MR2373954], for it holds

 infqh∈Qh∥q−qh∥Q≤Chk|q|k. (29)

We are ready to state the following convergence theorem.

###### Theorem 4.1.

Let be the solution of problem (4) and be the solution of problem (26). Then it holds

 ∥u−uh∥0≤Chk+1|u|k+1,and∥u−uh∥V≤Chk|u|k+1,∥p−ph∥Q≤Chk(|u|k+1+|p|k).
###### Proof.

We begin by remarking that as a consequence of the inf-sup condition with classical arguments (see for instance Proposition 2.5 in [fortin1991mixed]), there exists such that

 Π0,Kk−1(divuI)=divuI=Π0,Kk−1(divu)for all K∈Th, (30) (31)

Let us set . From (30) and (26), we have that and thus . Now, using (5), (22), (26) and introducing the piecewise polynomial approximation (28) together with (21), we have

 α∗α∥δh∥20≤α∗a(δh,δh)≤ah(δh,δh)=ah(uI,δh)−ah(uh,δh)=ah(uI,δh)+b(δh,ph)=ah(uI,δh)=∑K∈ThaKh(uI,δh)=∑K∈Th(aKh(uI−uπ,δh)+aK(uπ,δh))=∑K∈Th(aKh(uI−uπ,δh)+aK(uπ−u,δh))−a(u,δh)=∑K∈Th(aKh(uI−uπ,δh)+aK(uπ−u,δh))+b(δh,p)=∑K∈Th(aKh(uI−uπ,δh)+aK(uπ−u,δh))≤C∑K∈Th(∥uI−uπ∥0,K+∥u−uπ∥0,K)∥δh∥0,K≤C(∥uI−uπ∥0+∥u−uπ∥0)∥δh∥0

then

 ∥δh∥0≤C∥uI−uπ∥0+∥u−uπ∥0.

The -estimate follows easily by the triangle inequality. It is also straightforward to see from (4) and (26) that

 b(u−uh,qh)=0for all qh∈Qh,

than we get for all and therefore

 ∥div(u−uh)∥0=∑K∈Th∥divu−Π0,Kk−1(divu)∥0,K≤Chk|divu|k≤Chk|u|k+1,

from which the estimate in the norm. We proceed by analysing the error on the pressure field. Let , then from the discrete inf-sup condition (27), we infer:

 ~β∥ph−qh∥Q≤supvh∈Vhvh≠0b(vh,ph−qh)∥vh∥V=supvh∈Vhvh≠0b(vh,ph−p)+b(vh,p−qh)∥vh∥