A Family of C0 Quadrilateral Plate Elements Based on the Refined Zigzag Theory for the Analysis of Thin and Thick Laminated Composite and Sandwich Plates

The present work focuses on the formulation and numerical assessment of a family of C0 quadrilateral plate elements based on the refined zigzag theory (RZT). Specifically, four quadrilateral plate elements are developed and numerically tested: The classical bi-linear 4-node element (RZT4), the serendipity 8-node element (RZT8), the virgin 8-node element (RZT8v), and the 4-node anisoparametric constrained element (RZT4c). To assess the relative merits and drawbacks, numerical tests on bending (maximum deflection and stresses) and free vibration analysis of laminated composite and sandwich plates under different boundary conditions and transverse load distributions are performed. Convergences studies with regular and distorted meshes, transverse shear-locking effect for thin and very thin plates are carried out. It is concluded that the bi-linear 4-node element (RZT4) has performances comparable to the other elements in the range of thin plates when reduced integration is adopted but presents extra zero strain energy modes. The serendipity 8-node element (RZT8), the virgin 8-node element (RZT8v), and the 4-node anisoparametric constrained element (RZT4c) show remarkable performance and predictive capabilities for various problems, and transverse shear-locking is greatly relieved, at least for aspect ratio equal to 5 × 102, without using any reduced integration scheme. Moreover, RZT4c has well-conditioned element stiffness matrix, contrary to RZT4 using reduced integration strategy, and has the same computational cost of the RZT4 element.


Introduction
In the last decades, multi-layered composite and sandwich structures have been increasingly used in different engineering fields (aerospace, automotive, marine, nuclear, and civil), thanks to their advantages over the traditional metallic structures due to their high specific stiffness and strength and their tailoring capabilities. Commonly built up with a sequence of layers of composite materials arranged thorough thickness, such structures typically exhibit a transverse shear deformability higher than their metallic counterpart, thus their structural analysis requires adequate models to take into account these effects.
In the framework of the displacement-based beam, plate, and shell theories, equivalent single layer (ESL) theories, where the displacement field is assumed smoothly continuous through-thethickness, and layer-wise (LW) theories, where a local displacement field is assumed for each layer, are commonly used; see the recent review of Abrate and Di Sciuva [1,2]. Among the ESL theories, for engineering applications, the classical plate theory (CPT) is widely used especially for thin plates. In thick laminates, the shear deformability is no longer negligible and must be taken into account. Contrary to the CPT, the first-order shear deformation theory (FSDT), first formulated by Mindlin and extended to the heterogeneous anisotropic plates by Whitney and Pagano [3], allows the rotation of the normal to the reference surface to be independent of the slope of the reference surface. In the FSDT, the thickness-wise distribution of the transverse shear stresses in a multi-layered plate is a piecewise constant function, thus due to the step-wise variation of the transverse shear moduli along the thickness, it violates the continuity requirement of the transverse stresses at the interfaces and the vanishing transverse shear stresses conditions on the top and bottom surfaces of the plate. It also needs an appropriate shear correction factor due to the assumption of linear transverse shear stresses distribution along the plate thickness. Generally speaking, these models are not capable to predict the stress distributions along the thickness with acceptable accuracy, the range of their applicability being limited to compute the global responses, like maximum deflection, natural frequencies, and buckling loads [4]. To improve the accuracy of the estimates of the in-plane displacements and stresses distribution along the thickness, the higher-order shear deformation theories (HOSDT) are used. In the framework of HOSDT, polynomials, trigonometric, exponential, and hyperbolic expansions of the axial displacement are commonly used; see Abrate and Di Sciuva [1,2]. These theories are generally more accurate than FSDT and do not require the use of ad-hoc shear correction factors but are unable to satisfy the transverse shear stresses continuity condition at the layer interfaces. Moreover, from a computational point of view, finite elements based on such theories require C 1 continuity shape functions for the transverse displacement [5][6][7][8][9], while FSDT elementbased require only C 0 continuity [4]. In LW theories, each layer has its own displacement fields, each independent from the others, which leads to a quasi-3D description of the displacement and stress fields. This assures a high degree of accuracy, but for laminates with several numbers of layers or complex structures, the computational cost of these theories can become prohibitive [10][11][12][13]. A compromise between the computational efficiency of ESL theories and the accuracy of LW theories is represented by the so-called zigzag theories (ZZT). In order to simulate the distortion of the normal and the step-wise distribution of the transverse shear stresses along the thickness typical of multilayered structures, in the ZZT the in-plane displacements are written as a superposition of a throughthe-thickness polynomial distribution (global contribution such as in the ESL theories) and a piecewise linear continuous function, (local or zigzag contribution). Unlike the LW theories, the kinematics of ZZT has a fixed number of kinematic unknown variables, independent on the number of material layers. Generally, the ZZ theories provide more accurate response predictions for relatively thick laminated composite and sandwich structures than the corresponding ESL theories, and generally comparable to those given by LW theories, at a lower computational cost. The first to investigate displacement-based zigzag theories was Di Sciuva [14][15][16]. From a computational point of view, the main drawback of the classical ZZT is that the finite element formulation requires C 1 shape functions, according to Di Sciuva [17][18][19][20]. In order to develop a zigzag theory free of the drawbacks of classical Di Sciuva's ZZ theory, Tessler et al. [21][22][23][24] developed the refined zigzag theory (RZT) for beams, plates, and shells. The reader interested in a comparison of the performances and range of applicability of various zigzag functions is invited to read the paper of Gherlone [25]. From a computational point of view, the formulation of finite elements based on the RZT requires only C 0 continuity of shape functions [26][27][28][29][30][31][32][33][34][35][36][37][38]. Although C 0 -continous elements are the most-widely used in commercial finite element software due to their computational efficiency, it is well known that overstiff unrealistic numerical results are obtained when using the conventional low-order isoparametric C 0 finite elements (bi-linear shape functions for both deflection and bending rotation) in the analysis of thin and very thin plates; this phenomenon is known in the literature as the transverse shear-locking effect, according to Oñate [38][39][40].
There are many approaches to overcome this difficulty. Generally, using higher-order isoparametric shape functions alleviates the shear-locking problem, but the computational cost increases. Therefore, higher-order elements are not well-suited for the nonlinear analysis (post-buckling analysis, transient nonlinear analysis, crash simulation, etc.) and failure analysis [41,42].
One of the more successful strategy of overcoming this problem in the 4-node bi-linear isoparametric plate element is to use the procedure commonly known as the reduced/selective integration technique to compute the transverse shear strain energy [38,[43][44][45][46].
Unfortunately, there is an undesirable effect associated with this scheme, i.e., there may be instances where extra zero strain energy modes will appear.
Other strategies investigated in the open literature are based on the penalty constraints [38]. Among these, the discrete penalty constraints, such as the discrete Kirchhoff theory (DKT) [47,48], where the transverse shear strain energy is neglected in the formulation and the Kirchhoff assumption of zero transverse shear strains are imposed at specific discrete points in the element are commonly used; the mixed interpolation of tensorial components (MITC) strategy [49,50], in which ad hoc chosen strain fields are tied to the displacement-based strains; B-bar method [51]; and F-bar method [52,53], which relies on the concept of deviatoric/volumetric split and the replacement of the compatible deformation gradient with an assumed modified counterpart.
Another strategy is the assumed natural strain (ANS) method [49,54], which consists of interpolating the strain field at a set of chosen points, known as tying points, where the strain values coming from the quadrature points are replaced, in a weighted manner, by assumed strains. A synergic use of the DKT and ASN strategies can be found in References [55,56]; the discrete shear gap (DSG) method [57], where the discretized shear gap field is determined by the interpolation of the nodal shear gaps. The anisoparametric interpolation strategy was proposed by Tessler et al. [58] for beam finite element and by Tessler [59] for plate elements. It requires that the shape functions for the transverse deflection be a complete polynomial one degree higher than the shape functions used to interpolate the bending rotations. This element was named anisoparametric virgin element [58]. To achieve simple nodal patterns, the constraint condition of constant effective transverse shear strain along the edges of the element is imposed. This allows one to eliminate the extra degree of freedom in the transverse deflection. The resulting element is named constrained anisoparametric, and it is simple, variationally consistent, and, more interestingly, computationally efficient.
In the present work, the attention is focused on the formulation, implementation, and numerical assessment of a family of C 0 quadrilateral plate elements based on RZT, with special emphasis on their numerical performances not only for thick, but also for thin and very thin plates. Using full integration Gaussian quadrature scheme the stiffness and mass matrices and nodal loads vector are computed. This work is an extension to the plate of the approach successfully applied to the beam in Gherlone et al. [26] and Di Sciuva et al. [34].
The first two elements of the family (RZT4 and RZT8) are isoparametric quadrilateral plate elements. RZT4 is a 4-node element with bilinear shape functions; RZT8 is an 8-node element with quadratic serendipity shape functions. Each node has seven degrees of freedom (dof's): Three displacements (u1, u2, w), two bending rotations (θ1, θ2), and two amplitudes of the zigzag functions, (ψ1, ψ2). The third one, named RZT8v (v stands for virgin), is an anisoparametric quadrilateral plate element, i.e., it is an RZT4 plate element for the six kinematic variables (u1, u2, θ1, θ2, ψ1, ψ2) and a RZT8 plate element for the transverse displacement, w. The last one, named RZT4c, is a 4-node anisoparametric constrained element. It is obtained from the RZT8v element by imposing the constraint condition of constant effective transverse shear strain along the edges of the element.
To numerically assess the relative merits of these four C 0 plate elements, in particular, of the RZT4c element, several cases are examined. Convergences studies and parametric investigations, under various boundary conditions and transverse pressure distribution, with regular and distorted meshes, are carried out by considering static and free vibration problems.
It is concluded that the bi-linear 4-node element (RZT4) has performances comparable to the other elements in the range of thin plates when reduced integration is adopted but presents extra zero strain energy modes. The serendipity 8-node element (RZT8), the virgin 8-node element (RZT8v), and the 4-node anisoparametric constrained element (RZT4c) show remarkable performance and predictive capabilities, at a reduced computational cost, for various problems, and the transverse shear-locking is greatly relieved without using any reduced integration schemes, for thin plate (aspect ratio equal to 5 × 10 2 ), at least for 8 × 8 mesh. Moreover, RZT4c has well-conditioned element stiffness matrix, contrary to RZT4 using reduced integration strategy, and the same computational cost.

Theoretical Basis of Refined Zigzag Theory (RZT) and Finite Element Modeling
We consider a rectangular multilayered flat plate of length a, width b, and total thickness h, assumed constant ( Figure 1). It is composed of a finite number N of orthotropic elastic layers of uniform thickness perfectly bonded together. The material of each layer is assumed to have a plane of elastic symmetry parallel to the reference surface (here chosen to be the midsurface of the plate) and principal orthotropy directions arbitrarily oriented with respect to the laminate reference frame (x1, x2).
The points of the plate are referred to an orthogonal Cartesian co-ordinate system is the set of in-plane co-ordinates on the reference plane, and ≡ is the co-ordinate normal to the reference plane ( Figure 1); the origin of the reference frame is fixed at the center of the middle-plane of the plate. So, ∈ (− , + ), ∈ (− , + ), and ∈ (− , + ).
If not otherwise stated, the superscript (k) is attached to quantities of the kth layer (k = 1, N), whereas the subscript (k) is attached to quantities evaluated at the kth interface (k = 1, N−1) between the k and (k + 1) layer. Also, we use the subscript b and t to indicate the top and bottom surfaces of to the coordinate i x ; a dot over a quantity refers to the partial derivative with respect to time of that quantity. Finally, the Einsteinian summation convention over repeated indices is adopted, with Latin indices ranging from 1 to 3, and Greek indices ranging from 1 to 2. The displacement field at time t in the RZT theory can be written as is the vector of the displacement components of the generic point belonging to the plate and ( , ) is the vector of the in-plane displacement components. As usual in the plate theory, the transverse displacement is assumed to be constant along the thickness, i.e., 3 In the refined zigzag theory, the in-plane kinematics is based on the superposition of a global (G) firs-order kinematics (that of the FSDT plate theory) and a local (L) layer-wise correction of the in-plane displacements (see, Figure 2). Thus, the in-plane displacement field at time t is assumed to be (see Figure 2) gives the contribution, which is continuous with its first derivatives with respect to the z-coordinate (the global (G), contribution of the first-order kinematics) and gives the local (L), contribution to the in-plane displacement, which is continuous with respect to z, but with jumps in the first derivative at the interfaces between adjacent layers. In Equations (4) and (5), where and are the displacements along the − and − axis of a point belonging to the middle plane of the plate; and are the bending rotation of the normal to the middle surface along the directions + and − , respectively; and represent the spatial amplitudes of the zigzag functions ( ) and ( ) , respectively. It should be noted that FSDT is a special case of the RZT, i.e., RZT reduces to FSDT when ( ) = (see Equations (1) and (2)). By taking into account these definitions, in Equation (1) For the th layer of thickness ℎ ( ) , the refined zigzag functions have the following expressions [34]: where ( ) and ( ) is the transformed reduced transverse shear stiffness modulus of the kth layer (see Equation (22)). Equations (9) and (10) show that the refined zigzag functions ( ) are a-priori known piecewise continuous functions of z, depend only on the thickness and transverse shear mechanical properties of the constituent layers, and vanish on the bottom = − and top = + surfaces of the plate (see Figure 2). In Equation (8) and in the body of the paper, I stands for the identity matrix and 0 for the null rectangular matrix, which dimensions follow from the rule of the matrix product and partitioning.

Virtual Work Principle and Discretized Equations of Motion
The discretized equations of motion are obtained via the dynamic version of the principle of virtual displacements (D'Alembert principle) where (δ is the variational operator) U The virtual variation of the work done by the internal stresses (strain energy) reads In Equation (12), σ  and ε  are the vectors of stress and linearized strain components, respectively. By taking into account the assumed displacement field, Equation (1), and assuming, as usual in plate theory, 33 0 is the vector of the stress and moment stress resultants for unit length, 23 12 , , The constitutive equations for a generally orthotropic layer are 1,1 In Equation (22), ( ) (i, j = 1, 2, 6) and ( ) (i, j = 4, 5) are the plane stress transformed stiffness moduli of the kth layer, that are functions of the z-coordinate. Substituting Equations (20) and (22) into Equations (16) and (17) yields the plate constitutive equations in terms of the generalized displacements The virtual variation of the work done by the inertia forces reads Substituting Equation (1) into Equation (26) yields In writing Equation (30), it is assumed that the plate is subjected to a transverse pressure z q applied on the top surface of the plate, and to boundary transverse loads per unit length applied on the edge parallel to -axis. Obviously, for free vibration.
In the previous equations, 3 ( ) x ρ is the material mass density and an overbar indicates the prescribed value of a quantity. The other symbols have been defined above. Due to difficulty to obtain closed form solutions, we search for an approximate solution.

Finite Element Formulation
As is well known, the accuracy of the finite element solution depends on the ability of the assumed shape functions to accurately model the deformation modes of the structure.
In the following sections, several C 0 quadrilateral plate elements derived using the displacement approach in conjunction with the RZT summarized in the previous Section will be described.

General Equation
Let us consider a general quadrilateral plate element with NN surrounding nodes. A generic 8node (NN = 8) quadrilateral plate element is shown in Figure 3. The reference plane is the physical plane ( , ). From a computational point of view, each element can be mapped in a conventional square element on the natural plane ( , ) where ( , ) ∈ (−1, 1), thus the numerical integration using the Gaussian quadrature can be easily implemented. As usual in the finite element method, let us expand the vector of unknown generalized displacements (Equation (7)) of points belonging to element e in the form where ( ) ( ) e t q is the vector of the nodal degrees-of-freedom of the element e and ( ) ( , ) e ξ η N is the matrix of the element shape functions in the natural co-ordinates ( , ) ξ η .
In Equations (26) and (39), J is the determinant of the Jacobian matrix of the co-ordinate transformation from the physical plane to the natural plane.

Shape Functions and Constrained Technique
From the variational statement used to derive the element governing equation, Equation (11), only the first derivatives of the kinematic variables appear in the virtual variation of the strain energy expression, Equation (13). Like FSDT, in formulating finite elements, refined zigzag theory (RZT) requires only C 0 shape functions, which makes it attractive in formulating beams, plate, and shell elements [38].
In this Section, the attention is focused on the formulation of the shape functions of a family of C 0 quadrilateral plate element based on RZT. Table 1 summarizes the topology and the shape functions of the developed family of C 0 quadrilateral plate elements. The first two (RZT4 and RZT8) are isoparametric quadrilateral plate elements. RZT4 is a 4-node element with bilinear shape functions; RZT8 is an 8-node element with quadratic serendipity shape functions. Each node has seven dof's: Three displacements (u1, u2, w), two bending rotations (θ1, θ2), and two amplitudes of the zigzag functions, (ψ1, ψ2) (see, Figure 1). For these elements, the same set of shape functions is used for approximating all the generalized displacements vector, As well known, a drawback of low-order C 0 finite elements (bi-linear shape functions for both deflection and bending rotation) coupled with the full Gaussian quadrature to compute exactly the integrals of the element matrices, is the shear-locking effect. To eliminate shear-locking, reduced integration strategy is often used. Unfortunately, there is an undesirable effect associated with this scheme, i.e., there may be instances where extra zero strain energy modes will appear.
Tessler and Dong [58], for the Timoshenko beam element, identified the reason of shear-locking effect in the Kirchhoff constraints of the shear strain measure for very thin structures. As a possible solution, they proposed to use the anisoparametric interpolation, where the deflection w is approximated by a polynomial with one order of degree higher than the others used for the remaining kinematic variable (bending rotation, θ). These elements, named Virgin elements by Tessler and Dong [58], have one extra node in the middle beam point having only the w-dof. This approach has been extensively used to formulate C 0 beam elements based on the RZT in Gherlone et al. [26] and Di Sciuva et al. [34]. When the anisotropic interpolation scheme is applied to quadrilateral plate elements, the shape functions are bilinear for all kinematic variables except for the transverse displacement, for which the 8-node serendipity shape functions are used, i.e., the extra node in the middle point of each edges have only the w-dof (see , Table 1 and Figure 4). This plate element is named RZT8v (v stands for virgin) in Table 1: It is an RZT4 plate element for the six kinematic variables (u1, u2, θ1, θ 2, ψ1, ψ2) and an RZT8 plate element for the transverse displacement, w. Therefore, the shape function matrices for the RZT8v element read (see Table 1 The vector of the nodal transverse displacements can be split in two subvectors Constrained Technique The last one, named RZT4c in Table 1, is a 4-node anisoparametric constrained element (RZT4c). It is obtained from the RZT8v element applying the constrained-edge-strategy, as suggested by Tessler and Dong [58] and Tessler and Hughes [60,62], i.e., by imposing the constraint condition of constant effective transverse shear strain along the edges of the element.
For the FSDT, it is possible to refer either to the transverse shear strains ( and ) or to the transverse shear stress resultants ( and ) using the following relations.   (63) and (64) 3, Figure 4. Eight-node virgin quadrilateral plate element. Geometry, physical co-ordinates (x1, x2), and natural co-ordinates (ξ, η). Boundary co-ordinate system (n, s). Definition of the parameters aij and bij entering in Equations (56) and (63) Considering the RZT, the previous relations (Equation (45)) have an additional term (see Equations (17), (24), and (25) It has been shown by Gherlone et al. [26] and by Di Sciuva et al. [34] that these two nondimensional transverse-shear-parameters have a range of values , ( 1,0) L T r r ∈ − , with , = 0 for homogeneous material properties through-the-thickness (when ( ) = 0 = 1,2), and , → −1 for highly heterogeneous cross-sections, such as laminated composites and sandwiches. As said, to generate a constrained element from the virgin type element defined before, following Tessler and Dong [58] and Tessler and Hughes [60,62], we enforce along the edges of the element the following condition , are the tangential edge rotation and the tangential zigzag amplitude, respectively; see Figure 5.
is the angle between the direction s and axis, The edge-coordinate is ∈ (− , ) where is the length of the edge defined by nodes i an j (see Figure 5). In the natural plane, the non-dimensional coordinate is or , in particular for the edges along the direction = and for the edges along the direction = .     )  12  12  12  12  , , ,  )  23  23  23  23  , , , where ( Figure 5) are the projections onto the axes of the edge's lengths (see Figure 5). Solving Equation (53) Finally, substituting Equation (58) into Equation (44) Therefore, the interpolation for the transverse deflection reduces to It can be shown that the first term of row matrix in Equation (60) and that the second term can be rewritten as follows  N a N a N a N a N a N  a N   b N b N b N b N b N N a N a N a N  a N   r  b N b N b N b N b N b N Thus, the RZT4c element has only the four corner nodes with the shape functions matrix The stiffness and mass matrices and the nodal force vector of the RZT4c, with respect to the RZT4 and RZT8, must be taken into account of the new correction introduced with the Equation (64). Substituting Equation (42) into Equation (34) and comparing the result with Equation (64), it is evident that RZT4c element has the same topology of RZT4 bilinear element (28 dof's), but has better behaviour in thin regime, due to the presence of the i S matrices in the third row, alleviating the shear locking problem.

Numerical Results
In order to assess the relative merits (convergence characteristics and numerical accuracy) of the developed family of C 0 quadrilateral plate elements when applied to the linear elastodynamic analysis of thick, thin, and very thin laminated composite and sandwich plates, several numerical tests have been carried out. The investigation was conducted on sample problems for which analytical solutions are available.
The static and free vibration behaviour of laminated composite and sandwich square plates simply supported on all edges is investigated, by comparing the global and local behaviour as estimated by Finite Element Method (FEM)with the analytical solutions. The results presented in this Section refer to centre non-dimensional deflection,  Table 2 gives the mechanical characteristics and material density of the unidirectional layers constituting the laminated composite and sandwich plates investigated in this study, whose acronyms (S stands for Sandwich, L for Laminate composite, and I for Isotropic), number and thickness of the layers, lamina orientation of the principal material direction 1, and material of the single layer, are given in Table 3. In order to compare results from analytical computations and those coming from the finite elements, some results are given using the following ratios:

RZT Analytical Performances
The accuracy and reliability of the RZT plate theory in predicting the static, dynamic, and buckling behaviour of thick laminated composite and sandwich beams, plates, and shells have been assessed elsewhere [21][22][23][24]35,36]. Thus, no special emphasis will be given here to analytical results.
We consider a symmetric thick (span-to-thickness ratio, a/h = 6) sandwich square plate with thin two-layered cross-ply face-sheets (S) simply supported on all edges, under bi-sinusoidal transverse pressure of intensity 0 q . As an example of the accuracy of RZT model, Figure 6 shows the thicknesswise distribution of the non-dimensional in-plane displacement , maximum transverse displacement, , bending stresses, and , and transverse shear stress, . For the purpose of comparison, the exact elasticity solution obtained using Pagano approach [68] is also plotted. Note that transverse shear stress distribution along the thickness is obtained by integration of the local 3D elasticity equilibrium equations in the absence of body forces. Notice the great accuracy of the RZT. It is also important to remark that in RZT no shear correction factor is used, contrary to the FSDT.

Convergence Analysis
In this Section, the static convergence characteristic of each element of the family (RZT4, RZT8, RZT8v, and RZT4c) is investigated.
Due to geometry, boundary conditions, and loading symmetry with respect to the reference axes, if not otherwise specified, all FEM results have been obtained meshing only one-quarter of the plate with a regular mesh, i.e., all the elements in the physical plane are square elements (see Subsection 4.5). Thus, the number of elements indicated below refers to the number of elements along the side of one-quarter of the full-plate.
Thick (a/h = 10) and very thin (a/h = 10 3 ) square plates, simply supported on all edges, under uniform pressure of intensity 0 q , are considered for this analysis: A sandwich plate with thin twolayered cross-ply face-sheet (S), a five-layered cross-ply plate (L), and an isotropic plate (I). The reference (analytical) solutions (Table 4) have been obtained using the Ritz method [69] for the deflection. For the implementation and the convergence analysis of this approximate method, see Reference [69]. In particular, as it is shown in Reference [69], 10 orthogonal polynomials in both x1 and x2 directions are enough to assure the converged Ritz solution. It is clear that when the numerical analyses are performed using the RZT (analytical or with FEM) for the isotropic plate (I), the zigzag terms vanish. Thus, the finite elements formulated on the RZT, for this particular case, coincide with the FSDT ones when the classical shear correction factor is equal to 1 (k 2 = 1) and these results are quoted with FSDT4, FSDT8, FSDT8v, and FSDT4c*.
In the RZT4c element, two non-dimensional transverse-shear-parameters ( , ) appear. For this element, three cases are considered: In the first one, and are computed as defined by Equation (47) and results are quoted with the label RZT4c ; in the second one, only for the sandwich (S) and laminated (L) plates, we assume the material highly heterogeneous with = = −1 and results are quoted with the label RZT4c ; in the third one, we assume = = 0, such as the isotropic case (I), also for the sandwich (S) and laminate (L) plates and results are quoted with the label RZT4c * . All the FEM results have been obtained using full integration strategy. Examination of the results plotted in Figures 7-9 reveals that the convergence for the most of all elements is monotonic in character, with the exception of the constrained element when used in very thin regime (a/h = 10 3 ) and coarse meshes (2 × 2 or 4 × 4). Generally, the convergence rates are different.    Among the four elements, RZT8 is the element with the best convergence rate, followed by RZT8v and RZT4c elements (with different values for non-dimensional transverse-shear-parameters) that are very close to each other and to RZT8 element. On the other hand, RZT8 element has a large number of nodes when compared to the other elements (see Table 5), so its use in the nonlinear static and transient dynamic and failure analyses can result in a very high computational cost, also if meshes with few elements are used. RZT8v has bilinear shape functions for all kinematic variables, with the exception of the deflection, for which quadratic serendipity shape functions are used. Therefore, the number of degrees of freedom is less than the RZT8 element, but higher than RZT4 and RZT4c elements (see Table 5). RZT4c has the same number of nodes as RZT4, but its convergence characteristics are comparable to that RZT8v element. Only the classical four-node RZT4 with bilinear shape functions has a low convergence rate for maximum deflection. Considering the isotropic case (I) the FSDT8v and FSDT4c elements are equivalent in thick and very thin regimes. Except for the RZT4 element, for the 16x16 mesh, the error in the FEM results is not higher than 3%. The convergence rates of RZT4c are very similar to the RZT8v, without a great influence of the non-dimensional transverse-shear-parameters ( , ). In particular, this element guarantees the same accurate results of the RZT8v with the same dof's of the RZT4 element. Due to this, the RZT4c element is very interesting from a computational point of view. In Table 5 are reported the total dof's for each mesh discretization using the various elements: It appears that the RZT4c gives the best compromise between accuracy and low computational cost.
Thus, in the following analyses, if not otherwise specified, an 8 × 8 regular mesh for a quarter of plate is used.

Effect of Non-Dimensional Transverse Shear Parameters
In the derivation of the RZT4c element from the RZT8v one, there is an element-material dependency by the presence in the element formulation of the two non-dimensional transverse shear parameters and . In the previous figures (Figures 7-9), the effect of these parameters is shown in the convergence results. For the isotropic case, these parameters are no relevant because there is not the zigzag effect. For the cross-ply (L) and sandwich (S) plates, this effect become important in the FEM results. When these parameters are enforced to be = = 0, such as the isotropic case, also for the non-homogeneous case, the convergence to the analytical solution is lower than the other cases. This condition is not the best to characterize this element. Now we consider the RZT4c and RZT4c cases. Table 6 shows the values for and computed for the cross-ply (L) and the sandwich (S). Table 6. and values for cross-ply (L) and sandwich (S).
Laminate L -0.1841 -0.1841 The approximation, used for RZT4c cases, is good if we consider the sandwich case, but some differences are present for the cross-ply plate. However, if we consider the results for the cross-ply, there are very few differences between considering and exactly computed or use = = −1 . Furthermore, it is preferable not having a material dependency in the element transverse shear strains. This could be a problem if we consider two or more adjacent elements where the material properties are different from one element to another, because at that point there will be two different transverse shear strains.
By taking into consideration these aspects and the good accuracy of the results for the RZT4c case, in the remaining part of the paper, the RZT4c element is used only considering the case = = −1.

Shear Locking Phenomenon
In order to assess the behaviour of the four plate elements with respect to the shear-locking effect, Figure 10 shows results for the maximum deflection (rw) ratio as a function of the span-tothickness ratio, a/h, for a simply supported five-layered symmetric cross-ply plate (L) under bisinusoidal transverse pressure, considering a 4 × 4 and an 8 × 8 mesh. The results quoted as RZT4 refer to the 4-node quadrilateral element using full integration. As is well known, the shear-locking effect is pathological for the full integrated RZT4 elementas is well known, the use of stabilisation techniques could avoid the appearance of hourglass modes; other code-specific techniques that have been also proposed (typically based on EAS or ANS), do not introduce any spurious mode, for both 4 × 4 and 8 × 8 mesh patterns. RZT4 is suitable only for the analysis of thick and moderately thick plates. RZT8 shows very good performances over a very large range of span-to-thickness ratios (from 6 to 10 3 ), with some initial effects of transverse shear locking appearing for very thin plate (a/h = 10 3 ). A good compromise between the accuracy of the RZT8 and the computational efficiency of the classical four node appears to be RZT8v. RZT4c perform very similar to RZT8v, but with a lower computational cost, equal to that of RZT4 element.
As we said, in the literature and in the FEM commercial codes, the shear locking effect is typically addressed using the reduced integration technique when computing the shear terms in the stiffness matrix. Figure 11 shows the behaviour of the RTZ4 element using full and reduced integration. It is easily concluded that the reduced integration strategy improves the performances of RTZ4. Yet, one typical problem affecting the reduced integration technique is the singularity of the stiffness matrix   Table 7 gives the number of rigid body modes for the various plate elements. As expected, we have three rigid translations along the axes and three rigid rotation around the axes for all the elements, with the exception of the RZT4 with reduced integration, where an extra zero-energy mode (spurious mode) appears. The other elements, RZT8, RZT8v, and RZT4c, are free from this unlikely effect. In the remaining part of the paper, the full integration is used for all the C 0 family elements.

Distorted and Regular Mesh Analysis
As is well known, the performances of the plate finite elements depend also on their skewness ratio, i.e., on the degree of the distortion of their plan form geometry. Figure 12a,b respectively gives the regular and distorted 16 × 16 mesh for the whole plate used in this study.
To assess the relative performances of the RZT elements when used in distorted geometry, a cantilever square sandwich plate (S) and multi-layered cross-ply plate (L) are investigated. Thick (a/h = 8) and very thin (a/h = 10 3 ) plates are considered. Table 8 gives the percent errors of FEM results on the maximum non-dimensional deflection under three types of load distribution (in Figure 13 are shown: Uniform distribution, constant uniform edge load, and two vertical opposite tip forces), using regular and distorted 16 × 16 meshes. For simplicity, the load cases are named as follows: C1 for the uniform transverse load pressure, C2 for the constant uniform edge load, and C3 for the two vertical opposite tip forces. For the same plates and meshes, Table 9 gives the percent errors of FEM results on the undamped fundamental frequency parameter. The percent error is defined as The analytical solutions have been obtained using the Ritz method [69]. In particular, as it is shown in Reference [69], 10 orthogonal polynomials in both x1 and x2 directions are enough to assure the converged Ritz solution.   As expected, regular meshes give results closer to the analytical ones. As before, RZT8 has the best performances also for distorted meshes, but at expense of higher computational costs. Again, RTZ4c appears to be a good compromise between accuracy and computational cost.
To further assess the performance of the RZT elements, the bending behaviour of a sandwich (S) cantilevered plates is investigated, comparing the bending moment and the transverse shear force distributions along the midline of the plate with the analytical ones. Two square plates (S) are considered, both loaded with the load distribution C1 (uniform transverse load) and C2 (uniform constant edge load) and meshed with 17 × 17 elements: (i) A thick plate (a/h = 8; a/b = 1) and (ii) a thin plate (a/h = 10 2 ; a/b = 1). The analytical solutions have been obtained using the Ritz method [69] with 10 orthogonal functions in x1 direction and 10 orthogonal functions in x2 direction. Plots of the distributions of shear force and bending moment along the midline of the plates parallel to x1 axis are given in Figure 14 (thick sandwich plate) and Figure 15 (thin sandwich plate) for the C1 load condition and in Figure 16 (thick sandwich plate) and Figure 17 (thin sandwich plate) for the C2 load condition. The shear force and bending moment were computed using the finite element correspond to the centroid of each elements on the midline of the plate. The shear force and bending moment are normalized as follows for the C1 load condition: Also, in this test, RZT4c element shows its good accuracy in estimating the shear force also near the clamped edge for thick sandwich plates.

Stress Distributions
In order to show the ability of FEM solutions to capture also the local behaviour (distribution of the stresses along the thickness), we consider bending under transverse uniform pressure of simply supported sandwich plate (S) for a/h = 8. Figure 18 compares the distributions of in-plane stresses at the centre of the plates and the transverse shear stress at the midpoint of the left edge, computed using the Ritz method [69] and the family of C 0 RZT plate elements. For this study, in order to have a better distribution of the stresses, a regular mesh 16 × 16 for a quarter of plate has been adopted. These results confirm in general the good accuracy of these elements to estimate the thicknesswise distributions of the in-plane stresses, while some difference when the transverse shear stresses are considered. As expected, for both sandwich (S), the RZT8 shows the best performances, followed by the RZT8v and the RZT4c. By taking into account that RTZ8 is computationally more expensive and that RTZ4 with reduced integration is penalized by the zero energy strain modes, it can be concluded that RZT4c represents the best compromise between accuracy and computational costs.

Concluding Remarks
A family of four C 0 quadrilateral plate elements based on the refined zigzag theory has been developed and numerically tested.
The numerical assessment given in the previous Section reveals that: 1. The convergence of the FEM solution for all the elements of the family is in general monotonic in character, although with different rates. 2. The bi-linear 4-node element (RZT4) is strongly affected by the shear locking effect in thin plates when full integration is used. 3. The bi-linear 4-node element (RZT4) has performances comparable to the other elements in the range of thin plates when reduced integration (RI) is adopted, but presents extra zero strain energy modes, if no stabilization technique is used. 4. The serendipity 8-node element (RZT8), the virgin 8-node element (RZT8v), and the 4-node anisoparametric constrained element (RZT4c) show good performance and predictive capabilities for various problems, and the transverse shear-locking is greatly relieved without using any reduced integration schemes, for thin plate (aspect ratio equal to 5 × 10 2 ) at least for mesh 8 × 8. 5. All elements, except for the conventional bi-linear RZT4 element with full integration, are adequate to predict the bending global response (transverse displacement, distribution of the transverse shear resultant, and bending moment) and undamped frequencies. The same holds for the predictive capabilities of the local behaviour (thickness-wise distributions of the bending stresses and transverse shear stresses). 6. RZT4c has well-conditioned element stiffness matrix, contrary to RZT4 when reduced integration strategy is used. 7. The computational cost of RTZ4c is much lower than that of RTZ8 and RTZ8v elements. 8. The accuracy of all the elements is sensitive to their degree of planform distortion.
Therefore, due its good numerical accuracy (generally comparable to that of RZT8 and RZT8v elements) and low computational cost (equal to that of RZT4 element), RTZ4c quadrilateral plate element appears to be a good candidate for the FEM analysis of laminated composite and sandwich plates, mainly when nonlinear or transient linear and nonlinear analyses are involved. As a general concluding remark, the presented results show clearly that RZT8 is the most robust element among those considered in the present investigation. Moreover, as it is well known, the quadrilateral isoparametric 8-node FE is not an optimal element. The final conclusion is that further studies based on different approaches should be carried out to propose a locking-free RZT plate FE. Therefore, following Tessler and Hughes [60,62] suggestions, in a work in progress, the concept of an appropriate element shear correction factor is adopted in conjunction with the anisotropic constrained strategy to further alleviate the shear locking effect (extend the range of the plate aspect ratios until 10 6 ).