On two simple virtual Kirchhoff-Love plate elements for isotropic and anisotropic materials

The virtual element method allows to revisit the construction of Kirchhoff-Love elements because the C1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C^1$$\end{document}-continuity condition is much easier to handle in the VEM framework than in the traditional Finite Elements methodology. Here we study the two most simple VEM elements suitable for Kirchhoff-Love plates as stated in Brezzi and Marini (Comput Methods Appl Mech Eng 253:455–462, 2013). The formulation contains new ideas and different approaches for the stabilisation needed in a virtual element, including classic and energy stabilisations. An efficient stabilisation is crucial in the case of C1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C^1$$\end{document}-continuous elements because the rank deficiency of the stiffness matrix associated to the projected part of the ansatz function is larger than for C0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C^0$$\end{document}-continuous elements. This paper aims at providing engineering inside in how to construct simple and efficient virtual plate elements for isotropic and anisotropic materials and at comparing different possibilities for the stabilisation. Different examples and convergence studies discuss and demonstrate the accuracy of the resulting VEM elements. Finally, reduction of virtual plate elements to triangular and quadrilateral elements with 3 and 4 nodes, respectively, yields finite element like plate elements. It will be shown that these C1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C^1$$\end{document}-continuous elements can be easily incorporated in legacy codes and demonstrate an efficiency and accuracy that is much higher than provided by traditional finite elements for thin plates.


Introduction
The necessity to accurately and efficiently model plates has a long history based on the fact that plates are installed as structural members in almost every building, in many machines and in airplanes. This has led early on to the development of plate elements in most discretisation schemes. In the finite B P. Wriggers wriggers@ikm.uni-hannover.de B. Hudobivnik hudobivnik@ikm.uni-hannover.de O. Allix olivier.allix@ens-paris-saclay.fr 1 element method associated development started more than 50 years ago with the work of [8,17,24]. It was obvious that general forms of finite elements were not able to fulfill C 1 -continuity when using Kirchhoff-Love plate theory as starting point. First elements that met C 1 -continuity were developed within the TUBA series of elements in [4]. Composite elements, consisting of four triangles were provided by [12,18] which then allowed to use such composite element as general quadrilateral other formulations for strictly rectangular elements can be found e.g. in [14]. All these formulations have the disadvantage that besides deflection w and rotations w ,x , w ,y additional higher order kinematical quantities, e.g. the curvatures w ,x x , w ,yy , appear as nodal unknowns which impedes the use of such element in an engineering environment.
The latter complication was circumvented by the introduction of Reissner-Mindlin elements, starting with the work of [20,39]. Due to some disadvantages in the thin plate limit, many different variants were considered in the following years, see the textbooks e.g. [19,30,31]. One of these variants is based on the discrete, point-wise formulation of C 1 -continuity, which leads to the discrete Kirchhoff triangles (DKT), see [6], and quadrilaterals (DKQ), see [7].
As shown in the work of [13], the virtual element method (VEM) can be applied to develop C 1 -continuous Kirchhoff-Love plate elements. VEM actually has the advantage that no higher order kinematical quantities have to be introduced and thus these elements can be easily integrated in engineering codes which is due to the simplicity of the definition of the ansatz functions, which are defined only at the element boundary. Hence the virtual element method allows to revisit the construction of Kirchhoff-Love elements. This was illustrated for two simple plate elements in [16] which includes error estimates. [26] employed the virtual element methodology to dynamic plate problems and also showed that these elements can be used for the buckling analysis of plates, see [25,27]. At the same time other virtual plate elements using Reissner-Mindlin kinematics were developed, see [11] and non-conforming discretisations were considered in [3] to model plate problems with VEM.
Here we present two simple VEM elements suitable for conforming Kirchhoff-Love plates which were already discussed for isotropic materials in [16]. We extend this formulation to isotropic and anisotropic material responses, provide a detailed description of the matrix forms and new stabilisations. The reason for restricting ourselves to low order ansatz functions was motivated by the fact that virtual elements of low order are the most suitable choice for non-linear extensions, see e.g. [1,10,15,38].
The main difficulty in VEM concerns the stabilisation parts of the stiffness matrix. This question is even more crucial in the case of C 1 -continuous elements than in the case of C 0 -continuous elements, because the rank deficiency of the stiffness matrix associated to the projected part of the ansatz function is larger. Therefore this paper aims-besides providing a formulation in engineering terms-at introducing different possibilities for the stabilisation of virtual plate elements. Compared to the approach in [26], this work aims to improve stabilisation and additionally introduces an energy based stabilisation which is based on a sub-triangulation of the VEM elements, as was proposed in [38] for hyperelastic solids. In the latter stabilisation, different simple conforming and non-conforming plate elements, see [6,29], are compared regarding their contribution to the stiffness matrix, their precision and efficiency of the resulting VEM element.
Numerical examples will be considered to study the convergence behaviour of both developed elements and to illustrate their applicability to engineering problems. We first deal with problems with analytical solution and then examine the case of rhombic plates where its is known that higher degree elements are preferable, see [5]. Examples including anisotropic material behaviour show the applicability of the method to composite plates. Since the developed virtual plate elements only introduce deflections and rotations as nodal unknowns, despite being C 1 -continuous, they can be easily incorporated in finite element software when they are reduced to triangular elements with three/six nodes or quadrilateral elements with four/eight nodes. Based on classical benchmark examples it can be concluded that these reduced virtual element have a much higher efficiency and accuracy than comparable traditional finite plate elements.

Description of the plate and of the constitutive relation
The geometry of the plate is described through its mid surface m with thickness, 2h, which will be assumed to be constant. A point X within the plate is given by the orthogonal projection X m onto the mid surface and its coordinate Z along the normal E Z of m (see Fig. 1).: where (E X , E Y , E Z ) is a orthonormal direct cartesian basis.
In the Kirchhoff-Love theory the displacement u(X m , Z ) of point X can be described by the deflection w of the plate where ∇ m , simply denoted by ∇, is the gradient operator with respect to X m and w belongs to H 2 ( m ). Therefore the rotation θ of the normal segment is equal to The in-plane strain ε associated with the displacement u(X m , Z ) is given by where the curvature operator χ is equal to In what follows we will assume that the plate is symmetrical with respect to its mid-surface. Introducing the classical Hooke's law under plane stress condition yields the constitutive tensor K cp . With that one obtains a constitutive relation linking the moment of the in-plane stress to the curvature where a defines an integral of quantity a over the thickness t = 2h of the plate For a plate that is homogeneous through the thickness the integral (6) yields The local bilinear form associated with the Kirchhoff-Love model can be written as where the moment follows from (8) for an isotropic homogeneous plate with Young's modulus E and Poisson ratio ν With the definition of the bending stiffness the bilinear form a(w, v) of the Kirchhoff-Love theory over which leads to the bilinear form characterising the internal virtual work of the Kirchhoff-Love problem.
Alternatively the strain energy density can be introduced In what follows the symbolˆis associated to the Voigt notation. For example, for the moment and the curvature one introduceŝ The constitutive relation may be expressed through the matrix D that is defined with (11) bŷ Thus the strain energy in (14) can be rewritten in Voigt notation as For a laminated plate it is necessary to model the anisotropic behaviour, see e.g. [32]. A laminate consists of a stack of composite plies that have a fiber direction which points in the direction of the angle φ to the X -axis. Let us note that we consider in this paper only stacking sequences that are symmetric with respect to their mid-plane in order to avoid coupling between tension and bending. For each ply k the orthotropic constitutive matrix may be defined with respect to the orthotropic basis (the notation¯refers to the orthotropic basis, while not using a bar refers to the global basis): where E 1 relates to the stiffness in fiber direction and the value E 2 is the stiffness perpendicular to the fiber direction. The Poisson ration of the ply is given by ν 12 and the shear modulus is G 12 . The Poisson ratio ν 21 = E 2 E 1 ν 12 is a dependent variable. This matrix has to be transformed to the X − Y coordinate system by the transformation with (20) Now an integration over the thickness has to be performed considering all n pl plies. This leads, analogous to (7), by employing a sum over all plies to the global material stiffness matrixD Ĝ where the integration is executed for each ply over the ply thickness h k = Z k − Z k−1 . Now the strain energy for the laminated plate is given by TD Gχ (22) and the moments follow fromM =D Gχ .

Basis of the virtual element method for Kirchhoff plates
Most of the mathematical theory of the virtual element method (VEM) for Kirchhoff plates was established in [13]. The theory has been extended to vibrations in [26] and to buckling analysis notably in [25,27]. We briefly recall some aspects of the theory useful for this paper. The domain m is partitioned into non-overlapping polygonal elements which need not be convex but basically star-shaped. Following the framework of the VEM, the space of a virtual Kirchhoff-Love element has to be C 1 continuous at the inter-element level. A zoo of elements may be defined by four integers, see [13]: • The first one, r , defines the regularity of the polynomial chosen for the deflection w on the contour of each element (and thus of its tangentiel derivative w, t ), • the second one, s, defines the regularity of the polynomial chosen for its normal derivative w, n , • the third one, m, defines the polynomial order of the bilaplacian inside the element and • the fourth one, k, possibly the most important, corresponds to the order of the method.
Not all possible choices of the quadruplet r , s, m, k are acceptable.
In what follows E denotes an element of a given mesh and e any of the corresponding edges of the element. The space of the polynomial P m is given by the degree, less or equal to m for the corresponding geometrical entity V h = {w ∈ H 2 ( m ), 2 w ∈ P m (E), w, t ∈ P r (e), w, n ∈ P s (e)} (23) Among all relevant choices of r , s, m, k some appear more natural than others. Since the local equation of equilibrium involves the fourth order derivative of w, the choice of m = k − 4 is the one minimising the number of internal d.o.f.s for a given order. We adopt this condition which avoids the introduction of internal nodes for the chosen elements studied in this paper. We make use of the convention that P m (E) = {0} for m < 0. The condition 2 w ∈ P m (E) is fundamental for the element to be uni-solvent that is ensuring the uniqueness of w inside the element knowing its degrees of freedom. The order of the method is ensured by the following consistency conditions where A h is the approximation of the bilinear form associated with the strain energy of the model: A method used in many papers on VEM to the recover stability along with consistency, see [9] for solids and [13] for plates, is the following: Let k be the energy projector onto the set of polynomials of degree k where k corresponds to the degree of projection of k (w) which can uniquely be prescribed knowing the degree of freedom of the element.
k is defined as follows: Here A h ( p, q) follows from (12) and (13). Generally the projection in (25) is not sufficient for defining a proper stiffness matrix because it leads to a rank deficient matrix. Using the following property of k : we note that the first term is computable as a function of the degrees of freedom of the virtual elements but the second one is not because (v h , w h ) are not defined by the ansatz (23) within the element. The stability of the method may nevertheless be ensured by replacing by an equivalent partial bilinear form which scales as the original one and is defined by the value of the element unknowns.
Let us note that, for a given element, there exist several ways to complement the stiffness matrix following from (25) in a stable manner. All possible stable methods behave asymptotically in the same manner but may lead to rather different results for coarse meshes. This is precisely one of the issues addressed in this paper.

Formulation of the chosen virtual element method for Kirchhoff-Love plates
By using the specific strain energy (9) the virtual element for the Kirchhoff-Love theory can be derived based on where X i are the coordinates of possible loaded points in the domain and [|M nt (X i )|] corresponds to the jump of M nt at X i computed following the curvilinear abscissa of the frontier of the domain. Let us denote by V 1 h the space corresponding to element 1 with k = 2 One obtains from (28) the consistency condition (24) for element 1 Since w, n and w are defined on the boundary the left hand side of the equation is known and the consistency condition is easily enforced. It leads to three conditions. Let us now denote the space V 2 h corresponding to element 2. For element 2, k = 3 and thus , w, t ∈ P 2 (e), w, n ∈ P 2 (e)} = {w ∈ H 2 (E), 2 w = 0, w, t ∈ P 2 (e), w, n ∈ P 2 (e)} . Again relation (28) yields the consistency condition ∀ p ∈ where w and w, n are known at the boundary. Hence the left hand side of the equation is known and the consistency condition is easily enforced.
Let us now consider the projector k . For element 1, 2 is defined as Thanks to the consistency conditions the right hand side of the first relation defining 2 is known. The second and third condition give rise to three additional conditions which completely define 2 ( p). It depends on 6 scalars as function of the degrees of freedom of the element. For the same reason, when deriving element 2, the projector 3 , involving 10 constants defining 3 ( p), is also perfectly defined.
In what follows we denote the number of elements by n T . For each element n V denotes the number of vertices which is equal to the number n E of edges, n N is the number of nodes and n D describes the number of degrees of freedom. The number of unknowns associated with the projection is given by n . Table 1 summarizes the data for element 1 and element 2 and provides the rank deficiency n R associated with the use of the sole projection for the definition of the stiffness matrix of the element. Let us note that the normal rank deficiency of the stiffness matrix in bending should be 3, due to rigid body modes.
In both cases no internal degrees of freedom are needed which simplifies possible extensions of the method to nonlinear cases. Nevertheless this choice implies some approximation when computing the virtual power of external body load for the element of V h .
The ansatz for both elements is based on Hermite cubic functions for w, ensuring the continuity of the normal component of the gradient along the edges. In each case, the cubic ansatz for the deflection w along the element edge is given for a boundary segment k of the virtual element, defined by the local nodes (1)-(2) (respectively (1)-(2)-(3)) with three unknowns for w and θ per vertex, see Fig. 2.
For element 1 the gradient of w along the edge is quadratic while its normal derivative is only linear. Element 1 is of order 2.
To construct a quadratic ansatz for the gradient along the edge one has to introduce an additional unknown per edge, for example the value of θ t at the mid node of the edge. This element is denoted element 2 ( Fig. 3).
It is interesting to compare the second element with the first one since it is of order 3, at the price of only one additional degree of freedom per edge. A second issue studied in the paper is related to the relative effective performance of elements 1 and 2 with respect to the type of stabilisation.
The segments of the elements are indexed by k = 1, . . . , n E . In order to use the same notation for the two element formulations we adopt the following local convention for an edge k of the element. The origin of the local coordinate system is situated in the node 1 (X k = 0) and node 2 is located in X k = L k (see Fig. 2). In element 2 a node 3 is added at the middle of the edge (X k = L k 2 ) with the unknown θ t . Moreover a non-dimensional local coordinate ξ k = X k L k is introduced to defined the ansatz function of both elements.
Hence all element nodes are placed at the vertices for w and its gradient in case of element one. The discrete space of test functions on m is denoted by V h , and for a conforming approach we require that V h ⊂ V . This requirement is met by the definition of the shape or basis functions in V h . The C 1 continuous functions of an element E include (but the space is larger than that) quadratic functions.
The tangential and normal vectors (T , N) change from segment to segment. In the 2D case, normal N and tangent T are given for a segment k as where L k is the length of the segment k and (X i , Y i ) with i = {1, 2} are the nodes of the vertices defining the segment.

VEM ansatz functions
The virtual element ansatz functions are defined within the local orthonormal basis associated with an edge is denoted by N, T , E z where N points outward. To ensure the C 1 continuity requirement of the ansatz functions the values of w its normal derivative ∇w, n have to match along the edge between two adjacent elements. The following relations hold for the derivatives in normal and tangential direction With the introduced notation the cubic ansatz for the deflection along the element edge is given for both elements on a segment k by where the basis functions are defined in terms of Hermite cubic splines The derivative θ n in (36) is then provided by which yields the explicit form for element 1 and 2.
The tangential rotation (θ ht ) k is defined by the ansatz

Derivation of the L 2 projector
Following the virtual element method, we split the ansatz space into a projected part (w h ) and a remainder.
As previously discussed, a main aspect of the virtual element method (VEM) is the choice of the projector of the deformation onto a specific ansatz space. This has now to be specified for the ansatz related to the two elements discussed above.

Some useful properties
One of the advantages of the virtual element method is, as we will see in the subsequent sections, that integrals have only to be computed at the boundary. But often polynomials are given as functions over the surface E of an element. Those polynomials can always be integrated as line integral over the element boundary. By using the divergence theorem, for example for a 2D domain one can write in general This theorem can now be used for polynomial functions which yields a simplified treatement Hence any polynomial can be computed exactly over a polygon with arbitrary shape, this even holds for a non-convex virtual element.

Derivation of element 1
Since element 1 is of order 2 due to the choice of a specific ansatz for θ t , in (41) the projection has to be quadratic. It is defined at element level by where a i are unknown parameters that have to be linked to the element unknowns by the projection procedure. The curvature of the projection follows from (46) as which is constant at element level. Let us recall that the L 2 projection is defined as follows, Since χ ( p) spans with (47) the space of constant curvature this condition is equivalent to The right hand side of (49) can be exactly integrated which yields the dyadic products T ⊗ N and N ⊗ N change from segment to segment. Hence by comparing (47) and (50) the unknowns a 4 to a 6 are obtained by inspection where all contributions related to the segments have to be added. We observe that a 4 to a 6 are defined as a linear combination of the nodal values w and θ t . The constants a 1 to a 3 still have to be determined. These follow from the two last conditions in (25) The term on the left hand side in (51) yields with (46) and involves only polynomials. The right hand side of (51) can be explicitly determined as From Eqs. (53) and (51) one derives with (54) Here the last term can be simply integrated over the boundary using relations (45). Thus (55) leads to a linear mapping of a 2 and a 3 in terms of the element unknowns. The last term that has to be connected to the nodal values of the ansatz is the coefficient a 1 . Using Eq. (52) the coefficient a 1 follows directly. With (46) and (52) one writes which yields a 1 immediately since a 2 to a 6 are known by the relations (50) and (55). With this result the L 2 projection (w h ) of the virtual element is completely defined in terms of the nodal unknowns of the virtual element. We can combine (50), (55) and (56) in the mapping using (34). The mapping in (57) will not be formed explicitly since it is not necessary when using the software AceGen, see e.g. [22], to automatically derive and code the stiffness matrix associated with element 1.

Derivation of the L 2 projector for element 2
Due to the choice of the ansatz for w and θ t in (38) and (42) element 2 is of order 3. Thus the projection is cubic and defined at element level by The curvature of the projection follows with (59) as a linear function in X and Y Let us consider again the L 2 projection which is, for element 2, defined as follows, The integral on the right hand side yields with the divergence theorem Note that χ( p) is symmetric. Hence the integrand of the first term of the right hand side can also be expressed as Again the integrand of the last term can be transformed using For the cubic ansatz in (59) the last term is zero. All integrals on the right hand side can be evaluated using the ansatz functions (38) and (42). Due to the fact that divχ ( p) is constant and the normal vector N k is constant at each edge k of the element we can rewrite (64) as The first integral of the right hand side can be evaluated analogously to (55). The only difference is that now some integral terms occur where the coordinates X and Y appear. The second integral is exactly the same as the right hand side in (56).
The additional conditions (51) and (52) have to be evaluated as well. The first one yields for the integral on the left side which can be simply integrated over the boundary using relations (45). The right hand side of (51) has the same form as and equalising (66) and (67) determines a 2 and a 3 of (59) as a linear combination of the the element unknowns. The explicit form for the coefficients a 2 and a 3 is equivalent to (55). Finally the coefficient a 1 can be computed as in (56) the only difference is that now on the left side the ansatz function (59) has to be inserted.
Also for element 2 a mapping can be introduced where is the vector of all nodal unknowns of element 2. Again the local rotations associated with specific edges can be transformed to global rotations using (58). The generation of the vectors and matrices that appear in the description of both elements can be unified, see Sect. A. This is especially useful when employing the software tool AceGen, see [21,22] to generate the software code.

Construction of the virtual stiffness matrix: different strategies
The basis for the element development is the total energy which can be obtained by assembling all element contributions for the n T virtual elements, In the following we will first discuss the formulation of the element part that stems from the projection, see last Section. Furthermore, different possibilities for stabilising the virtual plate elements 1 and 2 will be considered.

Part of the virtual element due to projection
With the strain χ of the plate and the loading q we obtain the potential energy of an element with the constitutive tensor D.
Note that the strain χ is constant for element 1 which leads to a trivial evaluation of the first term. For element 2 the strain χ is linear in X and Y , hence the integration involves terms up to second order in the coordinates which again can be performed over the boundary using (45).
The projection (w h ) is known in terms of the unknowns (w , θ X , θ Y ) k at the vertices for element 1 and 2 and at the unknown rotation θ t at the mid node of element 2, see (57) and (68). These mappings can be inserted into (70) which yields for element e = 1, 2 Once the integration over the element area is carried out the potential is just a function of the unknowns of the element. Then residual vector R e and stiffness matrix K e follow from for element e.

Classical stabilisation
Within the virtual element method the curvature of the projected part of the deflection w is approximated by a constant and linear part, depending on the element type, as discussed in the last Section. A construction of a virtual element which is based purely on this projection leads to a rank deficient element, see Table 1. Thus the formulation has to be stabilised, see [13,26]. In what follows we use first the stabilisation suggested in these papers. This leads to a stabilisation operator which addresses the error at all element nodes (vertices and midpoints): Here h e is the maximum diameter of the virtual element e thus h 2 e can be interpreted as the element area e . The function w h and the projection (w h ) have to be evaluated at the vertices X i .
The formulation (73) is now extended to an integral describing the total error on the edge instead of discrete values. This new stabilisation takes into account the distribution of each degree of freedom along the edge Again h 2 e can be used in this equation instead of e . The edge integral in (74) is evaluated numerically using Gauss quadrature where the Jacobian J ξ in the transformation d = J ξ dξ is the length of the k th edge J ξ = ∂ X ∂ξ = L k and ξ ∈ [0, 1] is the local coordinate.
When employing stabilisation (74) the resulting tangent matrix has full rank. Note, that stabilisation (73) is easier to implement than the second one.

Energy stabilisation
A stabilisation technique based on energy differences was first introduced for virtual solid elements in [38] for virtual elements. It is based on the idea of introducing the following stabilisation energy, see e.g. [23] The second term on the right hand side ensures consistency of the total potential energy since U stab (w h − (w h )) → 0 for h → 0 where h relates to the element size.
Often it is convenient to use the same strain energy ψ, see (17), for the consistency part U c ( (w h )) and the stabilisation part U stab (w h − (w h )). In that case a stabilisation parameter β is introduced to control the amount of stabilisation. By multiplying which can be reformulated as In this formulation it remains to compute U [ψ(w h )] which on a first glance seems to be impossible since w h is not known in the interior of a virtual plate element. The solution is to  [29] and the discrete Kirchhoff plate element (DKT), see [6]. In the following we will describe the main features of these elements.
Morley element It is based on six degrees of freedom and has six nodes with one degree of freedom at each node. The nodal values are related to the deflection w i at the vertices and the rotations θ k around the tangent of the element at the mid nodes. Thus the Morley element will match only part of the unknowns of the consistency part of element 1 and 2 and introduces extra degrees of freedom at the mid nodes in the interior of the virtual elements due to the submesh, see Fig.  4a.
The unknown nodal values can be combined in the vector p M T = {w 1 , w 2 , w 3 , θ 4 , θ 5 , θ 6 } T which yields where N I are quadratic shape functions that can be found explicitly in [37]. This ansatz leads to a non-conforming plate bending element with constant curvature and moments. With the ansatz (80) the plate curvatures can be computed from (15) leading tô where B M includes the derivatives of the ansatz functions N I . The strain energy of one element i k can then be defined following (17) with (16) as: Note that the integrand is constant. Hence the strain energy can simply be written as where i k is the area of one of the internal triangles. Since no integration is needed this yields a very efficient code for the stabilisation part U (w h ).
DKT element This element uses a shear deformable theory, but enforces the Kirchhoff constraint pointwise. It has three nodes with 3 degrees of freedom per node which are the deflection w and two rotations γ x and γ y around the X and Y axis, respectively. Hence the element has the same degrees of freedom as the virtual element 1 but misses the rotation at the mid side used in the virtual element 2. The unknowns are at the vertices, see the submesh in Fig. 4b. The DKT stabilisation does not introduce new unknowns compared to the consistency part. The interpolation leads to the vector of unknowns in the element p D T = {w 1 , γ x1 , γ y1 , w 2 , γ x2 , γ y2 , w 3 , γ x3 , γ y3 } T . Based on the introduced constraints special quadratic shape functions were derived in [6] which yield compatible rotations β x and β y in terms of the nodal unknowns p D T . The derivative of these rotations with respect to X and Y leads to the curvatures within the element which can be used to formulate the strain energy of one element following (17) and (16) Note that the integrand is of quadratic order and thus has to be numerically integrated using an adequate Gauss quadrature. Both elements used within the energy stabilisation need a local assembly of all triangular elements in the submesh of the virtual element.

Comparison of the different stabilisations and numerical examples
In this section we first study and compare the two classical stabilisations and then the energetic ones. This allow us to select one which be used subsequently all along the paper. The performance of the proposed virtual plate elements will be illustrated by means of numerical examples that are related to applications in engineering. Isotropic material behaviour as well as anisotropic materials and composites will be investigated. All numerical solutions are compared with analytical solutions where available.

Notation used in the examples
In general the following element formulations for plates will be analysed: The following types of stabilisation will be employed for the virtual plate elements: • st-1a/b: classical, nodal error stabilisation (73), where "a" refers to using the element size h 2 e and "b" to employing the element area e , both are used in the denominator. In general the following types of elements and meshes will be considered • Q1: quadrilateral mesh with 4 edges per element, • T1: triangle mesh with 3 edges per element, • VO-U: Voronoi-type mesh, with uniformly distributed cell seeds, • VO-R: Voronoi-type mesh, with randomly distributed cell seeds.
This means that Q1 meshes consist of quadrilateral elements which are virtual elements with four edges and 12 unknowns (element 1) or 16 unknowns (elements 2). In the same way a T1 mesh consists of virtual elements with triangular shape having 9 unknowns (element 1) or 12 unknowns (element 2). The actual number of nodes in a Voronoi type mesh depends on the shape of the Voronoi cells.
In our examples the number of edges in a Voronoi-type mesh varies from three to ten.

Study and comparison of the different stabilisations
The first example is used to compare results related to the different stabilisation techniques. Additionally the general convergence characteristics of element 1 and 2 are reported. A fully clamped square plate of size B × H = 8 m × 8 m is used for this purpose since an analytical solution is available. The plate is subjected to a uniform loadq see Fig. 5a. The material parameters are given in Table 2.
The analytical results of the clamped plate problem can be found in [34]. Table 3 provides the reference solutions for various quantities for the given data.
The study is based on the meshes depicted in Fig. 5b-d for a discretisation with square, elongated rectangular and Voronoi elements, respectively.

Study and comparison of classical stabilisation variants
The convergence study for virtual element 1 (VE-1) and virtual element 2 (VE-2) using the two types of classical stabilisation, is shown for square elements in Fig. 6, for rectangular elongated elements in Fig. 7 and for regular and unstructured Voronoi meshes in Fig. 8. The convergence study is made regarding the deflection at the center of the plate and the energy. One can observe that the choice of stabilisation does not have a big influence for the quadrilateral element VE-1. We will discuss the influence of the stabilisation when using triangular shape elements in Sect. 7. The effect between stabilisation type "a" and "b" related to the normalisation with element size h 2 e (st-a) or the area e (st-b) is negligible. But generally, using the area yields slightly better results in all cases. There is a slight difference in stabilisations 1 and 2. For Q1 meshes stabilisation 1 is the best, but for Voronoi meshes, stabilisation 2 is better by a small amount. It seems that the best stabilisation depends on the mesh type. As expected the structured meshes perform better than unstructured Voronoi type meshes. It is interesting to note that in this example the Q1 meshes yield the best results which will be explored further in Sect. 7. Figure 9 depicts contour plots of the deflection w and the bending moment M X X , for the virtual element VE-2:st-1b. but graphically all mesh types produce similar results and demonstrate that the developed virtual plate elements are capable of computing meaningful engineering solutions. The contour plots of the structured meshes as well as the unstructured Voronoi meshes report minimum and maximum bending moments (max M X X = 3.285 and min M X X = −1.466) that match exactly the analytical results in Table  3.

Study of energetic stabilisation
While studying the energy stabilisation using Morley and DKT element different issues where encountered. We will therefore first discuss the stabilisation of VE-1 using the Morley element before studying the stabilisation of VE-2 using the DKT element.
Energetic stabilisation of VE-1 using the Morley element Elements VE-1 and Morley are based on ansatz functions that lead to constant curvature within the element. Thus the Morley element was selected for stabilisation. One aspect of the energetic stabilisation is the parameter β which controls the influence of the stabilisation part of the energy. Figure  10 shows convergence curves obtained for the clamped plate under uniform load depending on β.
It appears that values of β ranging from 0.4 to 0.8 yield similar results. For smaller values of β the convergence behaviour deteriorates. In total it can be concluded that a stabilisation with the Morley element is not satisfactory since the best results is achieved with the classical stabilisation.
Energetic stabilisation of VE-2 using the DKT element Figure 11 depicts the convergence of the DKT element and the VE-2 element for a Q1 mesh. We first note that the convergence of DKT is lower than two. The convergence curves for a stabilisation with β = 0.2 and β = 0.05 increase the coarse mesh accuracy when compared to the DKT solution. A smaller value for β yields a better convergence. Up to a certain level of mesh refinement both stabilisations even lead to the second order rate of convergence, as is achieved with the classical stabilisation (st-1b). This in fact is not surprising since using a non vanishing value of β for small element sizes will prevent the second order convergence rate expected for VE-2 in the limit h e → 0. This observation led to a recipe in which β is scaled by the size of the element e in the refined mesh with respect to the element size 0 in the initial coarse mesh with 1 element leading to the formula β e = 1 2 e 0 . This scaling allows to recover the expected rate of convergence for VE-2 which is parallel to the curve using the classical stabilisation.
In conclusion the use of the proposed classical stabilisation is more simple and more efficient than the energetic one and therefore will be used in the rest of the paper. This conclu-  Deflection w (a-c) and Moments M X X (d-f) for structured and unstructured meshes, note that M Y Y is same but flipped over the domain diagonal Fig. 10 Convergence study related to the parameter β on a Q1 mesh for VE-1 sion could be different in the case of material non-linearities, a case for which it is not obvious to propose an efficient classical stabilisation. The results for stabilisation 1b and 2b are quite similar but stabilisation 1b seems generally to give the best results. In what follows we make use of stabilisation 1b systematically without referring to it anymore.

Dealing with singularities, the rhombic plate
Another example is the simply supported rhombic plate, see Fig. 12a, which is more complex due to a singularity. The plate is simply supported and loaded by a uniform load with H = B = 8 m. The angle α in Fig. 12a is chosen as α = 30 o which yields an internal, obtuse angle of 150 o . The material parameters were taken from Table 2.
Near the obtuse angles the solution is singular which leads to a solution belonging to H γ − , ∀ > 0 with γ = 2 − α π −α . Therefore the convergence is governed by the singularity for uniform meshes.
One of the advantages of the VEM is that it is easily possible to adaptively refine a mesh, see Fig. 12b, in a non uniform manner, see Fig. 12c allowing to recover a better convergence rate. The refinement is based on bisecting h e as shown in Fig. 12c. This can be simply done by adding a new nodes at the middle of the edge of the element since the number of nodes can be arbitrary in the virtual element scheme and avoids the burden of hanging nodes. Five refinement steps were executed towards the obtuse corner. Due to the local refinement one can observe the drastic increase of convergence rate compared to uniformly refining of the mesh by a factor of 2 as shown in Fig. 13. The used refinement however is not sufficient to recover the optimal convergence rate, see [5], but it shows the sensitivity of the solution with respect to the mesh refinement at the singularity. Here the adaptive refinement is performed in a way, that the difference of the energy between element nodes is a criteria for refinement. Thus to keep the difference small, the elements at the obtuse corners are refined, as depicted in Fig. 12c. In total 5 refinement steps were used to obtain the solution in Fig. 13. In this convergence study, additionally, the energy of the previous refinement step is used as an indicator for the refinement in the next step.
An analytical solution can be found in [28] for different plate angles, however it is only provided up to three digits. Using finite element convergence studies convergence results for the rhombic plate were reported in [5,36] with 3 digit accuracy, thus the obtained results related to the adaptive refined mesh do not converge linearly. We note, that the asymptotic convergence rate of 0.2 (related to the singularity) is recovered for all element formulations with a uniformly refined mesh which is clearly depicted in Fig. 13.
Finally, Fig. 14 shows the distribution of the bending moments M X X and M Y Y for the rhombic plate using a uniform mesh with 32 × 32 Q1 virtual elements.

Rectangular orthotropic plate
To investigate the convergence behaviour of the new plate elements for orthotropic material behaviour a rectangular plate is considered for which the analytical solution is provided in [35]. The orthotropic directions coincides with the main axis of the plate.
Using the same notation as in [35] the problem is characterised by a simplified material stiffness matrix, see (21), where t is the thickness of the plate whose length is denoted by a = 2 mm and its width by b = 1 mm. The plate is simply supported and subjected to a uniform load q. In the example the following geometric data are chosen as shown in Table  4. By noting that H = D 1 +2D xy , the deflection, see  We consider for the convergence analysis the deflection at the center (X , Y ) = ( a 2 , b 2 ) of the plate, denoted by w c . The value of w c is computed for the given data with Mathematica   Figure 16 depicts the error in the maximum deflection in the middle of the plate. We observe a quadratic convergence rate for element 2 for Q1 and Voronoi meshes VO-R and VO-U. As in the previous example, Q1 element performs best, however the convergence rate is more erratic. This example underlines that the Kirchhoff-Love virtual elements converge for these more intricate materials with the same order as for isotropic materials.

Plate with anisotropic material as cantilever beam
Next we deal with double cantilever plate specimens made of composite material. This setup was used in [2] to investigate the value of the critical energy release rate depending of the relative angle of a stack of plies adjacent to the considered interlaminar interface.  For the test setup shown in Fig. 17 bending solution for three stacking sequences are computed. The plate has dimensions B = 50 mm and width H = 20 mm. This plate is assumed to be clamped at the left side while a line loadq is applied on the right part. The material under consideration is a made of M18/M55J unidirectional plies. In what follows 1, refers to the direction of the fiber and 2 refers to the orthogonal direction inside the ply within the beam, see Table 5. We first compute the response for a solution with one ply whose orientation is at 22.5 o . In Fig. 18a the distribution of M X X is plotted with respect to an amplified deformed configuration. Since the orthotropic axis of the plate do not coincide with the global axis of the cantilever plate a coupled bending-torsion type of response is produced, as expected. This is the reason why those test are usually not performed. One prefers testing stacking sequences which avoid this kind of coupling as it is the case for stacking sequences of the type [φ, −φ, −φ, φ, −φ, φ, φ, −φ] s . Due to the symmetry of the stacking sequences decoupling between tension occurs with respect to the mid-plane but yields an orthotropic bending response. The total material matrix is computed from (21). For such stacking sequence and for φ = 22.5 o the distribution of M X X is plotted on an amplified deformed configuration, see Fig. 18b.

Use of the C 1 -continuous virtual element in a finite element environment
Since the developed virtual plate element can have arbitrary number of nodes it is possible to construct a triangular or quadrilateral C 1 -continuous plate element using element 1 or 2. Those elements can be implemented easily in conventional finite element codes. The advantage is that the construction of C 1 -continuous finite plate elements need a high ansatz order, like in the TUBA family of triangular elements in [4] or are-as quadrilaterals-a composition of four triangular elements, see [18]. The developed formulation provides the necessary C 1 -continuity for Kirchhoff plates with much simpler elements. It is interesting to investigate how the above derived virtual plate elements compare as triangular (3-noded) and quadrilateral (4-noded) elements with classical finite elements for Kirchhoff plates. Within two standard benchmark examples, see e.g. [30], we illustrate the performance of the derived virtual plate elements. The results are compared with finite elements that approximate that curvature in the same way as the derived virtual elements 1 and 2. For the virtual element 1 (VE-1) an equivalent Kirchhoff plate element with constant curvature is the nonconforming Morley element, already described in Sect. 5.3. For element 2 (VE-2) we select a non-conforming triangular plate element, developed by [33], and the DKT element, see Sect. 5.3, which rely on a linear curvature approximation.
Clamped plate under uniform load As first example we use a clamped plate under uniform loadq, already discussed in Sect. 6.2. The geometry and material data are the same as in in Sect. 6.2 and the analytical solution is given in Table 3. Results related to the convergence behaviour of a number of conforming and non-conforming plate elements are provided in [30] for the case of a uniform load.
We compared with the solutions obtained with finite elements from Morley and Specht and the DKT formulation which are plotted next to the results of the VE-1 and VE-2 elements used as standard triangles and quadrilaterals, see 19. For the pointwise stabilisation' 'st-1", see (73), Fig. 19a demonstrates that VE-1 has a better coarse mesh accuracy than the Morley element despite the fact that both elements have a constant approximation of the curvature. In the same way, VE-2 outperforms the DKT and Specht elements, although both elements approximate the curvature with a linear polynomial. It is interesting that the triangular element VE-1 has an extremely good coarse mesh accuracy when using the continuous stabilisation "st-2" in (74) as shown in Fig. 19b. The triangular element VE-1:T1 with stabilisation 2b is simple, efficient and has only three nodes. Thus it qualifies as a C 1 -continuous Kirchhoff plate element for legacy codes.
The rate of convergence is depicted in Fig. 20a for the constant curvature elements and demonstrates for all element formulations the same order of asymptotic convergence. Contrary to that, VE-2 has a higher convergence order when compared to the DKT and Specht elements as demonstrated  Fig. 20b. Thus the conforming virtual element 2 as a triangular well as as a quadrilateral element is performing extremely well. A distorted mesh, as depicted in Fig. 21 for quadrilateral and triangular element meshes with different densities, is used to investigate the behaviour of the quadrilateral and triangular virtual plate elements.
The convergence results are depicted in Fig. 22 for the continuous stabilisation "st-2". They show the same behaviour as the uniform mesh. Again the triangular element has better coarse grid characteristics than the quadrilateral element when using the continuous stabilisation, see Fig. 22a. Also the asymptotic convergence rate is maintained for the distorted mesh which is demonstrated in Fig. 22b.
Clamped plate under point load As a second example we apply the triangular and quadrilateral virtual element to the clamped plate under point load F = 64 kN. The geometrical and the material data are described in Sect. 6.2. The analytical solution is reported in [35] as w = −0.0896 F D . The plot in Fig. 23 depicts the convergence behaviour of the different element formulations.
As shown in Fig. 23a the pointwise stabilisation"st-1" yield the best results for VE-2 while Fig. 23b depicts the superior coarse mesh accuracy of element VE-1:T1 for the continuous stabilisation "st-2", as in the previous example. The asymptotic convergence behaviour can be observed in Fig. 24. Here again the coarse mesh accuracy of the VE-1:T1 element is demonstrated in Fig. 24a. Expectantly the same asymptotic convergence rate is achieved for VE-1: T1 and VE-1:Q1. Due to the point load the rate of convergence is lower for the higher order ansatz which can be seen in Fig.  24b. Here VE-2: Q1 depicts the best performance when using the pointwise stabilisation.
It is impressive to see in both examples, that the proposed virtual elements (VE-1:T1 st-2b and VE-2:Q2 st-1b) outperform the DKT element which, in the engineering literature, is known to be an excellent element. This qualifies both virtual plate elements as candidates for engineering software related to the analysis of thin plates. Especially the very simple triangular element VE-1:T1 st-2b is a good candidate since it fits easily in existing software packages, having only three nodes and the same number of unknowns at each node which is actually equivalent to using the Specht or DKT element.

Conclusion
The virtual element method is a simple and powerful tool to construct C 1 -continuous discretisations as needed in the case of Kirchhoff-Love theory for plates. This even works for elements with arbitrary shape and number of nodes. In this paper we have developed in a more engineering way of writing the theory presented in [13], extended its range of applications compared to [16] and also introduced and discussed several possible stabilisation techniques. Various examples and convergence studies were provided to demonstrate the ability of the virtual plate elements to solve engineering problems based on two simple element formulations. The second element is very appealing because despite its simple formulation it leads to second order accuracy and has, compared to existing plate elements, a very high coarse mesh accuracy. Moreover, both developed elements behave similarly well for isotropic and anisotropic cases. With its flexibility VEM is also well suited for non uniform mesh refinement which is required to recover the optimal rate of convergence in case of singularities encountered for example in rhombic plates.
Finally, it appears that restricting the element shapes to those of classical finite plate elements allows to enrich the set of elements available in legacy codes with simple and accurate C 1 -continuous plate elements that are simpler than the known C 1 -continuous finite plate elements and also much faster. This work could be extended to plate structures in the nonlinear regime, undergoing large deflections and rotations.
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/.

A Calculating the projection with AceGen
For the actual generation and implementation of the two elements and the associated stabilisations the software package AceGen, see [21,22] was employed. In the following the essential steps and matrices are summarized that were basis of the code derivation. The general way how to derivate the elements with the help of the automatic software generation tool AceGen is described in this appendix. Some steps of derivation of residual and tangent can be greatly simplified by using this software.
To obtain the unknown parameters a in (46) and (59), the following relations are considered where w = (w h ) and ∇w = ∇ (w h ) is used to simplify notation.