Virtual Elements for the Navier-Stokes problem on polygonal meshes

A family of Virtual Element Methods for the 2D Navier-Stokes equations is proposed and analysed. The schemes provide a discrete velocity field which is point-wise divergence-free. A rigorous error analysis is developed, showing that the methods are stable and optimally convergent. Several numerical tests are presented, confirming the theoretical predictions. A comparison with some mixed finite elements is also performed.


Introduction
The Virtual Element Method (VEM), introduced in [11,12], is a recent paradigm for the approximation of partial differential equation problems that shares the same variational background of the Finite Element Methods.The original motivation of VEM is the need to construct an accurate conforming Galerkin scheme with the capability to deal with highly general polygonal/polyhedral meshes, including "hanging vertexes" and non-convex shapes.Among the Galerkin schemes, VEM has the peculiarity that the discrete spaces consist of functions which are not known pointwise, but a limited set of information about them are at disposal.Nevertheless, the available information are sufficient to construct the stiffness matrix and the right-hand side.
In this paper, which may be considered as a natural evolution of our recent divergence-free approach developed in [15] for the Stokes problem, we apply the VEM to the Navier-Stokes equations in 2D.However, the non-linear convective term in the Navier-Stokes equations leads to the introduction of suitable projectors.These, in turn, suggest to make use of an enhanced discrete velocity space [43], that is an improvement with respect to that of [15].Instead, the pressure field is approximated by means of standard locally polynomial functions, without any continuity requirement across the elements.Furthermore, we consider two different discretization of the trilinear form arising from the convective term.The first one is the straightforward VEM version of the continuous trilinear form; however, the projector introduction causes a lack of skew-symmetry, despite the discrete velocity is divergence-free (up to machine precision).This leads to consider the second choice, which is simply the skew-symmetric part of the trilinear form mentioned above (cf.[32], for instance).We remark that we develop an error analysis focusing on this latter choice, but the numerical tests concern both alternatives.The outcome is a family of Virtual Elements, one per each polynomial order of consistency k, with k ≥ 2. To our best knowledge, this is the first paper where the VEM technology is addressed to the Navier-Stokes equations.
The main objectives of the present paper are the following.
• The development of a rigorous error analysis of the proposed methods.We highlight that our analysis provides some noteworthy element of novelty.Indeed, although we follow rather well-established lines for the error analysis of 2D Navier-Stokes Galerkin methods (see for example [32]), these need to be combined with new techniques that are peculiar to the VEM framework.In particular, the interpolant construction of Theorem 4.1 involves new arguments which might be useful even in different contexts (i.e. for other VEM spaces with different regularity requirement).
• A first but thorough assessment of the actual numerical performance of this new approach.
We provide a set of significant numerical tests, that highlight the features of our VEM approach.In addition to the important flexibility of dealing with general polygonal meshes, the presented scheme (we tested the case k = 2) displays the following favourable points.
1.The error components partly decouple: notably, the velocity error does not depend directly on the discrete pressures, but only indirectly through the approximation of the loading and convection terms.This is a consequence of the fact that our methods provide a discrete velocity which is point-wise divergence-free (the isochoric constraint is not relaxed).In some situations, e.g. for hydrostatic fluid problems, the partial decoupling of the errors induces a positive effect on the velocity approximation.Moreover, for the same reason, the VEM scheme seems to be more robust for small values of the viscosity parameter when compared with standard mixed finite elements.
2. Another advantage of the method is that, again due to its divergence-free nature, the same Virtual space couple can be used directly also for the approximation of the diffusion problem (in mixed form).This allows for a much easier coupling in Stokes-Darcy problems where different models need to be used in different parts of the domain.This observation adds up with the fact that, thanks to the use of polygons that allow hanging nodes, the gluing of different meshes in different parts of the domain is also much easier.
3. As in [15], the particular choice of degrees of freedom adopted for the velocity space yields a diagonal structure in a large part of the pressure-velocity interaction stiffness matrix.As a consequence, and without the need of any static condensation, many internal-to-element degrees of freedom can be automatically ignored when building the linear system.
We finally note that, nowadays, there do exist Galerkin-type finite element methods for the Stokes and Navier-Stokes equations that are pressure-robust (that is, the error on the velocity does not depend on the pressure, not even indirectly through the loading or convection terms).Some recent examples are [33,28].However, to our best knowledge, all the available schemes work only for standard simplicial/hexahedral meshes.Despite our method is not pressurerobust in the sense above, for arbitrary polygonal meshes it is the only conforming divergencefree scheme, a property which yields to important advantages, as outlined in points 1 and 2. Developing a conforming scheme which is both divergence-free and pressure-robust for general polygonal meshes, is currently an open problem.
A brief outline of the paper is the following.In Section 2 we recall the 2D Navier-Stokes problem, introducing the classical variational formulation and the necessary notations.Section 3 details the proposed discretization procedure.The approximation spaces and all the quantities that form the discrete problem, are introduced and described.Section 4 deals with the theoretical analysis, which leads to the optimal error estimates of Theorem 4.2 and bound (94).Finally, Section 5 presents several numerical tests, which highlight the actual performance of our approach, also in comparison with a couple of well-known mixed finite element schemes.

The continuous Navier-Stokes equation
We consider the steady Navier-Stokes Equation on a polygonal domain Ω ⊆ R 2 with homogeneous Dirichlet boundary conditions: with ν ∈ R, ν > 0, and where u, p are the velocity and the pressure fields, respectively.Furthermore, ∆, div, ∇, and ∇ denote the vector Laplacian, the divergence, the gradient operator for vector fields and the gradient operator for scalar functions.Finally, f represents the external force, while ν is the viscosity.We also remark that different boundary conditions can be treated as well.
Let us consider the spaces with norms We assume f ∈ [L 2 (Ω)] 2 and consider the bilinear forms Then a standard variational formulation of Problem (1) is: where It is well known that with the choices (3), we have (see for instance [32]): • the bilinear form b(•, •) and the space V and Q satisfy the inf-sup condition, i.e.
Let us introduce the kernel Then Problem ( 7) can be formulated in the equivalent kernel form Finally, by a direct computation it is easy to see that, if u ∈ Z is fixed, then the bilinear form c(u; Therefore we introduce, as usual, also the trilinear form c(•, 3 Virtual formulation of the problem

Virtual element space and polynomial projections
We outline the Virtual Element discretization of Problem (7).We will make use of various tools from the virtual element technology, that will be described briefly; we refer the interested reader to the papers [15,43].
Let { Ω h } h be a sequence of decompositions of Ω into general polygonal elements E with We suppose that for all h, each element E in Ω h fulfils the following assumptions: 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 investigated in [14].Using standard VEM notation, for k ∈ N, let us define the spaces the set of polynomials on E of degree ≤ k (with the extended notation ).For any n ∈ N and E ∈ Ω h we introduce the following useful polynomial projections: for all v ∈ V and for all • the L 2 -projection for scalar functions Π 0,E n : L 2 (Ω) → P n (E), given by with obvious extension for vector functions Π . In [15] we 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.In [43] we have presented an enhanced Virtual space, taking the inspiration from [1], to be used in place of the original one in such a way that the L 2 -projection can be exactly computable by the DoFs.In this section we briefly recall from [15,43] the notations, the main properties of the Virtual spaces and some details about the construction of the projections.
Let k ≥ 2 the polynomial degree of accuracy of the method.We introduce on each element E ∈ Ω h the (original) finite dimensional local virtual space [15] for some s ∈ L 2 (E) (16) where all the operators and equations above are to be interpreted in the distributional sense.
Then we enlarge the previous space Now we define the Virtual Element space V E h as the restriction of U E h given by where the symbol . From [9, 15, 43], we recall the following properties of the space V E h .The proof of the following result can be found in [43].Proposition 3.1 (Dimension and DoFs).Let V E h be the space defined in (17).Then the dimension and where n E is the number of vertexes of E. Moreover the following linear operators D V , split into four subsets (see Figure 1) constitute a set of DoFs for V E h : • D V 1: the values of v at the vertices of the polygon E, • D V 2: the values of v at k − 1 distinct points of every edge e ∈ ∂E, The proof of the following result can be found in [15] for Π ∇,E k and in [43] for the remaining projectors.

Proposition 3.2 (Projections and Computability). The DoFs D V allow us to compute exactly
in the sense that, given any v h ∈ V E h , we are able to compute the polynomials k−1 ∇v h only using, as unique information, the degree of freedom values D V of v h .Remark 3.1.Using the enhanced space V E h and following the same ideas of [15,43], it is possible to improve the results of Proposition 3.2 and compute exactly also the following higher order projections Π ∇,E k+2 : Moreover, given any polynomial q n of arbitrary degree n and any v ∈ V E h , an integration by parts shows that we can compute the moment For what concerns the pressures we take the standard finite dimensional space The corresponding degrees of freedom are chosen defining for each q ∈ Q E h the following linear operators D Q : • D Q : the moments up to order k − 1 of q, i.e.
Finally we define the global virtual element spaces as and with the obvious associated sets of global degrees of freedom.A simple computation shows that: where n P is the number of elements, n e , n V is the number of internal edges and vertexes in Ω h .As observed in [15], we remark that

Discrete bilinear forms and load term approximation
The next step in the construction of our method is to define a discrete version of the bilinear forms a(•, •) and b(•, •) given in ( 4) and ( 5) and trilinear form c(•; •, •) in ( 6).Here and in the rest of the paper the symbol C will indicate a generic positive quantity, independent of the mesh size (and of ν), but may depend on Ω and on the polynomial degree k.Furthermore, C may vary at each occurrence.First of all we decompose into local contributions the bilinear forms and i.e. as noticed in [15] we do not introduce any approximation of the bilinear form.We notice that ( 23) is computable from the degrees of freedom D V 1, D V 2 and D V 4, since q is polynomial in each element E ∈ Ω h .We now define discrete versions of the forms a(•, •) (cf.( 4)) and c(•; •, •) (cf.( 6)) that need to be dealt with in a more careful way.First of all, we note that for an arbitrary pair (u, v) following a standard procedure in the VEM framework, we define a computable discrete local bilinear form approximating the continuous form a E (•, •), and defined by for all u h , v h ∈ V E h , where the (symmetric) stabilizing bilinear form with α * and α * positive constants independent of the element E. It is straightforward to check that Definition (14) and properties ( 26) imply • stability: there exist two positive constants α * and α * , independent of h and E, such that, for all Remark 3.2.Condition (26) essentially requires that the stabilizing term where α E is a suitable positive constant.For example, in the numerical tests presented in Section 5, we have chosen α E as the mean value of the non-zero eigenvalues of the matrix stemming from the term (25).
Finally we define the global approximated bilinear form a h (•, •) : V h × V h → R by simply summing the local contributions: For what concerns the approximation of the local trialinear form c E (•; •, •), we set and note that all quantities in the previous formula are computable, in the sense of Proposition 3.2.As usual we define the global approximated trilinear form by adding the local contributions: We first notice that the form c h (•; •, •) is immediately extendable to the whole V (simply apply the same definition for any w, u, v ∈ V) .Moreover, we now show that it is continuous on V, uniformly in h.

Proposition 3.3. Let
Then C h is uniformly bounded, i.e. the trilinear form c h (•; •, •) is uniformly continuous with respect to h.

Proof. By a direct computation it holds
where the last inequality follows by using Hölder inequality.Let us analyse each term in the right hand side of (33).Employing the continuity of the projection Π 0,E k with respect the L 2 -norm we easily get For what concerns the second term (and analogously for the third one) we get (inverse estimate for polynomials) 34) and ( 35) in (33) we obtain Now applying Hölder inequality (for sequences) we get Finally, since H 1 (Ω) ⊂ L 4 (Ω), by Sobolev embedding it holds where the constant C h does not depend on h.
We can also define the local discrete skew-symmetric trilinear form with obvious global extension that is (obviously) still continuous and computable.
The last step consists in constructing a computable approximation of the right-hand side (f , v) in (7).We define the approximated load term f h as and consider: We observe that (41) can be exactly computed from D V for all v h ∈ V h (see Proposition 3.2).

The discrete problem
We are now ready to state the proposed discrete problem.Referring to ( 20), ( 21), ( 29), ( 39) and ( 23), we consider the virtual element problem: We point out that the symmetry of a h (•, •) together with (28) easily implies that a h (•, •) is continuous and coercive with respect to the V-norm.Moreover, as a direct consequence of Proposition 4.3 in [15], we have the following stability result.Proposition 3.4.Given the discrete spaces V h and Q h defined in (20) and (21), there exists a positive β, independent of h, such that: In particular, the inf-sup condition of Proposition 3.4, along with property (22), implies that: div The well-posedness of virtual problem ( 42) is a consequence of the coercivity property of a h (•, •), the skew-symmetry of c h (•; •, •) and the inf-sup condition (43).We have Theorem 3.1.Assuming that Moreover, as observed in [15], introducing the discrete kernel Problem ( 42) can be also formulated in the equivalent kernel form Remark 3.3.An alternative choice for the discretization ( 42) is to substitute the skew-simmetric form c h (•; •, •) with c h (•; •, •).With that choice, a theoretical analysis can be developed using the guidelines in [34] in connection with the same tools and ideas of Section 4. We here prefer to consider the choice (42), that allows for a more direct stability argument.Nevertheless, in the numerical tests of Section 5 we will investigate both possibilities.
Remark 3.4.An additional interesting consequence of property ( 46) is that, following [15,43], the proposed virtual elements can accommodate both the Stokes (or Navier-Stokes) and the Darcy problems simultaneously.Indeed, due to property (46), the proposed velocity-pressure couple turns out to be stable not only for the Stokes problem, but also for the Darcy problem.This yields an interesting advantage in complex flow problems where both equations are present: the same spaces can be used in the whole computational domain.As a consequence, the implementation of the method and the enforcement of the interface conditions are greatly simplified (see also Section 5.6).

Theoretical analysis 4.1 Interpolation estimates
In this section we prove that the following interpolation estimate holds for the enhanced space V h .Since the proof is quite involved, we divide it in three steps.
where the constant C depends only on the degree k and the shape regularity constants , c (see assumptions (A1) and (A2) of Section 3.1).
Proof.Step 1.Let w I the approximant function of v in the space W h (cf.(16) and Proposition 4.2 in [15]) then it holds that Now let v I ∈ V h be the interpolant of w I in the sense of the DoFs D V , so that Let us define ϑ := v I − w I , then for every element E ∈ Ω h the following facts hold.
• Since v I and w I are polynomials of degree k on ∂E, by definition of • Since div v I and div w I are polynomials of degree k − 1 in E, by definition of D V 4 and homogeneous boundary data (50), we get Then by definition of D V 3, we infer Now we recall that, for any v h ∈ V h , the quantity Π ∇,E k v h depends only on the values of D V (v h ), see Proposition 3.2.Therefore, using (49), we have that Π ∇,E . Thus, by ( 52) and (53), where • By definition of W E h and Collecting ( 50), ( 51), ( 54), (56) it follows that (ϑ, s, g) solves the problem Step 2. We now analyse the well-posedness of Problem (57).We consider [H 1 0 (E)] 2 and L 2 (E) endowed with the H 1 and the L 2 -norm, respectively, and G ⊕ k (E) endowed with the scaled norm Then for all ψ ∈ [H 1 0 (E)] 2 and h ∈ G ⊕ k (E) where the last inequality follows by a scaled Poincaré inequality.Therefore all the involved bilinear forms are continuous.By the theory of problems in mixed form [19], due to the coercivity of a E (•, •) the well-posedness of problem (57) will follow if we show an inf-sup condition for the form b for suitable uniform positive constants b 0 , c 0 .It is well known (see [19]) that for all q ∈ L 2 0 (E) there exists Now let T E ⊂ E be an equilateral triangle inscribed in the ball B E (cf.assumption (A1)).
Then for all polynomial p ∈ P k (E), it holds p 0,E ≤ C p 0,T E for a suitable uniform constant C. Let h ∈ G ⊕ k (E) and we define q := rot(h) and where b ∈ P 3 (T E ) denotes the standard cubic bubble in T E with unitary maximum value.Therefore, we get Since rot : ), a scaling argument for polynomials on the triangle T E yields rot(h) 0,T E ≥ h −1 E h 0,T E .Thus using (61) we find Moreover using an inverse estimate for the polynomials bq and h Therefore by ( 62) and ( 63) Recalling (59), let us set ϕ := ϕ 1 + ξ ϕ 2 (cf.( 60) and ( 64)) where ξ is a positive constant.Then, it is clear that Moreover, by (58) and since div curl = 0, we have for any positive real number ε.Finally, setting 65) and (66) we get (59).
Step 3. Since problem (57) is well-posed, the following stability estimate holds where Then, by the definition of χ (see (55)) and by the continuity of the L 2 -projection, we get where the last inequality is justified since, by definition (14), the function Π ∇,E k w I − w I has zero mean value.Noting that Π ∇,E k is a projection with respect the H 1 semi-norm and using a triangular inequality, from (48) we finally get The thesis now follows from (67) and again (48), by adding all the local contributions.For what concerns the L 2 estimate, for each polygon E ∈ Ω h , we have that ϑ = 0 on ∂E (see (50)).Hence, from (67) it holds , from which we easily infer the L 2 estimate.

Convergence analysis
First of all, let us recall a classical approximation result for P k polynomials on star-shaped domains, see for instance [20].

Lemma 4.1. Let E ∈ Ω h , and let two real numbers s, p with
with C depending only on k and the shape regularity constant in assumption (A1).Now we prove two technical lemmata.
Proof.First of all, we set ) then by definition (38) and (39) We now analyse the two terms.For the term µ(w) by simple computations, we have Now, by definition of L 2 projection Π 0,E k and by Lemma 4.1, we have Applying Hölder inequality (for sequences), we get and by Hölder inequality and Sobolev embedding H s−1 (Ω) ⊂ W s 4 (Ω) we infer By ( 73) and (74) we finally obtain For what concerns the term β(w) in (71) using Hölder inequality we have and thus, by the continuity of Π 0,E k with respect the L 4 -norm (cf.( 35)), Using again the continuity of Π 0,E k with respect the L 4 -norm, by ( 76) and (77) we infer Applying the Hölder inequality and Sobolev embeddings H 1 (Ω) ⊂ L 4 (Ω) and H s+1 (Ω) ⊂ W s 4 (Ω), we obtain For what concerns the term γ(w) in (71), using Hölder and the continuity of Π 0,E k , it holds Using again the Hölder inequality and Sobolev embedding we get By collecting (75), ( 78) and ( 80) in (71) we finally get For the second term µ 2 (w) we only sketch the proof since we use analogous arguments.First by definition, then by adding and subtracting terms, we obtain For the term δ(w) we have and applying the Hölder inequality (for sequences) we easily get The Hölder inequality and the Sobolev embedding H s+1 (Ω) ⊂ W s 4 (Ω) yield The terms ε(w) and ζ(w) can be estimated using the usual argument (Hölder inequality, continuity of Π 0,E k with respect to the L 4 -norm and Sobolev embeddings).We conclude that We infer the thesis by collecting ( 81) and ( 85) in (70).
Lemma 4.3.Let C h be the constant defined in (32).Then for all v, z, w ∈ V it holds Proof.Since c h (•; •, •) is skew-symmetric by simple computations we obtain The thesis follows by definition (32).
Furthermore, we state the following result concerning the load approximation, which can be proved using standard arguments [11].Lemma 4.4.Let f h be defined as in (40), and let us assume f ∈ H s+1 (Ω), −1 ≤ s ≤ k.Then, for all v h ∈ V h , it holds We now note that, given v ∈ Z, the inf-sup condition (43) implies (see [19]): which essentially means that Z is approximated by Z h with the same accuracy order of the whole subspace V h .In particular by Theorem 4.1, assuming v ∈ H s+1 (Ω) ∩ Z, 0 < s ≤ k, we infer inf Theorem 4.2.Under the assumptions (9) and (44), let u be the solution of Problem (12) and u h be the solution of virtual Problem (47).
where F and H are suitable functions independent of h.
Proof.Let u I be an approximant of u in the discrete kernel Z h satisfying (86), and let us define δ h := u h − u I .Now, by the stability and the consistency properties (cf.( 27) and ( 28)) of the bilinear form a h (•, •), the triangular inequality and (86) give where u π is the piecewise polynomial of degree k defined in Lemma 4.1.Now since u and u h are solutions of Problem (12) and Problem (47) respectively, from Lemma 4.4 we obtain Now we observe that The first term can be estimated by Lemma 4.2 The second term, recalling that δ h = u h − u I , is bounded by Lemma 4.3 Collecting ( 90) and ( 91) in (89), we get and then by Theorem 4.1 we infer We observe now that from ( 45) and ( 44), it holds Therefore and from (10), ( 9), ( 45) and ( 44) we finally obtain The thesis easily follows from the triangular inequality.
Remark 4.1.We observe that, due to the divergence-free property of the proposed method, the estimate on the velocity errors in Theorem 4.2 does not depend on the continuous pressure, whereas the velocity errors of classical methods have a pressure contribution.A numerical investigation of this aspect, also in relation to the presence of a higher order load approximation term in the right hand side of (87), will be shown in the next section.
Remark 4.2.From the discrete inf-sup condition (43) the pressure estimate easily follows by standard arguments.Let (u, p) ∈ V × Q be the solution of Problem ( 7) and (u h , p h ) ∈ V h × Q h be the solution of Problem (42).Then it holds: for a suitable function K(•; •, •) independent of h.
Remark 4.3.In Theorem 4.2 we have assumed u and f in H s+1 (Ω).However, it is easy to check that the same analysis can be performed if we only require: In such a case, the higher order Sobolev norms on u, f appearing in Theorem 4.2 (and in the other results of this section) are substituted with the corresponding element-wise broken Sobolev norms.

Numerical Tests
In this section we present six sets of numerical experiments to test the practical performance of the method.All the tests are performed with the second-order VEM, i.e. k = 2.We also consider suitable second order Finite Elements for comparison.In almost all cases, both options c h (•; •, •) and c h (•; •, •) (see Remark 3.3) yield very similar results; in such cases, only the first choice is reported.On the contrary, whenever the results between the two choices are significantly different, both outcomes are shown.
In Test 5.1 and Test 5.2, we consider two benchmark problems for the Stokes and Navier-Stokes equation.They share the property of having the velocity solution in the discrete space.However, classical mixed finite element methods lead to significant velocity errors, stemming from the velocity/pressure coupling in the error estimates.This effect is greatly reduced (or even neglected) by our VEM methods (cf.Theorem 4.2 and estimate (94)).In Test 5.3 we analyse the stability of the method with respect to the viscosity parameter ν.In Test 5.4 and Test 5.5 we study the convergence of the proposed method for the Navier-Stokes and Stokes equations, respectively.A comparison with the triangular P2-P1 and the quadrilateral Q2-P1 mixed finite element methods, see for example [19], is also performed.Finally in Test 5.6 we assess the proposed virtual element method for flows which are governed by the Stokes system on one part of the domain, and by the Darcy's law in the rest of the domain, the solutions in the two domains being coupled by proper interface conditions (see Remark 3.4).In order to compute the VEM errors, we consider the computable error quantities: where in the previous formula "nodes" denotes the set of internal edges nodes and internal vertexes (cf.D V 1 and D V 2).For the pressures we simply compute error(p, L 2 ) := p − p h 0 .
An example of the adopted meshes is shown in Figure 2. The distorted quadrilateral meshes are obtained starting from the uniform square meshes and displacing the internal vertexes (with a proportional "distortion amplitude" of 0.3 for Q h and 0.5 for U h ).The non-convex WEB-like meshes are composed by hexagons, generated starting from the triangular meshes {T h } h and randomly displacing the midpoint of each (non boundary) edge.For what concerns the disk Ω D we consider the sequences of polygonal meshes: • {T h } h : sequence of triangular meshes with h = 1/5, 1/10, 1/20, 1/40, • {V h } h : sequence of CVT Voronoi meshes with h = 1/5, 1/10, 1/20, 1/40.It is well known that the velocity error between the exact velocity u and the discrete velocity u h of standard mixed elements like the Q2-P1 element for the incompressible Stokes equations is pressure-dependent, i.e. has the form where C 1 , C 2 are two positive uniform constants, whereas for the virtual element scheme (see Theorem 4.2 and [15]) the error on the velocity does not depend by the pressure, i.e.
We observe that for both VEM and Q2-P1, the pressures p 1 and p 2 do not belong to the discrete pressure space.Therefore we expect that the discrete Q2-P1 velocities are polluted by the pressure approximation.Table 1 shows the results obtained respectively with VEM and Q2-P1 for the case of polynomial pressure p 1 and sequence of meshes Q h .We observe that the virtual element method yields an exact hydrostatic velocity solution, since f is a polynomial of degree two, while the Q2-P1 finite element method, in accordance with the a priori estimate (95), shows non-negligible errors in the velocity.
On the other hand we note that, due to the load approximation procedure, there is a load dependent term in the right hand side of (96).As a consequence, in the test with goniometric pressure p 2 (where the load f is not a polynomial) we expect a slight pollution of the velocity errors also for the VEM scheme, although much smaller than for the FEM case.In Figure 4 we plot the errors for the goniometric pressure p 2 and the same sequence of meshes Q h .In accordance with the a-priori estimates (95), (96) and the above observation, we obtain quadratic convergence rate for the Q2-P1 finite element method, and fourth order convergence rate for the VEM scheme for the H 1 -velocity (quadratic for the L 2 -pressure errors).In this test we consider two benchmark Navier-Stokes problems taken from [33] on the disk Ω D where we compare the results obtained with VEM discretization with those obtained with the standard P2-P1 element for the sequence of meshes T h .The solutions are chosen in such a way that the pressures balance the nonlinear convective term yielding a vanishing external load f = 0.
In the first example we take ν = 1 and the exact solution We notice that the velocity u 1 belongs to the discrete space for both VEM and P2-P1 schemes.
In Table 2  We observe that the VEM non−skew yields an exact solution u h = u 1 .Indeed, in this simple case it holds This property is not verified by the skew-symmetric trilinear form c h (•; •, •).The P2-P1 does not yield the exact velocity solution since the velocity error of the method is polluted by the approximation of the pressure.
In the second example we set ν = 1 and we consider the exact solution In Figure 5 we show the results.We are in a similar situation to the previous example (the velocity u 2 belongs to the discrete spaces, whereas the pressure p 2 does not) but with an important difference.Also in this case we observe that VEM non−skew provides a better performance than VEM skew , but now u h = u 2 .Indeed, in this case it holds and using similar steps as in the proof of Theorem 4.2, for VEM non−skew we can derive Instead, for VEM skew we can only obtain Finally, figure 6 displays the results obtained with VEM non−skew and VEM skew for the sequence of polygonal meshes V h (see Figure 3).
The aim of this test is to check the actual performance of the virtual element method for small viscosity parameters, in comparison with the standard P2-P1 mixed finite element method.Figure 7 shows that the solutions of the virtual element method are accurate even for rather small values of ν.Larger velocity errors appear only for very small viscosity parameters.The reason for this robustness is again that the "divergence free" property of VEM yields velocity errors that do not depend directly on the pressure (but only indirectly through the higher order load approximation term, see Theorem 4.2).On the contrary, for the P2-P1 element the pressure component of the error can become the dominant source of error also for the velocity field.In addition, we note that for ν = 10 −4 , 10 −5 the P2-P1 element does not even converge.In Figure 8 we show the results obtained for the sequence of triangular meshes T h , also compared with the P2-P1 element.We notice that the theoretical predictions of Sections 4 are confirmed.Moreover, we observe that the virtual element method exhibit smaller errors than the standard P2-P1 method, at least for this example and with the adopted meshes.Finally we test the virtual element method with the sequence of polygonal meshes W h , obtaining that the theoretical results are confirmed as well (note that the N dof behaves like h −2 ).We observe that for the (less deformed) quadrilateral meshes Q h , both the virtual element method and the Q2-P1 preserve the theoretical order of accuracy, but the Q2-P1 element yields better results.Instead, for the (more deformed) sequence of meshes U h , the behaviour is completely different.The virtual element approach maintains the optimal second order accuracy, whereas the Q2-P1 element clearly suffers from an evident sub-optimality of the convergence rates (the pressure does not even seem to converge).Therefore, we may conclude that the VEM seems to be more robust with respect to large distortions of the mesh.At the interface between Stokes and Darcy regions, the system (97) is coupled using the Beavers-Joseph-Saffmann condition (see [7,40] for further details).
We observe that in our test problem, we set free boundary conditions on the right boundary edge of the Darcy region.To test the performance of the virtual element method, we compute the unknown fluxes quantity (see Figure 12   In this experiment the computational domain Ω := [0, 2] 2 is partitioned using two sequences of polygonal meshes: • {Q h } h : sequence of square meshes with element edge length h = 1/4, 1/8, 1/16, 1/32 ; • {P h } h : sequence of meshes obtained by gluing a Voronoi decomposition on the domain Ω S , a triangular decomposition on Ω D1 and a square decomposition on Ω D2 , with edge length h = 1/4, 1/8, 1/16, 1/32.
In addition, we use the square mesh Q 1/64 as the basis for the reference solution.An example of the adopted meshes is shown in Figure 13.In Figure 14 we show the plot of the numerical velocity and pressure.Note that the purpose of mesh family {P h } h is to show the robustness of the proposed method when, by exploiting the flexibiliy of polygonal grids, completely independent meshes are glued together.In Table 3 we show the results obtained by using the sequences of meshes Q h and P h , compared with those obtained with the reference mesh.We observe that both sequences of meshes exhibit appropriate convergence properties, confirming that the proposed virtual element method can automatically handle non-conforming polygonal meshes and the coupling between Darcy and Stokes flow problems.[8] L. Beirão da Veiga, F. Brezzi, and L. D. Marini.Virtual Elements for linear elasticity problems.Siam.J. Numer.Anal., 51:794-812, 2013.

Figure 1 :
Figure 1: Degrees of freedom for k = 2, k = 3.We denote D V 1 with black dots, D V 2 with red squares, D V 3 with green rectangles, D V 4 with blue dots inside the element.

Figure 3 Figure 3 :Test 5 . 1 (
Figure 3 displays an example of the adopted meshes.For the generation of the Voronoi meshes

Figure 7 :
Figure 7: Test 5.3: Errors of VEM (dotted lines) and P2-P1 (solid lines), with different values of ν for the meshes T h .

Figure 8 :
Figure 8: Test 5.4: Errors with VEM and P2-P1 for the meshes T h .

Table 1 :
Test 5.1: Errors with VEM and Q2-P1 for polynomial pressure p 1 and meshes Q h .