A Hu–Washizu variational approach to self-stabilized virtual elements: 2D linear elastostatics

An original, variational formulation of the Virtual Element Method (VEM) is proposed, based on a Hu–Washizu mixed variational statement for 2D linear elastostatics. The proposed variational framework appears to be ideal for the formulation of VEs, whereby compatibility is enforced in a weak sense and the strain model can be prescribed a priori, independently of the unknown displacement model. It is shown how the ensuing freedom in the definition of the strain model can be conveniently exploited for the formulation of self-stabilized and possibly locking-free low order VEs. The superior performances of the VEs formulated within this framework has been verified by application to several numerical tests.


Introduction
The VEM, originally formulated for Poisson's problems and the Laplace operator [1], has been successively extended to linear elastostatics [2]. Among the several contributions appeared in the literature, we mention here the work in [3], with the VEM formulation in elastostatics on low order 3D polyhedra, the one in [4], where a detailed derivation of the VEM for linear 2D elastostatics is presented for the lowest order together with its extension to arbitrarily higher order on general polygons, and [5], where the finite-deformation case is considered.
While mixed stress-displacement variational formulations of the VEM for elastostatics have been proposed in the literature (see, e.g., the works [6][7][8], based on Hellinger-Reissner variational principle, and [9]), to our knowledge, more general strain-stress-displacement mixed VEM formulations have never been investigated. In this work, we propose a formulation of the VEM for linear 2D elastostatics based on the mixed Hu-Washizu variational principle [10]. Making use of Prager's notion of generalized variables [11,12], the stress model is obtained directly from the strain model and, in practice, is not anymore a primal unknown of the problem, which therefore reduces to a strain-displacement formulation, where the key ingredient turns out to be the compatibility matrix enforcing strain-displacement compatibility in weak form. Using the VEM paradigm, it is then shown how it is possible to compute the compatibility matrix and, hence, the stiffness matrix, on arbitrarily shaped polygons.
As it is well-known, in most cases the VE stiffness matrix requires a stabilization to avoid the development of zeroenergy hourglass modes. Alternative hourglass stabilization strategies have been proposed in [1] and in [13]. For the connection between the VEM and the finite element hourglass control techniques, see [14]. In this work we present an orig-inal derivation of the stabilization stiffness matrix, departing from Argyris' notion of natural strains [15,16]. The key ingredient is the hourglass matrix, which is constructed based on a decomposition of the discretized displacement modes into a rigid body and pure deformation part and into an hourglass part. The resulting stabilization hourglass matrix turns out to be identical to the one in [1]. Once the hourglass matrix has been derived, the stabilization matrix is built according to what proposed in [17].
Even though the stabilization has proved to be effective in guaranteeing the correct convergence order, its substantially empirical nature, being based on artificial stiffness coefficients, is viewed as a limitation of the method. For this reason, there is a significant interest in investigating the possibility to formulate self-stabilized VEs, i.e., not requiring any artificial stabilization. In view of the weak enforcement of compatibility, in the proposed Hu-Washizu formulation the strain model is independent of the displacement one. It is then rather natural to exploit this additional freedom for the formulation of self-stabilized VEs. In the classical VEM, if the displacement model contains polynomials of order k, the strain model is assumed to be polynomial of order k − 1. In this work, a strain model of order p > k − 1 is assumed, leading to a stiffness matrix of the correct rank. However, we remark that for arbitrary polygons and polynomial order, a satisfying theoretical investigation on how the strain model should be selected in order to have the correct rank is still missing. While in the standard VEM, only nodal Degrees Of Freedom (DOFs) are required in the case k = 1, the price to pay for the self-stabilization is that additional, moment-type, DOFs have to be introduced in the element formulation. Even though these new DOFs can be easily condensed out, since they are internal to the element and not shared with neighboring elements, their presence implies a small additional computational effort. To avoid this, an alternative approach, based on a projection of the unknown virtual displacement field, based on what very recently proposed in [18,19], has also been considered, leading to a self-stabilized element, not requiring additional DOFs. It should be noted that while this work was in progress, an approach substantially identical to the one considered here, though not based on a Hu-Washizu formulation of the VEM, has been published in [20], leading to the same stabilized stiffness matrix in the case k = 1.
Another interesting aspect to be investigated is the behavior of the developed VEs in the incompressibility limit. Wriggers et al. [13] presented a VEM for large strain elasticity, with excellent locking-free behavior in the nearly incompressible limit and fully locking-free behavior in the incompressible case, in this latter condition by means of a penalization of the incompressibility constraint using pressure as a Lagrange multiplier. For the case of isotropic elasticity, in Park et al. [21] proposed a decomposition of the elastic tensor in deviatoric and volumetric parts. The sta-bilization stiffness matrix is then constructed as in [17], but using only the deviatoric part of the elastic tensor. The resulting VEM model is shown to provide accurate results also in the nearly incompressible case. In this work, we show through numerical tests that the new VEs, enriched with additional moment DOFs and self-stabilized, seem locking-free also in the nearly incompressible case and for highly distorted and/or non-convex element shapes.
The paper is organized as follows. The Hu-Washizu variational formulation of the 2D linear elastic continuum problem is recalled first. Then, the corresponding three-field finite element discretization, together with the generalized variable assumption, is introduced. It is then shown how the VEM paradigm can be used, starting from the defined mixed framework, to construct the element matrices for elements of arbitrary polygonal shapes. The mixed variational nature of the proposed formulation naturally leads to the derivations of Sect. 3, where two types of low-order self-stabilized 2D VEs are formulated. Numerical tests for the compressible and nearly incompressible case show the excellent performances of the new self-stabilized elements.
Throughout this work, Voigt notation is adopted, so that stress components are gathered in the stress vector σ and strain components in the strain vector ε. Furthermore, the material elastic tensor is replaced by the material matrix of elastic moduli D.

Hu-Washizu variational formulation of the continuum problem
Let us focus on the 2D linear elastostatic continuum problem, under the usual assumptions of small displacements and strains. Inelastic and thermal effects are not considered for simplicity, though they could be easily incorporated into the theory.
The solid body is represented by a domain ⊂ R 2 , whose boundary ∂ is composed of a constrained part ∂ u and a free part ∂ p , with ∂ u ∩ ∂ p = ∅ and ∂ u ∪ ∂ p = ∂ . On the constrained part, imposed displacementsū are assigned; on the free part, surface tractions p are applied. The body is also subjected to body forces b. The two inplane displacement components are collected into the vector u. All the aforementioned data and unknowns depend on the position vector x with respect to a Cartesian reference system.
The starting point is the definition of the three-field Hu-Washizu functional, assuming as independent variables displacements u, strains ε and stresses σ : with u =ū on ∂ u . In (1), D is the matrix of elastic constants and S denotes the compatibility differential operator, defined as: where ∂ (·) denotes the partial derivative with respect to (·). Its transpose S T is the equilibrium differential operator. Hu-Washizu variational theorem states that among all solutions, the one satisfying equilibrium, compatibility and constitutive law, makes the three-field functional stationary: Integrating by parts the integral involving the term Sδu and thanks to the arbitrariness of the variations δu, δε, δσ , the weak form of the governing equations is obtained: ∀δσ compatibility (6) where N is the matrix of director cosines of the outward normal to the boundary: and n x and n y are the two components of the outward normal unit vector n.
The set {u, ε, σ } that makes stationary the Hu-Washizu functional is then the solution of the problem, corresponding to a saddle point for the mixed functional.

Mixed finite elements based on the Hu-Washizu principle
The problem domain is subdivided into n e finite elements, each occupying a domain e . Let ∂ p e be a possible part of e boundary coinciding with a portion of the body boundary ∂ p subjected to surface tractions. Let ξ be a vector of nondimensional local coordinates in 2D: where x G and y G are the coordinates of the element centroid and h e is the maximum diameter of the element, i.e. a measure of the element size.

Remark 1
It should be noted that, unlike in standard isoparametric finite elements, no geometry mapping is required in the definition of virtual finite elements; the intrinsic variables ξ and η defined here have therefore substantially different meaning than in classical isoparametric formulations. This aspect allows for strongly distorted elements which are forbidden in the isoparametric setting.
In the spirit of the Hu-Washizu formulation, an independent modelling of local displacements, strains and stresses is introduced: where N u , N ε , N σ are matrices of shape functions defined in e whose dimensions are respectively 2 ×n u , 3×n ε , 3×n σ , n u , n ε and n σ being the number of parameters used to define the corresponding discretized fields.û,ε andσ are vectors of parameters, in general not coinciding with nodal values and without a physical meaning. While u h (ξ ) is required to be C 0 continuous across elements, the interpolation functions contained in N ε and N σ are continuous inside each element, but may not be so across element boundaries. To simplify the notation, the subscript e, denoting the considered element e , will be omitted unless strictly necessary.
In the proposed formulation, strain and stress parametersε andσ are required to correctly represent the element energy, in the sense that: Parametersσ andε satisfying this condition are said to be generalized variables in the sense of Prager. 1 Equation (12) implies that n σ = n ε and that: where I is the n ε ×n ε identity matrix. Consequently, possible choices for the stress shape functions N σ are: where the square and invertible n ε × n ε matrices E (hereafter referred to as elasticity matrix) and G are defined as: The second choice (15) will be used in the remainder of this paper.

Remark 2
The choice (15) for N σ has been proposed by Corradi in [12] and used, e.g., in [22,23] in the framework of variational elastoplasticity in generalized variables. The choice (14) has been proposed in [23] and recently used in [20] for an enhanced VEM formulation.
The final expression of the discretized functional associated to the generic element e is: where (13) has been exploited and the following quantities have been introduced: • element compatibility matrix, enforcing compatibility in weak form 1 The concept of generalized variables was introduced by Prager with reference to frame structures (see, e.g. [11,22]).
• element equivalent nodal forces vector (20) where F b and F p are the contributions to the equivalent nodal forces vector coming from body forces b in e and surface tractions p on ∂ p e , respectively.

Remark 3
If choice (15) for N σ is adopted, the construction of the compatibility matrix C in (18) can be seen as resulting from an L 2 projection of the symmetric part of the displacement gradient onto the discretized strain space.
The governing equations in discretized form are obtained by enforcing the stationarity of the discretized functional (17) with respect toû,ε andσ : ∂ε h e = 0 ⇒σ = Eε constitutive law (22) ∂σ h e = 0 ⇒ε = Cû compatibility (23) Replacing (23) in (22) and (22) in (21), one obtains: where: is the element stiffness matrix, symmetric and positive semidefinite, consistent (the superscript c stands for 'consistent') with the displacement and strain models. If n u − n ε ≤ 3 and C has n u − 3 independent rows, K c has the correct degree of singularity, equal to 3, i.e., equal to the number of rigid body modes in 2D. In contrast, if n u − n ε > 3, K c has a surplus of rank deficiency, equal to n u − n ε − 3, and zeroenergy (hourglass) modes can arise. 2 Hourglass modes are spurious deformation modes without physical meaning and are pathological. Hourglass stabilization is a typical issue of mixed finite element formulations and it has the role of reestablishing the correct degree of singularity of the stiffness matrix without upsetting accuracy.

Remark 4
Once the displacement DOFsû have been computed, the proposed mixed variational framework offers a straightforward strategy for stress and strain recovery. Using the strain and stress models N ε and N σ , one immediately has:

Hourglass stabilization
The basic idea to stabilize the approximate solution is that of adding a fictitious stiffness to the element hourglass modes. In this way, the expression of the discretized element mixed functional becomes: The vectorû H contains combinations of displacement parameters defining hourglass modes and is a matrix of fictitious stiffnesses. A possible choice for the definition of (see, e.g., [3]) is (scalar-based stabilization): where I is the n u × n u identity matrix, tr(·) denotes the trace operator and the coefficient 1/2 has been proposed in [4]. Another choice, proposed in [17], is (diagonal matrix-based stabilization): The hourglass modesû H in (28) can be defined in terms of the displacements DOFs as: An original construction of the hourglass matrix H, based on the natural approach of Argyris [12,15,16], is detailed in "Appendix A". The resulting stabilization stiffness matrix is identical to the one originally proposed in [1,17]. It is also worth mentioning the alternative approach to VEM stabilization originally proposed in [13]. For a discussion on the connections between the stabilization of hourglass modes in the VEM and in the finite element method, see [14]. Defining the local stabilizing stiffness matrix K s : and using the definition (25) of the consistent stiffness matrix and (28) of the stabilized functional, the stabilized element stiffness matrix K can be expressed as the sum of two terms: This stiffness matrix is symmetric and has now the correct rank deficiency, equal to 3.

Virtual element formulation
The general formulation presented in Sects. 2.2 and 2.3 can be used to formulate different mixed finite elements. In particular, this framework will be exploited here to construct mixed finite elements based on the Virtual Element Method (VEM). Let e denote the domain of a finite element extracted from a subdivision of the body domain in non-overlapping polygons having straight edges and arbitrary shape. This element is characterized by: • an arbitrary polygonal shape • an arbitrary number N V of vertices and straight edges.
Element local coordinates are defined by the simple linear transformation (8), without a geometry mapping from a parent element.
In the VEM, displacement shape functions N u are not explicitly known inside the element and are therefore said to be virtual. These functions are known only on the polygonal boundary of the element ∂ e , where they are described by a polynomial of degree k (order of accuracy of the method). The displacement field inside a virtual element is assumed to contain a complete polynomial of degree k plus other functions that, however, are never required to be explicitly defined in the element interior and are not used in the computation.
For k > 1 and for the self-stabilized VEs with k = 1 that will be discussed in Sect. 3.1, the VEM degrees of freedom (DOFs) are not just nodal displacements. The DOFs in a virtual element are: • nodal displacements at the element vertices; • nodal displacements at the nodes inside the element edges (i.e., at positions different from the edge extrema); • moments (scaled with respect to the element area e , to be defined later) of order up to k − 2 of the unknown approximate displacement field u = N uû . These are present only if k > 1 and are often referred to as 'internal DOFs'.
Since a polynomial of degree k in 2D requires for its definition n k parameters with: the displacement field turns out to be defined by means of n u DOFs, with: is the number of parameters necessary to describe a polynomial of degree k − 2 inside the element.
With the previous assumptions on the element displacement DOFs, the matrix N u defining the virtual displacement model is of the type: The first N V 2 × 2 blocks are associated to vertex DOFs; the successive N V (k − 1) 2 × 2 blocks are associated to edge DOFs; the last k(k −1) 2×2 blocks are associated to internal moment DOFs. For the subsequent developments, it is convenient to collect the monomials ξ a η b , with a + b ≤ k, in the following row vector: With this notation, the entries of q k are denoted as q j . For instance, for j = 4 one has q 4 = ξ 2 , etc. The last n k−2 = k(k−1) 2 shape functions in (37) (associated to internal moment DOFs) are assumed to have value 0 on the boundary and have scaled moment equal to 1 in their DOF and equal to 0 in correspondence of all the other DOFs. These last conditions can be expressed as (for j = 1, . . . , where δ i j is the Kronecker's delta and q j are the monomials in the vector q k−2 (38).
The virtual shape functions contained in the first kN V 2 × 2 blocks, (i.e., relative to a boundary node) are unknown in the element interior but are assumed to be a Lagrangian polynomial of order k on the element edges, taking value 1 in their node of definition and 0 in all the other nodes. Furthermore, their moments up to order k − 2 are assumed to be equal to 0. The strain field model N ε (ξ ) in a Hu-Washizu virtual element is defined a priori. If the displacement model is assumed to contain at least a complete polynomial of degree k, in the standard VEM the strain model is assumed to be a complete polynomial of degree k − 1. Thus, the number of parameters required to describe the strain field in 2D is: The matrix N ε (ξ ) can be expressed as: The key operator in the present VEM formulation is the compatibility matrix C defined in (18), projecting the symmetric part of the displacement gradient field onto the space P k−1 of polynomials of degree up to k − 1 of the approximate strain field. The computation of C requires the computation of the symmetric and invertible matrix G, defined in (16), and of A, defined in (19). The matrix G is directly computable once the degree of accuracy k of the method is defined, by computing the integrals by means of a subtriangulation technique. For more general strategies of polynomial integration on polygons, see for instance [24]. According to this procedure, the element is subdivided into N V subtriangles obtained by connecting the element centroid to the element vertices. 3 Then, a standard Gaussian quadrature rule for each triangle is applied, in such a way that polynomials of degree 2k − 2 are exactly integrated. The computation of the matrix A is less straightforward, since its expression contains the matrix of displacement shape functions N u , unknown in the element interior. To overcome the problem, one can proceed integrating by parts: Matrix A 1 results from the boundary integral of known quantities, since N u i (ξ ) is assumed to be a polynomial of degree k on ∂ e , if i denotes a boundary node, and to vanish on ∂ e , if i denotes a moment DOF. This integration can be performed exactly for the polynomials of order 2k −1 resulting from the product of polynomials of degree k − 1, related to N ε , and polynomials of degree k, related to N u , by the Gauss-Lobatto quadrature rule, using the k + 1 edge nodes as integration points. Since the displacement shape functions assume just values 1 and 0 in correspondence of the boundary nodes, the integration turns out to be remarkably simple.
For what concerns A 2 , the first part of the integrand is: that can be also expressed in the scaled monomial basis q j (ξ ) as: where each of the n k−2 matrices M j contains some zero entries and other non-zero entries that are constant values and proportional to 1/h e . Replacing (47) in the expression of A 2 , one obtains: The integrals in (48) are the moments of order up to k − 2 of the displacement shape functions N u , easily computable exploiting (39)-(42).

Remark 5
From the expression (45) of A 2 , it is clear that the required number of internal DOFs is strictly related to the order of the polynomials assumed for the modelling of the strain field in N ε . In standard VEM, since the strain field is modelled by a complete polynomial of degree k − 1, the internal moments have to be assumed up to the order k − 2. This observation will be relevant in Sect. 3, when dealing with self-stabilized virtual elements, in which the enrichment of the strain field model will lead to a number of internal displacement DOFs higher than in standard VEM.
Once the compatibility matrix C and the elastic matrix E have been computed (the latter by means of the subtriangulation technique already described for the matrix G), one can immediately compute the consistent stiffness matrix K c from (25). As for the stabilization stiffness matrix K s in (32), 4 whether one chooses for the scalar-based stabilization or the diagonal matrix-based one, its computation is basically reduced to the computation of the hourglass matrix H. The computation of the element vector of equivalent nodal forces F can be carried out following the standard approach used in the VEM (see e.g. [4]) and it is briefly summarized in "Appendix B".

Self-stabilized virtual elements k = 1
One of the main limitations of the standard VEM is the need for a stabilization of the stiffness matrix. Focusing on virtual elements with k = 1, with 4 and 5 nodes (k = 1 triangles do not require to be stabilized), two methodologies are therefore proposed for the formulation of self-stabilized elements. The main difference between the two is related to the introduction or not of additional internal degrees of freedom. The first category of self-stabilized elements introduces additional internal moments as displacement DOFs. It is numerically shown how elements of this kind: • do not require the stabilization of the stiffness matrix • exhibit a superior accuracy with respect to standard VEM • are locking-free, in the sense that the approximated displacement field does not suffer volumetric locking in the presence of nearly incompressible materials.
The second category is based on vertex displacement DOFs only as the standard VEM with k = 1. It is numerically shown how these elements: • do not require the stabilization of the stiffness matrix • exhibit an accuracy similar to the standard VEM • are not volumetric locking-free.
In both cases, the key idea is to enrich the strain field with respect to standard VEM, exploiting the freedom offered by the mixed Hu-Washizu approach, where the strain field may be defined independent of the displacement field and, therefore, has not to be a polynomial of order k − 1 as in standard VEM. In the following, the maximum degree of interpolation of the strain field is indicated by p, while k still denotes the polynomial degree of interpolation of the displacements along the boundary of the virtual element.

Remark 6
The first approach proposed here is substantially identical to the 'Uncoupled Polynomial Representation' strategy just proposed in [20], while this work was in progress, the only difference resting in the choice of the stress model N σ . The option in (14) has been adopted in [20], leading to an energy projection of the displacement gradient, while the choice for N σ in (15) has been adopted in the present paper. The second approach proposed is basically an extension to plane elasticity of the strategy presented in [18] for the Poisson equation, similarly to what very recently done in [19].

Remark 7
The first self-stabilization procedure proposed here could be straightforwardly extended to higher order displacement models, though at the cost of significant additional computational burden. For an 8-node quadrilateral, with quadratic interpolation on the edges (k = 2), e.g., a cubic model ( p = 3) would be required for the strain field, with 30 strain parameters and 28 displacement DOFs (16 nodal displacements plus 12 internal moments) [20]. For this reason, only k = 1 elements will be considered here.
The second procedure could also be used in the presence of higher number of edges and/or higher order displacement models, always keeping the same number of DOFs as the corresponding standard VEM of the same order. However, as already mentioned in the Introduction, we remark that a clear and satisfying analysis on how to design stabilization free VEs in the general case (i.e. arbitrary polygons and polynomial order) is still missing.

Self-stabilized elements with additional internal degrees of freedom
First, let us focus on k = 1, 4-node elements with additional internal moment DOFs. The acronym VEM4SS is used to denote 4-node self-stabilized virtual elements with linear displacement interpolation along the edges. In standard VEM, a constant strain field (of order k − 1 = 0) is assumed and no moment DOFs are necessary, since matrix A 2 in (45) vanishes. As discussed in Sect. 2.2, a necessary condition to have a self-stabilized element is that n u − n ε ≤ 3. In standard VEM, for this type of elements one has n u = 8 and n ε = 3, hence the consistent stiffness matrix K c contains two hourglass modes and requires a stabilization. The need for a stabilization can be eliminated by enriching the strain field with linear terms (i.e., p = 1) in one of different ways, at the cost of two additional moment DOFs for the displacement model. Since the strain model now contains polynomial terms of degree p = 1, unlike in the standard VEM with k = 1, from (45) one has that the first moments of the displacement shape functions have to be considered for the computation of matrix A 2 , implying that two moment DOFs are required. Considering that N V = 4 and that there are 2 displacement DOFs per vertex and 2 internal moment DOFs, the final number of displacement DOFs is n u = 2 × 4 + 2 = 10. Hence, N u is a 2 × 10 matrix.
As a first attempt, the following 3 × 7 strain model is assumed: where the first three columns define a constant strain state (as for standard VEM with k = 1), the fourth and the fifth columns correspond to the two hourglass modes of the 4node element and the last two are necessary to define a complete first order polynomial for each strain component and to make unnecessary the stabilization of the consistent stiffness matrix K c . Therefore, in this case p = 1 and n ε = 7. This type of element will be denoted as VEM4SS7-10DOFs, emphasizing the fact that 7 strain parameters and 10 displacement DOFs are considered. Alternatively, also the following strain model can be considered: In this case p = 1 and n ε = 9, and this type of element is denoted as VEM4SS9-10DOFs, since 9 strain parameters and 10 displacement DOFs are still considered. For both the proposed elements, n u − n ε ≤ 3 and, if the rows of the compatibililty matrix C are independent, the element consistent stiffness matrix K c has rank deficiency 3, i.e., the element is self-stabilized. Indeed, denoting byû R i the generic vector of displacement parameters corresponding to one of the three rigid body modes, since SN uû R i = 0 by construction, one immediately has that: even in the case n u − n ε < 3. The whole procedure is analogous to the one in the standard VEM. However, it is worth observing that in this case, being k = 1 and p = 1, the integrands in the matrix A 1 are polynomials of degree k + p = 2 along each edge. In this case, the number n of integration points necessary for the exact integration of A 1 can be obtained from the condition: and if n is not an integer, the the smallest greater integer number has to be considered. In this case, n = 3, i.e., the two vertex nodes and an additional point at the middle of the edge are required. In conclusion, the construction of these elements require a boundary integration of higher order than in the standard VEM, while matrix A 2 is exactly computable thanks to the introduction of the two internal moment DOFs. The final difference with respect to the standard VEM is related to the construction of the equivalent nodal force vector F b due to the body forces b. The idea is that of projecting the body force vector b onto the space P p−1 of polynomials of degree up to p − 1. In this case, since p = 1, the vector b is projected onto the space P 0 of constant polynomials leading to a simple computation in terms of the two internal moment DOFs. It is worth noting that the computation of F b is basically the same as for standard VEM with k = 2 (see, e.g., [4]). The only difference is in the dimensions of the matrix N u . Indeed, for a standard VEM quad element with k = 2, N u is a 2 × 18 matrix, while in this case N u is a 2 × 10 matrix. In both cases, the only non-zero components of F b are those related to the two internal DOFs.
Exploiting the same idea illustrated for the 4-node element, also a 5-node element has been developed, with the acronym VEM5SS indicating that it is a 5-node selfstabilized virtual element. Also in this case, in the presence of a linear strain model p = 1, two internal DOFs are required for the computation of matrix A 2 . Since on the element boundary one has k = 1 and N V = 5, the number of displacement DOFs is n u = 2 × 5 + 2 = 12. To satisfy the condition n u −n ε ≤ 3, at least 9 strain parameters are needed and the model N ε in (50) with n ε = 9 has to be adopted.
Since n u − n ε = 3, the element consistent stiffness matrix K c has rank deficiency 3 and the element is self-stabilized if the rows of C are independent. Also in this case, the whole procedure is analogous to the one for the standard VEM. The same observations made for the VEM4SS-10DOFs elements are valid also in this case and will not be repeated here.

Remark 8
The self-stabilized VEs proposed above may seem somehow similar to the enhanced strain finite elements of Simo and Rifai [25]. There are however substantial differences. Enhanced strain finite elements are based on a nonlinear geometry mapping, as all isoparametric finite elements, and hence, they are unavoidably subject to distortion sensitivity when the geometry transformation Jacobian becomes singular. In contrast, in the proposed approach no geometry mapping between the master and the current element is required, leading to a self-stabilized VE significantly robust with respect to mesh distortion. Furthermore, in the enhanced strain approach, strains are composed of two terms: a compatible strain plus an enhanced, incompatible term, orthogonal to the stress field. Since the symmetric gradient of displacement is virtual in the VEM, no 'compatible' part of the strain can be explicitly defined in the assumed polynomial strain model. Finally, the proposed VE formulation is valid for general polygonal elements, not only for quadrilaterals, even though a systematic application to elements with more than five edges may be not computationally convenient.

Remark 9
The proposed self-stabilized VEs have some similarities also with the Pian-Sumihara (P-S) element [26]. The P-S element is based on a Hellinger-Reissner formulation, with an elementwise polynomial bilinear interpolation for the displacements and an elementwise 5-parameter modeling of the stress field. The stress can be eliminated at the element level, leading to a symmetric positive system in displacements only. The main difference with the proposed VEM model is the definition of the stress field. In the proposed Hu-Washizu formulation, the stress model is a consequence of the strain model, which requires at least 7 parameters to produce a stable element, while the standard P-S element uses 5 parameters for the stress definition. Moreover, while the P-S element is known to be less sensitive to mesh distortion with respect to standard quadrilateral elements, it is still based on a nonlinear geometry mapping. Finally, as already commented in the previous Remark, unlike the VEs, also the P-S formulation applies to quadrilaterals only.

Self-stabilized elements without additional internal degrees of freedom
These elements differ from those described in the previous Sect. 3.1 only for the computation of matrix A 2 in (45) and of the local equivalent nodal forces vector due to body forces F b . Matrix A 2 in (45) contains the integral of the unknown displacement shape functions N u and cannot be computed as it is. According to the proposed strategy, this term is computed by means of a suitable projection of the gradient of N u onto the gradient of known polynomial functions N 1 of order 1 [18,19]. Let N 1 (ξ ) be the 2 × 6 matrix of monomials up to the first order: Let u 1 (ξ ) be the approximate displacement field defined as where ∇ s denotes the symmetric gradient operator whose matrix representation is S, onto u 1 is defined as: from which one can write: The matrix in square brackets at the left hand side is obviously singular in correspondence of the three rigid body modes. Therefore, in order to solve system (55) forŝ, we need to add three other independent conditions, capable to fix the rigid body modes. A way to perform this step is the following. An identical term is added at both equation sides: where the matrix R 1 is such thatŝ R = R 1ŝ , whereŝ R are combinations of parameters defining rigid body modes. Based on this definition, wheneverŝ represents a pure deformation mode, only the first addend of matrix G ∇ s in (56) comes into play. In contrast, ifŝ represents a rigid body mode, only the second one intervenes and, in this way, the symmetric gradient projection is not modified. The construction of the matrix R 1 can be performed following the same procedure used for the computation of the hourglass matrix H illustrated in "Appendix A". One can set: whereŝ D are combinations of parameters defining pure deformation modes. It is worth noting that N 1 does not contain hourglass modes. The displacement parametersŝ can be expressed in terms of natural parametersp R 1 andp D 1 as (see "Appendix A", Eq. (A1)): The 6 × 3 matrix T R 1 of rigid body modes is given by: where the first two columns define the rigid body translations along ξ and η, respectively, and the third one defines the inplane rigid body rotation. By virtue of the orthogonality of rigid and deformation modes, one has: from which the 3 × 1 vector of parametersp R 1 can be derived: The 6 × 1 vectorŝ R can be eventually extracted from the total vectorŝ as follows: Hence, the 6 × 6 matrix R 1 required in (56) is fully computable by means of the matrix T R 1 , whose expression is given in (60). Note that wheneverŝ represents a pure deformation mode, one hasŝ R = R 1ŝ = 0.
The matrix G ∇ s in (56) is now invertible and, hence, the vector of parametersŝ can be finally computed as: where the matrix ∇ s 1 = (G ∇ s ) −1 A ∇ s defines the symmetric gradient projection operator.
The matrix G ∇ s can be easily computed, since it contains integrals over e of constant quantities, as well as integrals over ∂ e of known functions. On the other hand, the first part of the matrix A ∇ s is not immediately computable, since it contains the integral over e of the displacement virtual shape functions N u . However, one can integrate this term by parts: so that: where both terms are straightforwardly computable since N u is known on the element boundary. Adopting the usual Gauss-Lobatto rule, the first integral is exactly computable using 2 integration points over each edge (that are also boundary nodes of the element), while the second one requires an additional integration point at the middle of the edge, for a total of 3 integration points over each edge.
Once ∇ s 1 in (64) is computed, the term A 2 in (45) is rewritten as: and is immediately computable. It is notable that the computation of the last integral in (67) turns out to be very simple, due to the fact that: from the definition of centroidal coordinate system. It finally results: Once the matrix A is computed, the compatibility matrix C is immediately determined using Eq. (18) and the local consistent stiffness matrix can be also computed, following Eq. (25).
Having the same displacement DOFs as the standard VEM element with k = 1, the local equivalent nodal forces vector due to body forces F b can be computed as in the standard VEM.
The approach described above has been implemented for the 4-node element with both 7 and 9 strain parameters, using (49) and (50), and for the 5-node element, using (50). Consistently with the nomenclature used in 3.1, these three elements are respectively indicated by the acronyms VEM4SS7-8DOFs, VEM4SS9-8DOFs and VEM5SS-10DOFs. All these elements are self-stabilized since the condition n u − n ε ≤ 3 is satisfied.

Numerical tests
The mixed Hu-Washizu procedure for self-stabilized virtual elements presented in the previous Sections has been implemented into a MATLAB code for k = 1. For all the numerical tests, the diagonal matrix-based stabilization technique has been adopted in the case of the standard VEM. For all the three considered tests, the comparison between the standard VEM k = 1 and the self-stabilized elements k = 1, p = 1 is proposed, evidencing the better performances of the latters.
Units for the quantities in the examples will not be specified, though they have been taken in a consistent way (e.g. N/mm 3 for body forces; N/mm 2 for surface tractions, stresses, Young's modulus and Lamé constants; mm for lengths).

Convergence test with known analytical solution
The first application of the VEM is related to a classical 2D plane strain convergence test with known analytical solution. Specifically, the problem domain, depicted in Fig. 1, is a unit square = [0, 1] 2 with constrained displacements all over its boundary ∂ u ≡ ∂ , i.e., ∂ p = ∅.
The data of the problem are: • Lamé constants λ = 1 and μ = 1 (corresponding to E = 2.5 and ν = 0.25) • body forces in b x = −π 2 [−(λ + 3μ) sin(π x) sin(π y) +(λ + μ) cos(π x) cos(π y)] b y = −π 2 [−(λ + 3μ) sin(π x) sin(π y) +(λ + μ) cos(π x) cos(π y)] (70) • kinematic boundary conditions on ∂ u ≡ ∂ The analytical solution of the problem in terms of displacements in is given by: Different meshes have been tested for the assessment of the VEM convergence, each of them with an increasing number of elements: a square mesh, a mesh with convex distorted quadrilateral elements, a mesh with convex and non-convex quadrilateral elements and two pentagonal meshes with increasing level of distortion. All these meshes are depicted in Fig. 2. Convergence upon mesh refinement has been assessed in terms of the L 2 -norm of the strain error, defined as: where ε and ε h denote respectively the exact and the approximated strain field over the generic element of the virtual element mesh. The integrals in (73) are computed numerically by means of the usual subtriangulation technique, evaluating the exact and the approximate strains at the quadrature points. Convergence upon mesh refinement has been evaluated also in terms of a displacement error e u , defined as: where U(x v ) is the global vector containing the exact displacement solution evaluated at the mesh nodes v, whereas U h (x v ) is the VEM solution in terms of global nodal displacement DOFs. The results of VEM convergence analyses show that in all cases, irrespective of the level of element distortion, the slope of the error e L 2 tends to the order of approximation k of the method as the mean element size h decreases, when plotted in log-log scale as a function of h. Concerning the displacement error e u , the slope tends to k + 1 as h decreases, when plotted with the same scales used for e L 2 . The standard VEM with k = 1 has been compared to the self-stabilized elements presented in Sect. 3. Figure 3 shows the convergence curves of the strain error e L 2 for the different considered meshes. VEM4 and VEM5 refer to the standard quadrilateral and pentagonal VEM with k = 1, respectively. As can be seen, all the self-stabilized elements exhibit the expected convergence rate k = 1 of the standard VEM. The self-stabilized elements with additional internal DOFs always exhibit a higher accuracy than the standard VEM. Furthermore, in the case of the very simple square geometry, the self-stabilized quadrilateral element with 9 strain parameters and 10 displacement DOFs (VEM4SS9-10DOFs) shows superconvergent behaviour, with a doubled slope with respect to the theoretical value 1. As no clear theoretical explanation for this superconvergent behavior is available, one can think that this is the consequence of the combination of several ingredients: the symmetry of the analytical solution, the symmetry of the square mesh, 'aligned' with the analytical solution, etc.. The self-stabilized elements without additional moment DOFs exhibit an accuracy that is very similar to the standard VEM in all cases, but without requiring the stabilization.
The comparison of standard VEM and self-stabilized VEM in terms of displacement error e u is shown in Fig. 4. As can be seen, the expected convergence rate k + 1 = 2 is exhibited by all the considered elements.

Divergence-free convergence test
The second application of the VEM is related to a 2D plane strain convergence test with known displacementdivergence-free analytical solution. The problem domain is the same as for the first application, already depicted in Fig. 1, namely a unit square = [0, 1] 2 with constrained displacement all over its boundary ∂ u ≡ ∂ .
The data of the problem are: • Lamé constants λ = 9999 and μ = 1 (corresponding to E = 2.9999 and ν = 0.49995) • body forces in ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ b x = −4μπ 3 sin(π y) cos(π y)[cos 2 (π x) −3 sin 2 (π x)] b y = 4μπ 3 sin(π x) cos(π x)[cos 2 (π y) −3 sin 2 (π y)] (75) • kinematic boundary conditions on ∂ u ≡ ∂ The analytical solution of the problem in terms of displacements in is given by: u x = 2π sin 2 (π x) sin(π y) cos(π y) u y = −2π sin 2 (π y) sin(π x) cos(π x) The tested quadrilateral and pentagonal meshes are the ones already depicted in Fig. 2. Figure 5 shows the strain error curves related to the different meshes. As can be appreciated, the standard VEM exhibits severe volumetric locking behaviour for all the meshes. This may be not so for other VEM implementations, and in particular for k > 1 (see, e.g., [2]). Also the quadrilateral self-stabilized VEM without additional internal DOFs shows locking behaviour, while its pentagonal version exhibits the correct convergene rate. On the other hand, the self-stabilized VEM with additional internal DOFs is always volumetric-locking-free, able to keep the right convergence rate and more accurate than the other two approaches, as a result of the richer displacement field.

Cook's beam problem
The last numerical application concerns the classical Cook's beam problem. The geometry of the problem is shown in Fig. 6. It consists of a tapered cantilever beam, having the left end restrained in both directions and the right edge subjected to a uniform shear action along y. Also in this case plane strain conditions are assumed.
Since the closed-form solution of this problem is not available, convergence has been assessed by comparison of the vertical displacement u A y of the point A in Fig. 6 with its reference value taken from the literature. Both the compressible and the nearly incompressible cases have been analyzed.
The data of the problem are: • Young's modulus E = 70 • Poisson's ratio ν = 0.33 (compressible case), ν = 0.49995 (nearly incompressible case) • zero body forces in  • kinematic boundary conditions on ∂ u The considered quadrilateral and pentagonal meshes for the standard and the self-stabilized VEM with k = 1 are shown in Fig. 7. First, the compressible case (ν = 0.33) has been considered. In this case, the reference solution is u A y ≈ 0.0323. The corresponding results are shown in Fig. 8 in terms of the output parameter u A y as function of the mean element size h. As can be seen, the self-stabilized VEM with additional internal DOFs converges to the right value of the vertical displacement at A, always with enhanced accuracy with respect to the corresponding standard VEM. On the other hand, the selfstabilized VEM without additional moment DOFs always shows a slightly worse accuracy with respect to the standard VEM.
Even more interesting are the results in the nearly incompressible limit (ν = 0.49995), depicted in Fig. 9. In this case, the standard VEM with k = 1 always exhibits severe volumetric locking behaviour. Also the self-stabilized VEM without additional moment DOFs shows locking behaviour, more severe in the case of quadrilateral elements. Consistently with the results of 4.2, the self-stabilized VEM with additional moment DOFs is locking-free for all the different meshes and it rapidly converges to the right value of the tip displacement u A y ≈ 0.0277. Figure 10 shows the mean stress contour maps for the unstructured quad mesh with 498 elements and the different tested 4-node elements. The mean stress σ m is computed as (σ x + σ y + σ z )/3, where σ z = ν(σ x + σ y ) due to the plane strain hypothesis. The same quantity is reported in Fig. 11 for the distorted pentagonal mesh with 256 elements and the different tested 5-node elements. From these figures, one can note how only the self-stabilized VEM with additional moment DOFs seems to be able to completely remove the locking artifacts in the nearly incompressible situation.

Conclusions
A mixed Hu-Washizu variational formulation of the VEM has been presented. It allows in a straightforward way to cast the VE approach within the framework of mixed methods with a weak enforcement of compatibility, highlighting the role of the VEM for the computation of the compatibility matrix.
One of the main drawbacks of the VEM is that in most cases the VEs require a stabilization. While on one hand an original presentation of the stabilization technique, based on the natural approach of Argyris [15] and Corradi [12] has been proposed, on the other hand it has been shown how the formulation of the VEM as a mixed method quite naturally leads to the derivation of self-stabilized virtual elements, i.e., not requiring any stabilization. The basic idea is to increase the order of polynomial representation of the strain model, similarly to what has been proposed very recently in [20]. The resulting virtual element requires additional moment DOFs. To avoid this, an alternative technique has been proposed for the computation of the compatibility matrix, not requiring additional moment DOFs. The technique is based on a projection of the symmetric gradient of the displacement. A substantially identical technique has been very recently proposed in [18,19].
Quadrilateral and pentagonal self-stabilized, k = 1, p = 1, 2D VEs have been implemented in an in-house code and applied to a number of benchmark problems. The expected order of convergence has been obtained in all cases, with the elements stabilized by the addition of moment DOFs exhibiting a superior accuracy. In all cases the substantial distortion insensitivity of the VEM has been confirmed.
The new VEs have also been tested in the nearly incompressible limit. Also in this case, the new self-stabilized VEs with additional moment DOFs have provided superior performances, exhibiting an almost completely locking-free behavior. The new VEs, self-stabilized without additional moment DOFS, have instead shown performances very similar to those of the standard stabilized VEM, including a substantial deterioration of the convergence rate in the nearly incompressible limit. However, we highlight that a complete  theoretical analysis of the stabilization free VEs in a general framework (arbitrary polynomial order and arbitrary polygonal meshes) is not available. Nonetheless, for the present Hu-Washizu approach a stability and convergence study, concerning the case k = 1 on quadrilateral meshes, can be found in [27].
We conclude by noticing that possible interesting future developments may consider the extension of the proposed mixed variational formulation to three-dimensional VEs, to elastoplasticity and to elastodynamics.

Appendix A Computation of hourglass matrix H
To define the hourglass modesû H and the hourglass matrix H in (28) it is possible to proceed as follows (see, e.g [12,16]). According to the natural approach proposed by Argyris [15], the element nodal displacementsû can be expressed as a linear combination of n R rigid body modes and n D natural or straining modes, through a non-singular matrix T. In the case that n H hourglass modes are also allowed by the element kinematics, the vector of displacement parametersû can be expressed as: where T D+R u is a matrix of n D + n R columns, each one representing an independent deformation or rigid mode, T H u is a matrix of n H columns, each one representing an independent hourglass mode,p D+R u are natural parameters andp H u are hourglass parameters defining the amplitude of the corresponding hourglass mode. Deformation modes and hourglass modes are taken orthogonal to each other, i.e.: The orthogonality between deformation/rigid body modes and hourglass modes implies also that: Making use of this orthogonality property and of the decomposition (A1) ofû, one can write: From the previous expression it is possible to compute the natural parameters associated to deformation and rigid body modes: Substituting this expression in (A1) and solving forû H = T H up H u , one obtains: where the hourglass matrix H is defined as: and I denotes the n u × n u identity matrix. As it can be easily verified, one also has H T H = H. The problem of computingû H is then reduced to the construction of the hourglass matrix H and, hence, to the computation of the transformation matrix T D+R u associated to deformation and rigid body modes. In the standard VEM of order k, the approximate displacements locally contain all the polynomials of degree at most k; accordingly, the strains are initially modelled by projecting symmetric gradients of displacements onto polynomials of degree at most k −1. However, for arbitrary polygons the displacements are richer and other (typically non-polynomial) functions are required to cope with inter-element continuity. These additional functions are responsible for a surplus in the rank deficiency of K c with respect to the standard value 3. On the basis of this observation, the number of hourglass modes is equal to n u − 2n k , n k being the number of parameters required to define a complete polynomial of order k (35). Of course, n u − 2n k depends on the number of the element vertices, see (36). As a consequence, the approximate displacement without hourglass modes can be expressed as: where: Considering that: the following equality holds: Since by definition N u i (ξ i ) = 1, if ξ i are the coordinates of the ith boundary node and N u i is the corresponding shape function, considering (A11), one has, e.g., at the edge node 1: The matrix product at the right hand side of (A12) returns the first two rows of the matrix T D+R u . Repeating the same procedure for all the other DOFs, one obtains the whole n u × 2n k matrix: where: The first 2kN V rows contain the scaled monomials evaluated in correspondence of DOFs on the element boundary. The last k(k − 1) rows contain the following terms: with i = kN V + 1, kN V + 2, . . . , n u 2 and j = 1, 2, . . . , n k . All the other terms in the last k(k − 1) rows are zero.
Summarizing, the matrix T D+R u can be computed evaluating N k (ξ ) in correspondence of the element DOFs. Once T D+R u is computed, the hourglass matrix H can be computed by means of (A7) and finally K s can be evaluated.

Appendix B Computation of equivalent nodal forces vector
The procedure for the construction of the element equivalent nodal force vector is different for the cases k = 1 and k ≥ 2. Even though the idea is the same for k = 2 and k > 2, for the sake of clearness these two cases will be analyzed separately. Note that the term "nodal forces" is used even though some of the components of the force vector are not "nodal" when moments DOFs are present.
Looking at (20), the term F p is fully computable since N u is explicit on the boundary of the domain, where it is described by polynomial functions of degree k. Therefore, only the evaluation of the term F b will be discussed below.
In the case k = 1, the equivalent nodal force vector associated to the body forces b is computed in an approximate way, exploiting an integration rule associated to the element vertices. Basically, the load vector is uniformly distributed to each DOF of the virtual element. Let us introduce the 2 × 1 vector f: whose components are The equivalent nodal forces vector due to b is computed as: where N V denotes the number of vertices of the virtual element.
In the case k = 2, the basic idea is the same as for the case k > 2, namely that of projecting the body force vector b onto the space P k−2 of polynomials of degree up to k − 2.
The computation results to be very simple due to the presence of the internal moment DOFs. In the specific case k = 2, the vector b is projected onto a vector of constantsb h , that can be derived from the condition: Assembling all the terms from j = 1 to j = n k−2 , one gets the algebraic system: from which the vectorb h containing the parameters of the polynomial expansion can be derived. The symmetric 2n k−2 × 2n k−2 matrix Q has the following structure: The 2n k−2 × 1 vectorsb h and f contain all the 2 × 1 contributions: Finally, the element vector of nodal forces equivalent to b can be computed from: where the matrices in brackets are the n k−2 moments of the displacement shape functions N u , once again computable exploiting (39), (40), (41) and (42). Due to the fact that the moments of the displacement shape functions on the element boundary are zero, the only non-zero components of F b are those related to moment DOFs.