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

2.4 MB·

by  Giuseppe Vacca

An H1-conforming Virtual Element Methods for Darcy
equations and Brinkman equations
Giuseppe Vacca ∗
Abstract
ThefocusofthepresentpaperisondevelopingaVirtualElementMethodforDarcyand
Brinkmanequations. In[15]wepresentedafamilyofVirtualElementsforStokesequations
7
andwedefinedanewVirtualElementspaceofvelocitiessuchthattheassociateddiscrete
1
0 kernel is pointwise divergence-free. We use a slightly different Virtual Element space
2 having two fundamental properties: the L2-projection onto Pk is exactly computable on
the basis of the degrees of freedom, and the associated discrete kernel is still pointwise
b
divergence-free. TheresultingnumericalschemefortheDarcyequationhasoptimalorder
e
F of convergence and H1 conforming velocity solution. We can apply the same approach to
developarobustvirtualelementmethodfortheBrinkmanequationthatisstableforboth
3
the Stokes and Darcy limit case. We provide a rigorous error analysis of the method and

several numerical tests.
]
A
N 1 Introduction
.
h
t TheVirtualElementMethods(inshort,VEMorVEMs)isarecenttechniqueforsolving
a
PDEs. VEMs were recently introduced in [5] as a generalization of the finite element method
m
on polyhedral or polygonal meshes. In the numerical analysis and engineering literature there
[
hasbeenarecentgrowthofinterestindevelopingnumericalmethodsthatcanmakeuseofgen-

2 eral polygonal and polyhedral meshes, as opposed to more standard triangular/quadrilateral
v (tetrahedral/hexahedral)grids. Indeed,makinguseofpolygonalmeshesbringsfortharangeof
0 advantages, including for instance automatic hanging node treatment, more efficient approxi-
8
mationofgeometricdatafeatures,betterdomainmeshingcapabilities,moreefficientandeasier
6
adaptivity, more robustness to mesh deformation, and others. This interest in the literature is
7
0 also reflected in commercial codes, such as CD-Adapco, that have recently included polytopal
. meshes.
1
0 We refer to the recent papers and monographs [25, 11, 20, 21, 37, 42, 44, 43, 45, 50, 51, 29,
7 30,41,28]asabriefrepresentativesampleoftheincreasinglistoftechnologiesthatmakeuseof
1 polygonal/polyhedralmeshes. Wementionhereinparticularthepolygonalfiniteelements,that
:
v generalizefiniteelementstopolygons/polyhedronsbymakinguseofgeneralizednon-polynomial
i shape functions, and the mimetic discretisation schemes [38, 12], that combine ideas from the
X
finite difference and finite element methods.
r
The principal idea behind VEM is to use approximated discrete bilinear forms that require
a
only integration of polynomials on the (polytopal) element in order to be computed. The
resultingdiscretesolutionisconformingandtheaccuracygrantedbysuchdiscretebilinearforms
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].
∗Dipartimento di Matematica e Applicazioni, Università degli Studi di Milano Bicocca, Via Roberto Cozzi
55-20125Milano,Italy;E-mail: [email protected].
1
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 ex-
ploiting the flexibility of the Virtual Element construction in a new way. In particular, they
defineanewVirtualElementspaceofvelocitiescarefullydesignedtosolvetheStokesproblem.
Inconnectionwithasuitablepressurespace,thenewVirtualElementspaceleadstoanexactly
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
introducingaslightlydifferentvirtualspaceforthevelocitiessuchthatthelocalL2 orthogonal
projectionontothespaceofpolynomialsofdegreelessorequalthank(wherekisthepolynomial
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
BrinkmanequationsthatstemsfromtheaboveschemefortheDarcyproblem. Mathematically,
theBrinkmanproblemresemblesboththeStokesproblemforfluidflowandtheDarcyproblem
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
featuresandtheconvergenceanalysisoftheproblem. InSection5wedevelopastablenumerical
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
WeconsidertheclassicalDarcyequationthatdescribestheflowofafluidthroughaporous
medium. Let Ω⊆R2 be a bounded polygon then the Darcy equation in mixed form is
find (u,p) such that
K−1u+∇p=0 in Ω,
(1)
divu=f in Ω,
u·n=0 on ∂Ω,
where u and p are respectively the velocity and the pressure fields, f ∈ L2(Ω) is the source
term and 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 ∂Ω, the
external force f has zero mean value on Ω. We consider the spaces
(cid:26) Z (cid:27)
V:={u∈H(div,Ω), s.t u·n=0 on ∂Ω}, Q:=L2(Ω)= q ∈L2(Ω) s.t. qdΩ=0
0
Ω
2
equipped with the natural norms
kvk2 :=kvk2 +kdivvk2 , kqk :=kqk ,
V [L2(Ω)]2 L2(Ω) Q L2(Ω)
and the bilinear forms a(·,·): V×V→R and b(·,·): V×Q→R defined by:
Z
a(u,v):= K−1u·vdΩ, for all u,v∈V (2)
Ω
Z
b(v,q):= divvqdΩ for all v∈V, q ∈Q. (3)
Ω
Then the variational formulation of Problem (1) is

find (u,p)∈V×Q, such that

a(u,v)+b(v,p)=0 for all v∈V, (4)
b(u,q)=(f,q) for all q ∈Q,
where
Z
(f,q):= fqdΩ for all q ∈Q.
Ω
Let us introduce the kernel
Z:={v∈V s.t. b(v,q)=0 for all q ∈Q};
then it is straightforward to see that
kvk :=kvk for all v∈Z.
V [L2(Ω)]2
It is well known that (see for instance [24]):
• a(·,·) and b(·,·) are continuous, i.e.
|a(u,v)|≤kakkuk kvk for all u,v∈V,
V V
|b(v,q)|≤kbkkvk kqk for all v∈V and q ∈Q;
V Q
• a(·,·) is coercive on the kernel Z, i.e. there exists a positive constant α depending on K
such that
a(v,v)≥αkvk2 for all v∈Z; (5)
V
• b(·,·) satisfies the inf-sup condition, i.e.
b(u,q)
∃β >0 such that sup ≥βkqk for all q ∈Q. (6)
kvk Q
v∈Vv6=0 V
Therefore, Problem (4) has a unique solution (u,p)∈V×Q such that
kuk +kpk ≤Ckfk
V Q L2(Ω)
with the constant C depending only on Ω and 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 ω in R2 and k ∈N, we will denote
by P (ω) the polynomials of total degree at most k defined on ω, with the extended notation
k
P (ω)=∅. Let {T } be a sequence of decompositions of Ω into general polygonal elements
−1 h h
K with
h :=diameter(K), h:= sup h .
K K
K∈Th
We suppose that for all h, each element K in T fulfils the following assumptions:
h
3
• (A1) K is star-shaped with respect to a ball of radius ≥ γh ,
K
• (A2) the distance between any two vertexes of K is ≥ch ,
K
where γ 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 K is piecewise constant with respect to T on Ω.
h
Using standard VEM notation, for k ∈N, let us define the spaces
• P (K) the set of polynomials on K of degree ≤k,
k
• B (K):={v ∈C0(∂K) s.t v ∈P (e) ∀ edge e⊂∂K},
k |e k
• G (K):=∇(P (K))⊆[P (K)]2,
k k+1 k
• G (K)⊥ :=x⊥[P (K)]⊆[P (K)]2 with x⊥ :=(x ,−x ).
k k−1 k 2 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 H1 semi-norm projection onto the
space of polynomials can be exactly computed,
• thechoiceoftheVirtualspaceofvelocitiesandtheassociateddegreesoffreedomguarantee
thatthefinaldiscretevelocityispointwisedivergence-freeandmoregenerallythediscrete
kernel is contained in the continuous one.
Inthissectionwebrieflyrecallfrom[15]thenotations,themainpropertiesoftheVirtualspaces
andsomedetailsoftheconstructionoftheH1 semi-normprojection. Letk ≥2thepolynomial
degreeofaccuracyofthemethod,thenwedefineoneachelementK ∈T thefinitedimensional
h
local virtual space
(cid:26)
WK := v∈[H1(K)]2 s.t v ∈[B (∂K)]2,
h |∂K k
(cid:26) −∆v−∇s∈G (K)⊥, )
k−2 for some s∈L2(K) (7)
divv∈P (K),
k−1
where all the operators and equations above are to be interpreted in the distributional sense.
Itiseasytocheckthat[P (K)]2 ⊆WK,andthat(see[15]fortheproof)thedimensionofWK
k h h
is
dim(cid:0)WK(cid:1)=dim(cid:0)[B (∂K)]2(cid:1)+dim(cid:0)G (K)⊥(cid:1)+(dim(P (K))−1)
h k k−2 k−1
(k−1)(k−2) (k+1)k (8)
=2n k+ + −1.
K 2 2
The corresponding degrees of freedom are chosen prescribing, given a function v ∈ WK, the
h
following linear operators D , split into four subsets (see Figure 1):
V
• D 1: the values of v at the vertices of the polygon K,
V
• D 2: the values of v at k −1 distinct points of every edge e ∈ ∂K (for example we
V
can take the k−1 internal points of the (k+1)-Gauss-Lobatto quadrature rule in e, as
suggested in [7]),
• D 3: the moments of v
V
Z
v·g⊥ dK for all g⊥ ∈G (K)⊥,
k−2 k−2 k−2
K
4
• D 4: the moments up to order k−1 and greater than zero of divv in K, i.e.
V
Z
(divv)q dK for all q ∈P (K)/R.
k−1 k−1 k−1
K
Figure 1: Degrees of freedom for k =2, k =3. We denote D 1 with the black dots, D 2 with
V V
the red squares, D 3 with the green rectangles, D 4 with the blue dots inside the element.
V V
ForallK ∈T ,weintroducetheH1 semi-normprojectionΠ∇,K: WK →[P (K)]2,defined
h k h k
by
 Z
 ∇qk :∇(vh− Π∇k,Kvh)dK =0 for all qk ∈[Pk(K)]2,
K (9)
Π0,K(v − Π∇,Kv )=0,
0 h k h
where Π0,K is the L2-projection operator onto the constant functions defined on K. It is
0
immediatetocheckthattheenergyprojectioniswelldefinedanditclearlyholdsΠ∇,Kq =q
k k k
for all q ∈ P (K). Moreover the operator Π∇,K is computable in terms of the degrees of
k k k
freedom D (see equations (27)−(29) in [15] and the subsequent discussion).
V
3.2 The modified virtual space and the projection Π0,K
k
Let n a positive integer, then for all K ∈ T , the L2-projection Π0,K: WK → [P (K)]2 is
h n h n
defined by
Z
q ·(v −Π0,Kv )dK =0 for all q ∈[P (K)]2.
n h n h n n
K
It is possible to check (see Section 3.3 of [15] for the proof) that the degrees of freedom D
V
allow us to compute exactly the L2-projection Π0,K. On the other hand we can observe that
k−2
we can not compute exactly from the DoFs the L2-projection onto the space of polynomials of
degree ≤ k. The goal of the present section is to introduce, taking the inspiration from [1], a
new virtual space VK to be used in place of WK in such a way that
h h
• the DoFs D can still be used for VK,
V h
• [P (K)]2 ⊆VK,
k h
• the projection Π0,K: VK →[P (K)]2 can be exactly computable by the DoFs D .
k h k V
To construct VK we proceed as follows: first of all we define an augmented virtual local
h
space UK by taking
h
(cid:26)
UK := v∈[H1(K)]2 s.t v ∈[B (∂K)]2,
h |∂K k
(cid:26) −∆v−∇s∈G (K)⊥, )
k for some s∈L2(K)
divv∈P (K),
k−1
5
Now we define the enhanced Virtual Element space VK as the restriction of UK given by
h h
(cid:26) (cid:16) (cid:17) (cid:27)
VK := v∈UK s.t. v−Π∇,Kv, g⊥ =0 for all g⊥ ∈G (K)⊥/G (K)⊥ ,
h h k k k k k−2
[L2(K)]2
(10)
wherethesymbolG (K)⊥/G (K)⊥denotesthepolynomialsinG (K)⊥thatareL2−orthogonal
k k−2 k
to all polynomials of G (K)⊥. We proceed by investigating the dimension and by choosing
k−2
suitable DoFs of the virtual space V . First of all we recall from [9] the following facts
h
dim(cid:0)[B (∂K)]2(cid:1)=2n k, dim(P (K))= k(k+1), dim(cid:0)G (K)⊥(cid:1)= k(k+1) (11)
k K k−1 2 k 2
where n is the number of edges of the polygon K.
K
Lemma 3.1. The dimension of UK is
h
dim(cid:0)UK(cid:1)=2n k+ k(k+1) + (k+1)k −1.
h K 2 2
Moreover as DoFs for UK we can take the linear operators D and plus the moments
h V
Z
D : v·g⊥dK for all g⊥ ∈G (K)⊥/G (K)⊥.
U k k k k−2
K
Proof. The proof is virtually identical to that given in [15] for WK and it is based (see for
h
instance [24]) on the fact that given
• a polynomial function g ∈[B (∂K)]2,
b k
• a polynomial function h∈G (K)⊥,
k
• a polynomial function g ∈P (K) satisfying the compatibility condition
k−1
Z Z
gdΩ= g ·nds,
b
K ∂K
there exists a unique pair (v,s)∈UK ×L2(K)/R such that
h
v =g , divv=g, −∆v−∇s=h. (12)
|∂K b
Moreover,sincefrom[9],rot: G (K)⊥ →P (K)isanisomorphism,wecanconcludethatthe
k k−1
mapthatassociatesagivencompatibledataset(g , h, g)tothevelocityfieldvthatsolves(12)
b
is an injective map. Then
dim(cid:0)UK(cid:1)=dim(cid:0)[B (∂K)]2(cid:1)+dim(cid:0)G (K)⊥(cid:1)+(dim(P (K))−1)
h k k k−1
and the thesis follows from (11).
Proposition 3.1. The dimension of VK is equal to that of WK that is, as in (8)
h h
dim(cid:0)VK(cid:1)=2n k+ (k−1)(k−2) + (k+1)k −1. (13)
h K 2 2
As DoFs in VK we can take D .
h V
Proof. From (11) it is straightforward to check that
dim(cid:0)G (K)⊥/G (K)⊥(cid:1)=dim(cid:0)G (K)⊥(cid:1)−dim(cid:0)G (K)⊥(cid:1)=2k−1.
k k−2 k k−2
Hence, neglecting the independence of the additional 2k−1 conditions in (10), it holds that
dim(cid:0)VK(cid:1)≥dim(cid:0)UK(cid:1)−(2k−1)=2n k+(k−1)(k−2)+(k+1)k−1=dim(cid:0)WK(cid:1). (14)
h h K 2 2 h
We now observe that a function v∈VK such that D (v)=0 is identically zero. Indeed, from
h V
(9), it is immediate to check that in this case the Π∇,Kv would be zero, implying that all its
k
moment are zero, in particular, since v ∈ VK, all the moments D of v are also zero. Now,
h U
from Lemma 3.1, we have that v is zero. Therefore, from (14), we obtain that the dimension
of VK is actually the same of WK, and that the DoFs D are unisolvent for VK.
h h V h
6
Proposition 3.2. The degrees of freedom D allow us to compute exactly the L2-projection
V
Π0,K: V →[P (K)]2, i.e. the moments
k h k
Z
v·q dK
k
K
for all v∈V and for all q ∈[P (K)]2.
h k k
Proof. Let us set
q =∇q +g⊥ +g⊥.
k k+1 k−2 k
with q ∈ P (K)/R, g⊥ ∈ G⊥ (K) and g⊥ ∈ G⊥(K)/G⊥ (K). Therefore using the
k+1 k+1 k−2 k−2 k k k−2
Green formula and since v∈V , we get
h
Z Z
v·q dK = v·(∇q +g⊥ +g⊥)dK
k k+1 k−2 k
K K
Z Z Z Z
=− divvq dK+ v·g⊥ dK+ Π∇,Kv·g⊥dK+ q v·nds.
k+1 k−2 k k k+1
K K K ∂K
Now, since divv is a polynomial of degree less or equal than k−1 we can reconstruct its value
fromD 4andcomputeexactlythefirstterm. ThesecondtermiscomputablefromD 3. The
V V
third term is computable from all the D using the projection Π∇,Kv. Finally from D 1 and
V k V
D 2 we can reconstruct v on the boundary and so compute exactly the boundary term.
V
For what concerns the pressures we take the standard finite dimensional space
QK :=P (K) (15)
h k−1
having dimension
(k+1)k
dim(QK)=dim(P (K))= .
h k−1 2
The corresponding degrees of freedom are chosen defining for each q ∈QK the following linear
h
operators D :
Q
• D : the moments up to order k−1 of q, i.e.
Q
Z
qp dK for all p ∈P (K).
k−1 k−1 k−1
K
Finally we define the global virtual element spaces as
V :={v∈[H1(Ω)]2 s.t v·n=0 on ∂Ω and v ∈VK for all K ∈T } (16)
h |K h h
and
Q :={q ∈L2(Ω) s.t. q ∈QK for all K ∈T }, (17)
h 0 |K h h
with the obvious associated sets of global degrees of freedom. A simple computation shows
that:
(cid:18) (cid:19)
(k+1)k (k−1)(k−2)
dim(V )=n −1+ +2(n +(k−1)n )+(n +(k−1)n )
h P 2 2 V E V,B E,B
and
(k+1)k
dim(Q )=n −1,
h P 2
wheren isthenumberofelements,n ,n (resp.,n ,n )isthenumberofinternaledges
P E V E,B V,B
and vertexes (resp., boundary edges and vertexes) in T . As observed in [15], we remark that
h
divV ⊆Q . (18)
h h
Remark 3.1. By definition (16) it is clear that our discrete velocities field is H1-conforming, in
particular we obtain continuous velocities, whereas the natural discretization is only H(div)-
conforming. This property, in combination with (18), will make our method suitable for a
(robust) extension to the Brinkman problem.
7
3.3 The discrete bilinear forms
The next step in the construction of our method is to define on the virtual spaces V and
h
Q a discrete version of the bilinear forms a(·,·) and b(·,·) given in (2) and (3). For simplicity
h
we assume that the tensor K is piecewise constant with respect to the decomposition T , i.e.
h
K is constant on each polygon K ∈ T . First of all we decompose into local contributions the
h
bilinear forms a(·,·) and b(·,·), the norms k·k and k·k by defining
V Q
X
a(u,v)=: aK(u,v) for all u,v∈V
K∈Th
X
b(v,q)=: bK(v,q) for all v∈V and q ∈Q,
K∈Th
and
!1/2 !1/2
X X
kvk =: kvk2 for all v∈V, kqk =: kqk2 for all q ∈Q.
V V,K Q Q,K
K∈Th K∈Th
We now define discrete versions of the bilinear form a(·,·) (cf. (2)), and of the bilinear form
b(·,·) (cf. (3)). For what concerns b(·,·), we simply set
Z
X X
b(v,q)= bK(v,q)= divvqdK for all v∈V , q ∈Q , (19)
h h
K∈Th K∈Th K
i.e. as noticed in [15] we do not introduce any approximation of the bilinear form. We notice
that(19)iscomputablefromthedegreesoffreedomD 1,D 2andD 4,sinceqispolynomial
V V V
in each element K ∈T . On the other hand, the bilinear form a(·,·) needs to be dealt with in
h
a more careful way. First of all, by Proposition 3.2, we observe that for all q ∈[P (K)]2 and
k k
for all v∈VK, the quantity
h
Z
aK(q ,v)= K−1q ·vdK.
k k
K
is exactly computable by the DoFs. However, for an arbitrary pair (u,v) ∈ VK ×VK, the
h h
quantity aK(w,v) is clearly not computable. In the standard procedure of VEM framework,
h
we define a computable discrete local bilinear form
aK(·,·): VK ×VK →R (20)
h h h
approximating the continuous form aK(·,·) and satisfying the following properties:
• k-consistency: for all q ∈[P (K)]2 and v ∈VK
k k h h
aK(q ,v )=aK(q ,v ); (21)
h k h k h
• stability: there exist two positive constants α and α∗, independent of h and K, such
∗
that, for all v ∈VK, it holds
h h
α aK(v ,v )≤aK(v ,v )≤α∗aK(v ,v ). (22)
∗ h h h h h h h
Let RK: VK ×VK →R be a (symmetric) stabilizing bilinear form, satisfying
h h
c aK(v ,v )≤RK(v ,v )≤c∗aK(v ,v ) for all v ∈V such that Π0,Kv =0 (23)
∗ h h h h h h h h k h
with c and c∗ positive constants independent of h and K. Then, we can set
∗
(cid:16) (cid:17) (cid:16) (cid:17)
aK(u ,v ):=aK Π0,Ku ,Π0,Kv +RK (I−Π0,K)u ,(I−Π0,K)v (24)
h h h k h k h k h k h
for all u ,v ∈VK.
h h h
It is straightforward to check that Definition (9) and properties (23) imply the consistency
and the stability of the bilinear form aK(·,·).
h
8
Remark 3.2. In the construction of the stabilizing form RK with condition (23) we essentially
requirethatthestabilizingtermRK(v ,v )scalesasaK(v ,v ). FollowingthestandardVEM
h h h h
technique (cf. [5, 7] for more details), denoting with u¯h, v¯h ∈RNK the vectors containing the
values of the N local degrees of freedom associated to u ,v ∈VK, we set
K h h h
RK(u ,v )=αKu¯Tv¯ ,
h h h h
whereαK isasuitablepositiveconstantthatscalesas|K|. Forexample, inthenumericaltests
presented in Section 6, we have chosen αK as the mean value of the eigenvalues of the matrix
(cid:16) (cid:17)
stemming from the term aK Π0,Ku , Π0,Kv in (24).
k h k h
Finally we define the global approximated bilinear form a (·,·): V ×V → R by simply
h h h
summing the local contributions:
X
a (u ,v ):= aK(u ,v ) for all u ,v ∈V . (25)
h h h h h h h h h
K∈Th
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 (u ,p )∈V ×Q , such that
 h h h h
a (u ,v )+b(v ,p )=0 for all v ∈V , (26)
h h h h h h h
b(u ,q )=(f,q ) for all q ∈Q .
h h h h h
We point out that the symmetry of a (·,·) together with (22) easily implies that a (·,·) is
h h
(uniformly)continuouswithrespecttotheL2 norm. Moreover,asobservedin[15],introducing
the discrete kernel:
Z :={v ∈V s.t. b(v ,q )=0 for all q ∈Q },
h h h h h h h
it is immediate to check that
Z ⊆Z.
h
Thenthebilinearforma (·,·)isalsouniformlycoerciveonthediscretekernelZ withrespectto
h h
the 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 V and Q defined in (16) and (17), there exists
h h
a positive β˜, independent of h, such that:
b(v ,q )
sup h h ≥β˜kq k for all q ∈Q . (27)
kv k h Q h h
vh∈Vhvh6=0 h V
Inparticular, thetheinf-supconditionofProposition3.3, alongwithproperty (18), implies
that:
divV =Q .
h h
Finally we can state the well-posedness of virtual problem (26).
Theorem 3.1. Problem (26) has a unique solution (u ,p )∈V ×Q , verifying the estimate
h h h h
ku k +kp k ≤Ckfk .
h V h Q 0
9
4 Theoretical results
WebeginbyprovinganapproximationresultforthevirtuallocalspaceV . Firstofall, let
h
us recall a classical result by Brenner-Scott (see [22]).
Lemma4.1. LetK ∈T ,thenforallu∈[Hs+1(K)]2 with0≤s≤k,thereexistsapolynomial
h
function u ∈[P (K)]2, such that
π k
ku−u k +h |u−u | ≤Chs+1|u| . (28)
π 0,K K π 1,K K s+1,K
We have the following approximation results (for the proof see [16]).
Proposition 4.1. Let u ∈ V∩[Hs+1(Ω)]2 with 0 ≤ s ≤ k. Under the assumption (A1) and
(A2) on the decomposition T , there exists u ∈W such that
h int h
ku−u k +h |u−u | ≤Chs+1|u| .
int 0 K int 1,K K s+1,K
where C is a constant independent of h.
For what concerns the pressures, from classic polynomial approximation theory [22], for
q ∈Hk(Ω) it holds
inf kq−q k ≤Chk|q| . (29)
h Q k
qh∈Qh
We are ready to state the following convergence theorem.
Theorem 4.1. Let (u,p) ∈ V×Q be the solution of problem (4) and (u ,p ) ∈ V ×Q be
h h h h
the solution of problem (26). Then it holds
ku−u k ≤Chk+1|u| , and ku−u k ≤Chk|u| ,
h 0 k+1 h V k+1
kp−p k ≤Chk(|u| +|p| ).
h Q k+1 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 [24]), there exists u ∈V such that
I h
Π0,K(divu )=divu =Π0,K(divu) for all K ∈T , (30)
k−1 I I k−1 h
ku−u k ≤C inf ku−vk and ku−u k ≤C inf ku−vk . (31)
I 0 0 I V V
vh∈Vh vh∈Vh
Let us set δ =u −u . From (30) and (26), we have that divδ =0 and thus δ ∈Z . Now,
h I h h h h
using (5), (22), (26) and introducing the piecewise polynomial approximation (28) together
with (21), we have
α αkδ k2 ≤α a(δ , δ )≤a (δ , δ )=a (u , δ )−a (u , δ )
∗ h 0 ∗ h h h h h h I h h h h
=a (u , δ )+b(δ ,p )=a (u , δ )
h I h h h h I h
= X aK(u , δ )= X (cid:0)aK(u −u , δ )+aK(u ,δ )(cid:1)
h I h h I π h π h
K∈Th K∈Th
= X (cid:0)aK(u −u , δ )+aK(u −u, δ )(cid:1)−a(u, δ )
h I π h π h h
K∈Th
= X (cid:0)aK(u −u , δ )+aK(u −u, δ )(cid:1)+b(δ ,p)
h I π h π h h
K∈Th
= X (cid:0)aK(u −u , δ )+aK(u −u, δ )(cid:1)
h I π h π h
K∈Th
X
≤C (ku −u k +ku−u k )kδ k
I π 0,K π 0,K h 0,K
K∈Th
≤C (ku −u k +ku−u k )kδ k
I π 0 π 0 h 0
then
kδ k ≤Cku −u k +ku−u k .
h 0 I π 0 π 0
10