First-order VEM for Reissner–Mindlin plates

In this paper, a first-order virtual element method for Reissner–Mindlin plates is presented. A standard displacement-based variational formulation is employed, assuming transverse displacement and rotations as independent variables. In the framework of the first-order virtual element, a piecewise linear approximation is assumed for both displacement and rotations on the boundary of the element. The consistent term of the stiffness matrix is determined assuming uncoupled polynomial approximations for the generalized strains, with different polynomial degrees for bending and shear parts. In order to mitigate shear locking in the thin-plate limit while keeping the element formulation as simple as possible, a selective scheme for the stabilization term of the stiffness matrix is introduced, to indirectly enrich the approximation of the transverse displacement with respect to that of the rotations. Element performance is tested on various numerical examples involving both thin and thick plates and different polygonal meshes.

Notwithstanding the importance of the shear deformable plate in common structural applications, at best of authors knowledge, very few virtual elements for the analysis of Reissner-Mindlin plates have been presented in the literature. Beirão da Veiga et al. proposed in [23] a virtual element that assumes the shear strains and the transverse displacement as independent variables while the rotations are obtained by post-processing. This stratagem allows to avoid shear locking, while considering for each vertex 5 degrees of freedom plus 1 degree of freedom per side, i.e. basically 6 degrees of freedom per vertex in a first-order approximation framework. Moreover, Chinosi [24] presented a virtual element for Reissner-Mindlin plates which makes use of the mixed interpolation of tensorial components (MITC) philosophy to overcome shear locking phenomenon. In a first-order approximation, for each vertex/side this approach considers 3 rotational degrees of freedom, 1 transverse displacement degree of freedom, and 1 shear stress degree of freedom, i.e. basically 5 degrees of freedom per vertex.
The present paper aims at proposing simple and efficient first-order virtual elements for Reissner-Mindlin plates based on a standard displacement-based variational formulation, assuming only the transverse displacement and the rotations as independent variables in the vertexes, i.e. saving degrees of freedom with respect to existing approaches. The consistent term of the stiffness matrix is developed assuming uncoupled polynomial approximations for the generalized strains, possibly adopting different polynomial degrees for bending and shear parts. To mitigate potential shear locking in the thin-plate limit while keeping the formulation as simple as possible, a selective stabilization scheme is used for the stabilization term of the stiffness matrix. In particular, the idea is to indirectly enrich the approximation of the transverse displacement with respect to that of the rotations by assuming the stabilization order for the transverse displacement higher than the one used for the rotations. Numerical tests are conducted on a number of polygonal meshes, involving thin and moderately thick plates under different loading and support conditions to test the performance of the method.
The paper is organized as follows. The basic plate equations and the displacement-based variational formulation are presented in Sect. 2. Section 3 is devoted to the description of the VEM formulation. The main assumptions with regards to the calculation of the consistent and stabilization terms of the stiffness matrix are discussed in Sect. 4. Numerical results are presented in Sect. 5. Some final considerations end the paper.

Basic equations
A homogeneous plate P with mid-surface Ω and constant thickness h is considered, resulting P = Ω × h. A Cartesian reference coordinate system (O, x, y, z) is introduced so that x, y describe the plane of the mid-surface Ω, while the zaxis is in the thickness direction, with −h/2 ≤ z ≤ h/2, as illustrated in Fig. 1.
The Reissner-Mindlin plate theory is herein adopted to model the response of the plate. Thus, the displacement field is represented in the form: where θ x and θ y represent the rotation of the fiber orthogonal to the mid-surface of the plate in the x − z and y − z planes, respectively, while w is the transverse displacement, i.e. the deflection, of the plate. The two rotations are organized in the vector θ = θ x θ y T , useful for the next developments.
The generalized strains for the Reissner-Mindlin plate theory are the bending curvatures χ = χ x χ y χ xy T and the shear strains γ = γ xz γ yz T , defined as: where the bending and shear compatibility operators are introduced as: The equilibrium equations of the plate are: where vectors M and S collect the bending moment and shear resultants, respectively: while q and m = m x m y T are the prescribed transverse loads and couples acting on Ω.
A linear elastic response is assumed for the plate; thus, the constitutive equations are written in the form: where C b and C s are the matrices of bending and transverse shear moduli. In particular, for the isotropic case, the elasticity matrices are: being E the Young's modulus, ν the Poisson's ratio and k a correction factor to account for non-uniform distribution of shear stresses through the thickness. Because of the homogeneity of the plate in the thickness, it is set k = 5/6. The equilibrium problem of the linear elastic shear deformable plate is ruled by Eqs. (2), (4) and (6) in Ω, together with appropriate boundary conditions on ∂Ω.
The variational formulation of the above problem in terms of generalized displacements w and θ stems from the virtual work principle: ∀ δw ∈ U 0w and δθ ∈ U 0θ (9) where U w and U θ are the spaces of the admissible generalized displacements, U 0w and U 0θ the spaces of their variations and δ ext the virtual work done by external (body and boundary) loads.

VEM formulation for Reissner-Mindlin plates
To solve the problem in Eq. (9) through the VEM scheme, the domain Ω is discretized by means of nonoverlapping polygons with straight edges. Each polygon can be denoted with Ω E , while ∂Ω E indicates its boundary. In the VEM, the approximated generalized displacements in the element interior Ω E , denoted as w h and θ h in the following, are assumed to be not explicitly known and, thus, are referred to as virtual.
On the contrary, the approximated generalized displacements on the element boundary ∂Ω E , denoted in the following bỹ w h andθ h , are assumed to be explicitly known and written as: where N V w and N V θ are the matrices of the approximation functions on ∂Ω E andṼ is the vector collecting the generalized displacements of the nodes on the element boundary. VectorṼ has mk w + 2mk θ components, where m is the number of vertexes of the polygonal element, and k w and k θ are the degrees of the polynomial representations assumed over the element boundary forw h andθ h , respectively.

Projection operator and consistent term
As neither w h and θ h nor their gradients are explicitly computable in the element interior, the key aspect of the VEM consists in the introduction of a projection operator, which represents the approximated generalized strains associated with the virtual generalized displacements, called projected generalized strains and denoted by χ P and γ P in the following. Accordingly, given the virtual displacements w h and θ h , χ P and γ P can be defined as the unique functions that minimize (11) where subscript norm indicates the type of norm used. Following the approach proposed in [25], an energy norm can be used and, thus, the minimization of Eq. (11) leads to: find where P p b (Ω E ) and P p s (Ω E ) denote the spaces of polynomials of degrees (up to) p b and p s defined on Ω E , and are the spaces of the approximation functions adopted for the representation of χ P and γ P , respectively, which are assumed as: N P b and N P s being the matrixes of the polynomial approximation of the consistent generalize strains of the plate andε is the vector of the coefficient of that approximation. Using Eqs. (13) and (2), Eq. (12) can be rewritten as: Integrating by parts and using Eq. (10) yields: where and being n x and n y the components of the unit outward normal on ∂Ω E . The last three integrals in the r.h.s. of Eq. (15) collect the moments of the virtual generalized displacements (up to order p s − 1 for w h and up to order max( p b − 1, p s ) for θ h ) and their evaluation would require w h and θ h explicitly known in the interior of the element. Note that this differs from the plane elasticity case, in which only moments up to k − 1 arise. To overcome this problem, in the VEM, the moments of the virtual displacement are assumed as internal degrees of freedom of the element, in addition to the external ones associated with the nodes on the element boundary. Hence, denoting byV the vector collecting the local internal degrees of freedom, Eq. (15) can be rewritten as: wherẽ Over the generic element, using Eqs. (13) and (18), the bilinear form in (9) can be written in terms of projected strain as: where In Eq. (20), the consistent part of the stiffness matrix K c can be recognized, which takes the form: where

Stabilization term
A suitable stabilizing term may be needed to assure the rank sufficiency of the stiffness matrix and avoid zero-energy modes. Indeed, the presence of a stability term is standard for VEM [9][10][11]26,27], even though some recent advancements [13,28] highlighted the possibility to device formulations which could avoid the need for stabilization. Denoting by s w and s θ the stabilization orders for w and θ , respectively, the generalized displacements within the virtual element can be written both in terms of vector-valued polynomials up to degree s w and s θ , respectively, and also in terms of the virtual basis of the admissible generalized displacements of the element. In other words, the stabilization orders s w and s θ represent the maximum degree of the polynomials which can be found in the virtual fields of w and θ , respectively. By way of example, s w = 1 means that all quadratic and higher order terms are penalized in the transverse displacement virtual field. Accordingly, D w and D θ are set as the matrices associated to the change of basis from the polynomial space to the virtual functions space, i.e. the components of the matrices D w and D θ are given by the evaluation of the polynomials on the degrees of freedom w and θ , respectively, of the virtual element.
Consequently, denoting with A w a suitable collocation operator which matches the degrees of freedom related to w to the global ones, and with A θ a suitable collocation operator which matches the degrees of freedom related to θ to the global ones, the stabilization matrix can be written as: with being t w and t θ scaling factors which serve to guarantee the correct scaling of the energy with respect to the element size and material properties and are further specified in the following. Finally, the local stiffness matrix results in It should be noted that the definitions of the stabilization terms in Eq. (25) takes inspiration from the one discussed in [11]. In particular, for further convenience, it is worth to remark that such definitions of the stabilization terms imply that the choice of s w (s θ ) leads to stabilize all the virtual transverse displacement (rotations) fields of order > s w (> s θ ). In other words, only the virtual transverse displacement (rotations) fields of order ≤ s w (≤ s θ ) are free to develop within the element from the stabilization point of view.

First-order virtual elements
Basing on the VEM formulation presented in the previous section, here, the development of first-order virtual elements for Reissner-Mindlin plates is detailed. In particular, the main assumptions for the consistent term, the stabilization term and the loading term are discussed.

Consistent term
First-order approximation functions are assumed for bothw h andθ h , i.e. k w = k θ = k = 1. Then, similarly to what is usually done in plane elasticity, the polynomial degrees p b and p s of the approximations for the generalized strains are chosen related to the order k according to the compatibility equations, see Eq. (2), which suggest p b = k − 1 and p s = k. Hence, assuming uncoupled polynomial approximation of each generalized strain component, the following expressions for N P s and N P b are obtained: This choice for N P b and N P s leads to the origin of 1 local internal degree of freedom for w (i.e. its mean value) and 6 local internal degrees of freedom for θ (i.e. mean and first order moments), see Fig. 2, as it can be checked from Eq. (19). On this regard, it is worth to note that a standard static condensation [29] can be used to express the internal degrees of freedom in terms of those on the element boundary, so reducing the final degrees of freedom only to the latter, as carried out e.g. in [25].
According to the above assumptions, the virtual element space W h for the approximation of the transverse displace-ment w is made of vector valued polynomials of degree ≤ 1 on each side of the element (which are also globally continuous on the element border). The dimension of W h is: Note that there is one internal degree of freedom for w. The virtual element space h for the approximation of the rotations θ is made of vector valued polynomials of degree ≤ 1 on each side of the element (which are also globally continuous on the element border). The dimension of h is: Note that there are three internal degrees of freedom for both θ x and θ y . A schematic description of the degrees of freedom over a generic virtual element is shown in Fig. 2.

Stabilization term
As discussed in Sect. 3.2, the main task of the stabilization term is to assure the rank sufficiency of the stiffness matrix. As in the standard VEM formulation, a stabilization scheme with s w = s θ = k = 1 would suffice to guarantee that the stiffness matrix has the correct rank and no zero-energy modes are present. However, unfortunately, the resulting element would suffer shear locking phenomena in the thin-plate limit, as it is clearly highlighted in "Appendix A". As it is well-known, in displacement-based finite element formulations, locking phenomena can be alleviated by enriching the representation of the transverse displacement with respect to that of the rotations (i.e. k w > k θ ), for example using different nodal schemes for w and θ or by employing linked interpolations [30]. Although adopting a representation of w richer with respect to that of θ in Eq. (10) appears potentially doable also in the context of the virtual element formulation described in the previous section, a different and simpler strategy is herein followed. Indeed, the concept is to pursue the enrichment of the approximation of w with respect to that of θ by operating on the stabilization term, i.e. in a somehow indirect way. Particularly, the idea is to use a selective stabilization scheme as described in the following.

Rotations
Inspired by what is commonly done in two-and threedimensional elasticity, the stabilization order for the rotations θ is assumed equal to the order of the representation assumed for θ on the element boundary, i.e. s θ = k = 1. As noted in Sect. 3.2, this choice implies to stabilize all the virtual rotations fields of order greater than 1 or, equivalently, to allow the rotations fields up to order 1 (i.e. the linear rotations fields) to develop within the element. Accordingly, the  25), is chosen as: Note that, with this assumption, the rotations field within the element is governed by 6 parameters. Accordingly, the matrix D θ results: where the short symbol − q denotes the integral average of any scalar polynomial q on the virtual element domain, − x q and − y q denote the scaled moments of the scalar polynomial q with respect to x and y, respectively, on the virtual element domain, and where the i-th vertex of the element has coordinates (x i , y i ) being i = 1, 2, . . . , m.

Transverse displacement
The stabilization order for the transverse displacement w is assumed s w > s θ (being s θ = 1 in this case, see above), which means that some higher-order terms in the virtual transverse displacement field are allowed to develop within the element and, hence, w can rely on an approximation richer that that of θ . As it can be easily argued, an excessive stabilization release of the transverse displacement field can render the stabilization not effective in avoiding zero-energy modes and, hence, s w has to be set while paying attention to preserve the rank sufficiency of the stiffness matrix. Of course, the higher is the side numerousness m (i.e. the richer is the element boundary kinematics), the higher is the order of the virtual transverse displacement field that can be activated within the element in absence of stabilization. Hence, the effect that the choice of the stabilization order s w has on the virtual transverse displacement field, and, hence, on the rank sufficiency of the stiffness matrix, depends on the element side numerousness m. On this regards, it could be observed that we expect that the rank sufficiency of the stiffness matrix is preserved if there is a balance between the strain approximation assumed within the element and the number of degrees of freedom governing the representations of the rotations and transverse displacement fields. In other words, we expect the stiffness matrix to have the correct rank when the number of generalized strain modes plus the number of rigid body motions is greater or equal to the number of parameters governing the rotations and transverse displacement fields within the element (necessary condition). In our case, on one hand the number of generalized strain modes is equal to 9, see Eqs. (27) and (28), and the number of the rigid body motions is equal to 3. On the other hand, the number of the parameters governing the rotations field is equal to 6 (see Eq. 31). Therefore, the above balance can be fulfilled when the number of the parameters governing w is less than or equal to 6. In this way, the number of deformation modes plus rigid body motions is greater or equal to the number of parameters governing the generalized displacements. With this in mind, the following two strategies are proposed.
1. No stabilization is introduced for the transverse displacement if its degrees of freedom, that are m + 1 (vertex degrees of freedom plus 1 internal degree of freedom, see Sect. 4.1), are less than or equal to 6. This means that no stabilization is introduced for w if m < 6 (i.e. in triangular, quadrilateral, and pentagonal elements). Accordingly, it is herein adopted that when m < 6 the stabilization matrix in Eq. (24) becomes K s = A T θ K s θ A θ . 2. For m ≥ 6, the transverse displacement field is stabilized so to allow the development of a transverse displacement field within the element governed by 6 parameters, for fulfilling the above discussed balance. This can be easily achieved by assuming the polynomial space at the base of the construction of the matrix D w as (s w = 2): Note that, once the stabilization term is activated, it does not depend on the number of sides of the element, but only on the assumed polynomial space for D w . Accordingly, the matrix D w results: For m = 6, it should be pointed out that zero-energy modes would arise if the transverse displacement field is not stabilized.
All the above choices preserve the rank sufficiency of the stiffness matrix. Indeed, 3 null eigenvalues are found in the stiffness matrix after stabilization for any element side numerousness. More details on the stability of the virtual elements herein discussed are given in the following section.

Scaling
In standard displacement-based VEM for plane elasticity [11], the scaling of the stabilization term is classically carried out trough the scalar τ tr(K c ), where τ is a positive real number typically chosen equal to τ = 1/2, and the trace term guarantees the scaling of the energy with respect to the element size and material properties. Here, it should be pointed out that the consistent stiffness matrix for the plate problem collects terms related to both transverse displacement and rotations which, beyond having different units, scale differently depending on the plate thickness, see e.g. Eqs. (7)- (8). Accordingly, it is proposed herein to differentiate the scaling between transverse displacement and rotations terms through the factors t w and t θ , see Eq. (25), which can be defined as:

Loading term
The loading term, being k = 1, can be approximated by applying an integration rule based on element vertexes [11] as in standard VE technology. In other words, the load is evaluated in the center of gravity of the polygon and is distributed to the polygon vertexes.

Stress recovery
In this work, the stress recovery follows the classic procedure of standard displacement-based VEM [12], i.e. the constitutive law is applied directly on the projected generalized strains to obtain the generalized stresses. It should be herein pointed out that it would be possible, in general, to set up more advanced stress recovery procedures, see e.g. [16,25,31,32], also for the case of plates.

Numerical results
In this section, the performance of the proposed first-order plate virtual elements is tested numerically. In particular, a stability assessment based on stiffness matrix eigenvalues is discussed in Sect. 5.1, while two different load cases (i.e. Test 1 with clamped boundary conditions and Test 2 with simply supported boundary conditions) are presented and discussed in Sects. 5.2 and 5.3, respectively. Finally, some specific considerations on triangles are discussed in Sect. 5.4.

Stability assessment
In order to assess the stability of the plate virtual elements described in Sect. 4, an eigenvalue analysis on the stiffness matrix of three different virtual elements is herein conducted for three values of plate thickness, i.e. h/L = 0.1, h/L = 0.001 and h/L = 0.00001 (being L = 1 the length of the plate). In particular, a right triangle with unit legs, a square with unit sides, and a regular hexagon with side length 1/4 are considered. The Young's modulus and Poisson's ratio of the material are set to be E = 10 3 /h 3 and ν = 0.3, respectively, while τ = 1/2 is assumed as in [11]. This stability assessment follows a classical procedure also utilized in the FEM context for Reissner-Mindlin plates (see e.g. [30,33]), that consists in comparing the eigenvalues by changing the thickness of the plate while keeping constant the bending stiffness. Stable elements will show a substantially constant ratio between the eigenvalues themselves while reducing thickness, whereas unstable elements, i.e. characterized by locking phenomena, will be characterized by one or more eigenvalues which show a significant change with respect to the others, highlighting a solution which changes while reducing the thickness. Table 1 collects the non-null eigenvalues after static condensation of the stiffness matrix for the three considered virtual elements, computed by keeping constant the bending stiffness while decreasing the thickness.
The obtained results reveal two fundamental points: • the proposed strategy for the stabilization works fine, as in all cases only 3 null eigenvalues are found ensuring the rigid body motions and, in the meantime, all the elements are free from spurious (zero energy) modes; • all the elements present almost constant eigenvalues for significant different values of the side-to-thickness ratio, confirming that elements do not suffer from shear locking.
It can be also remarked that, while for triangle and square no-stabilization is introduced for the transverse displacement field w, the stabilization according to Eq. (33) is adopted for w in the hexagon. Finally, this first numerical investigation ensures the correct features of the virtual elements.

Test 1: Clamped boundary conditions
A load case with clamped boundary conditions is herein considered on a unit square domain Ω = (0, 1) 2 , by assuming Dirichlet's conditions over the entire boundary. Two values of plate thickness are considered, i.e. h/L = 0.1 and h/L = 0.001, being L = 1 the length of the plate. The Young's modulus and Poisson's ratio of the material are set to be E = 1.0 and ν = 0.3, respectively, while τ = 1/2 is assumed as in [11]. Numerical tests are performed on the meshes shown in Fig. 3: structured regular quadrilaterals (QUAD, Fig. 3a), structured distorted quadrilaterals (QDIS, Fig. 3b), convex/concave quadrilaterals (RHOM, Fig. 3c), structured regular hexagons (HEXA, Fig. 3d), random convex polygons (POLY, Fig. 3e), convex/concave 6-vertex polygons (WEBM, Fig. 3f), structured regular 6-node triangles (TRI6, Fig. 3g), and structured regular 3-node triangles (TRI, Fig. 3g). In particular, this last case is discussed specifically in Sect. 5.4. In order to test the convergence properties of the approach, the following relative error in discrete L 2 -like norm (based on the vertex values) is introduced to have a global estimation on the generalized displacement error: where N tot is the global number of vertexes, w E is the exact transverse displacement, θ x E and θ y E are the exact rotations, w h is the approximated transverse displacement, and θ x h and θ y h are the approximated rotations. Also, a classical relative error in energy norm is introduced to have a global estimation on the generalized stress error: In the following, "log" denotes logarithm with base 10, and "dofs" denotes the total number of degrees of freedom that results after static condensation. Test 1 refers to a square plate clamped on the whole boundary, for which the analytical solution [34] is explicitly known. Indeed, by assuming the external load as: with m x (x, y) = 0 and m y (x, y) = 0, the analytical solution in terms of generalized displacements results as: 5 show the displacement and stress error convergence curves, respectively. Particularly, the meshes collected in Fig. 3 are tested for both h/L = 0.1 (see Figs. 4a and 5a) and h/L = 0.001 (see Figs. 4b and 5b). As it can be noted, the convergence rates are in agreement with the expected ones for both displacement and stress errors. Addi-tionally, it should be pointed out that the errors do not show significant differences when passing from a moderately thick plate (h/L = 0.1) to a thin plate (h/L = 0.001), thus showing a locking-free nature of the method.

Figures 4 and
Finally, a direct comparison of the numerical results versus the ones presented in [23] is collected in "Appendix B". As a result, the present approach shows a slightly better accuracy with respect to the results presented in [23], for both moderate thick and thin plates. This outcome appears remarkable, as the present approach is characterized by 3 degrees of freedom per vertex, while 6 degrees of freedom per vertex are supposed in [23].

Test 2: Simply supported boundary conditions
Test 2 refers to a simply supported square plate. Here, the same material of Test 1 is adopted, as well as the same meshes   Fig. 3) and error definitions of Sect. 5.2. The analytical solution is explicitly known also for Test 2. Indeed, by assuming the external load as with m x (x, y) = 0 and m y (x, y) = 0, the analytical solution in terms of generalized displacements results as: where the values of the coefficients W , R x , and R y for the present case are collected in Table 2. Figures 6 and 7 show, respectively, the displacement and stress error convergence curves. Analogously to Sect. 5.2, the meshes collected in Fig. 3 are tested for both h/L = 0.1 (see Figs. 6a and 7a) and h/L = 0.001 (see Figs. 6b and 7b). As it can be noted also in this case, the convergence rates are in agreement with the expected ones for both displacement and stress errors, and, analogously to the previous load case, the errors do not show significant differences when passing from a moderately thick plate (h/L = 0.1) to a thin plate (h/L = 0.001), thus confirming the locking-free nature of the method.

Considerations on triangles
The case of triangles m = 3 is herein discussed separately as it presents some peculiar critical aspects. It should be noted that it could be somehow expected that the selective sta- bilization scheme presented in Sect. 3.2, which indirectly aims at enriching the approximation of w with respect to θ , can be not optimal for triangular elements given their poor kinematics. In fact, triangles are characterized by the poorest transverse displacement field that can be activated in absence of stabilization on w. This can be noted in Figs. 8, 9, 10 and 11, which collect the convergence curves of displacement and stress errors for Test 1 and Test 2. Indeed, by setting τ = 1/2 as in the previously presented results, convergence rates show unsuitable agreement with the expected ones for the thin plate case (h/L = 0.001), although the moderately thick case (h/L = 0.1) does not show any specific drawback with this assumption. It is interesting to observe that the performance of triangles in the thin-plate limit could be recovered by the reduction of the coefficient τ which rules, in this case, the scaling of the stabilization matrix K s θ . Indeed, the weight of the stabilization matrix with respect to the local stiffness matrix decreases while reducing τ . With this in mind, Figs. 8, 9, 10 and 11 show also the convergence curves of displacement and stress errors for Test 1 and Test 2 by varying τ , i.e. for τ = 1/7, τ = 1/10, and τ = 1/20, in both h/L = 0.1 and h/L = 0.001 cases. Moreover, for completeness and in the same fashion of Tables 2, 3 collects the non-null eigenvalues for the triangle with τ = 1/20 for different thicknesses.
By focusing on the thin plate case h/L = 0.001, Figs. 8, 9, 10 and 11 show that the reduction of τ produces a positive effect in terms of convergence rates for both Test 1 and Test 2, i.e. displacement and stress error convergence rates tend to the expected ones while reducing τ . At the same time, it is also observed a progressive reduction in accuracy in the displacement error while decreasing τ , although no accuracy loss is observed in stress error curves. Anyway, it should be pointed out that if the rotations field is not stabilized, i.e. no stabilization is considered over a triangular element, zero-energy modes would arise in the local stiffness matrix. Summing up, in all cases, the convergence performance of the triangular elements is progressively enhanced while reducing τ .

Conclusions
In this paper, a first-order VEM for Reissner-Mindlin plates has been presented, based on a standard displacement-based variational formulation assuming transverse displacement and rotations as independent variables.
The consistent term of the stiffness matrix has been developed assuming uncoupled polynomial approximations for the generalized strains. Particularly, a linear polynomial approximation has been assumed for the shear response, whereas a constant approximation has been adopted for the bending one.  With the purpose of mitigating shear locking in the thinplate limit while keeping the element formulation as simple as possible, a selective scheme has been introduced in the calculation of the stabilization term of the stiffness matrix. In particular, the stabilization order for the transverse displacement is assumed higher than the one used for the rotations, so to achieve an indirect enrichment of the approximation of the transverse displacement with respect to that of the rotations. On the one hand, standard first-order stabilization has been adopted for the rotations (i.e. rotations fields up to order 1 are allowed). On the other hand, no stabilization has been introduced for the transverse displacement when the polygon side numerousness is lower than 6, while a quadratic stabilization order has been introduced for the transverse displacement when the polygon side numerousness is greater than or equal to 6. Following this choices, the resulting stiffness matrices have been found to be stable in the thin-plate limit and to preserve rank sufficiency.
A number of numerical tests have been conducted on several polygonal meshes, involving thin and moderately thick plates under different loading and boundary conditions (e.g. clamped and simply supported). Very good agreement of the convergence rates with the expected ones has been found for both thin and thick plates for meshes with polygons with more than 3 sides, thus showing the locking-free nature of the method. As expected, the case of triangles appeared not fully optimal given their poor transverse displacement field even when no stabilization is introduced. Anyway, the performance of triangles has been found to recover while reducing the weight of the stabilization matrix in the local stiffness matrix.
To conclude, the main objective of the paper, i.e. introducing simple and efficient virtual elements for Reissner-Mindlin plates, appears achieved as the proposed approach allows to save degrees of freedom with respect to existing approaches while keeping comparable (or even better) results. Further future developments could regard new strategies to develop the stabilization term of the stiffness matrix to optimally treat also triangles.
Funding Open access funding provided by Alma Mater Studiorum -Università di Bologna within the CRUI-CARE Agreement.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.

Appendix B
In this appendix, a results comparison is conducted between the present approach and the approach presented in [23]. Particularly, the comparison is pursued in terms of error on the transverse displacement, as it is the only primary variable that the two approaches have in common. Indeed, in [23] the rotations are obtained in post-processing. Accordingly, the following relative error e w in discrete L 2 -like norm (based on the vertex values) is introduced according to [23] to estimate and compare the error on transverse displacement: where N Ω E is the number of vertexes of the element Ω E , and |Ω E | is the area of the element Ω E .
The comparison of the transverse displacement error convergence curves between the present approach and the approach in [23] is shown in Fig. 14 for h/L = 0.1 and h/L = 0.001 (this time, given the data available in [23], in terms of the refinement parameter, i.e. the maximum element diameter of the mesh). Test 1 with clamped boundary conditions and the convex/concave 6-vertex polygons mesh (WEBM, see Fig. 3f) is herein considered.
As it can be noted (Fig. 14), the present approach shows a convergence rate in agreement with the expected one and the one shown in [23]. In addition, the present approach shows a slightly better accuracy with respect to [23] for both moderate thick (h/L = 0.1) and thin (h/L = 0.001) plates (Fig. 14). This appears remarkable, as the present approach is characterized by 3 degrees of freedom per vertex, while 6 degrees of freedom per vertex are supposed in [23]. Comparison of the transverse displacement error convergence curves between the present approach and the approach in [23] for a h/L = 0.1, and b h/L = 0.001. The horizontal axis of the graphs goes from coarse (left) to finer meshes (right). In this case, the expected convergence rate is 2 as the curves are in terms of the refinement parameter