Free and Forced Vibration of Laminated and Sandwich Plates by Zig-Zag Theories Differently Accounting for Transverse Shear and Normal Deformability

A number of mixed and displacement-based zig-zag theories are derived from the zig-zag adaptive theory (ZZA). As a consequence of their different assumptions on displacement, strain, and stress fields, and layerwise functions, these theories account for the transverse shear and normal deformability in different ways, but their unknowns are independent of the number of layers. Some have features that are reminiscent of ones that have been published in the literature for the sake of comparison. Benchmarks with different length-to-thickness ratios, lay-ups, material properties, and simply supported or clamped edges are studied with the intended aim of contributing toward better understanding the influence of transverse anisotropy on free vibration and the response of blast-loaded, multilayered, and sandwich plates, as well as enhancing the existing database. The results show that only theories whose layerwise contributions identically satisfy interfacial stress constrains and whose displacement fields are redefined for each layer provide results that are in agreement with elasticity solutions and three-dimensional (3D) finite element analysis (FEA) (mixed solid elements with displacements and out-of-plane stresses as nodal degrees of freedom (d.o.f.)) with a low expansion order of polynomials in the in-plane and out-of-plane directions. The choice of their layerwise functions is shown to be immaterial, while theories with fixed kinematics are shown to be strongly case-sensitive and often inadequate (even for slender components).


Introduction
Fibre-reinforced laminated and sandwich composites are fundamental to obtain a faster speed, longer range, larger payloads, less engine power, and a better operating economy of land, sea, and aerospace vehicles.Indeed, they possess excellent specific strength and stiffness properties along with many other advantages.However, their behaviour may be strongly influenced by local effects in their multi-phase structure, which can cause a relevant loss of strength and stiffness, and even lead to a possible catastrophic failure in service.To avert this, a three-dimensional description of their displacement and stress fields is required in the computer simulations.Displacements are no longer C 1 -continuous as in homogeneous non-layered materials, but instead they should be C 0 -continuous (zig-zag effect), this being the only way to guarantee the continuity of transverse shear and normal stresses and of the transverse normal stress gradient that is necessary to satisfy the fulfillment of local equilibrium equations in layered materials.
A multitude of theories is currently available for laminated and sandwich composites, which account for layerwise effects at varying degrees of accuracy and computational costs.The papers by Carrera et al. [1,2] and by Demasi [3] are cited wherein a broad discussion of this matter can be found.Theories can be summarily categorized into equivalent single-layer (ESL), discrete-layer (DL), and zig-zag (ZZ) formulations, and further into displacement-based and mixed theories, because displacements, strains, and stress fields can be chosen separately from one another.
Given the limited accuracy offered by ESL even to predict overall response quantities for certain loading, material properties, and stack-up [4,5], they still remain widespread by virtue of their low computational cost (see Burlayenko et al. [6] and by Jun et al. [7]), but nowadays, DL and ZZ are becoming more and more widespread in their applications.Indeed, they have the merit of accurately predicting the displacement and stress fields irrespective of lay-up, layer properties, and loading or boundary conditions.However, DL could overwhelm the computational capacity when structures of industrial interest are analysed (unless their use is limited to critical areas) due to having too many variables.Instead, ZZ having intermediate characteristics between ESL and DL can strike the right balance between accuracy and cost-saving, allowing designers' demand of theories to be met in a simple, already accurate form.In particular, ZZ accounting for the transverse normal deformability effect has been proven to be suitable for carrying out ply-level stress analyses at a lower cost than DL.
Relevant examples of dynamic studies carried out by DL and ZZ that highlight their superior performance include the papers by [1][2][3]6,7] Boscolo and Banerjee [8], Khdeir and Aldraihem [9], Sayyad and Ghugal [10], Kazanci [11], Lin and Zhang [12], Vescovini et al. [13], and in addition to the previous ones, [ [14][15][16][17][18][19] are also cited.Moreover, papers [1][2][3][4][5]8,13,20] are examples that prove the limitations of ESL and the importance of the transverse normal deformability under static localised loading for certain boundary conditions and lay-ups [21] and even for accurately predicting the first free-vibration mode frequencies.Although ZZ theories are ultimately finding an ever-increasing number of applications thanks to their accuracy at an affordable cost, the literature shows that undeservedly Di Sciuva's zig-zag-like fuction DZZ is less used than Murakami's zig-zag-like fuction MZZ (for a definition of acronyms, see Table 1), despite their better accuracy with a low order of expansion of analytical solutions, as also demonstrated by the numerical results of this paper.
To contribute to the dissemination of DZZs and their further development, hereafter the typologies and specific characteristics of available ZZs are briefly reviewed.The intended aim is to explain the different behaviour of DZZ and MZZ as a consequence of their different layerwise contributions, namely Di Sciuva's [22] (DZZ) or Murakami's [23] (MZZ) zig-zag functions.
DZZ incorporates layerwise contributions as the product of assumed zig-zag functions and zig-zag amplitudes, the latter being computed to a priori satisfy the interfacial compatibility of out-of-plane stresses (for this reason, they are referred in literature as physically-based ZZ).Instead, MZZ assumes zig-zag functions that a priori feature a periodic change of the slope of displacements at interfaces, irrespective of the stack-up.As a consequence, stresses are assumed apart from the kinematics within the framework of Hellinger-Reissner variational theorem.For this reason, they are called kinematic-based ZZ.So, the development of MZZ is easier, and an efficient C 0 formulation can more easily be generated, but they could not appropriately represent the displacement fields when the orientation angle of layers is aperiodic, as shown in the literature.According to findings by [1,2,5,[24][25][26], it can be stated in general that MZZ are accurate only with a rather high expansion order of variables across the thickness and in the in-plane directions.As shown in [27][28][29][30][31][32][33], instead, DZZ accurately predict displacements and (very often but not in all cases) stresses with a lower expansion order.
As shown in the papers by Carrera et al. [1,2,[24][25][26], Demasi [3], Mattei and Bardella [27], Icardi et al. [28,29], Li and Liu [30], Zhen and Wanji [31,32], and Shariyat [33], the transverse normal deformability, which is usually neglected because it is erroneously considered of little importance, on the contrary should be accounted for.Indeed, his contribution becomes of primary importance when boundary conditions other than simply-supported edges are considered, under a strong variation of mechanical properties across the thickness and for localised loading.However, a rather complex formulation of theories, i.e., with a sufficient number of parameters that can be defined in order to satisfy all of the relevant physical constraints, is required for these cases.Although very often power and Taylor's series expansions are used across the thickness, hierarchic polynomials, and trigonometric and exponential functions, a combination of both or radial basis functions could be required to more efficiently represent variables [34][35][36][37][38][39].
The idea of separately hiring kinematics and stresses, so as to account for the transverse normal deformation without having to represent in a piecewise way the transverse displacement, was applied by Barut et al. to DZZ [40].Other mixed DZZ that are based on a similar idea have a rather simple kinematics, but yet, the capability to quite accurately predict stresses has been formerly developed in HR form by Kim and Cho [41] and Tessler et al. [42] Despite the notable contributions mentioned up to now, further research is required on ZZ because recent studies by Zhen and Wanji [43], Gherlone [44], and Groh and Weaver [45] have shown that MZZ could be less accurate than RZT [42], assuming the same degree of representation across the thickness, while [43] shows just the opposite for different cases.So, it is necessary to ascertain even whether three-dimensional (3D) DZZ can be as accurate as MZZ with a minor computational burden, so as to be effectively suitable for industrial purposes, considering also other theories in addition to those already examined in [44,45], especially those with a piecewise representation of the transverse displacement.
A massive amount of free vibration studies is currently found in the literature for cross-ply laminates, but other boundary conditions and lay-up that are equally interesting for practical applications need to be considered when testing the accuracy of theories (see, e.g., Li et al. [46]).It should be noted that for clamped edges, an incorrect vanishing transverse shear force resultant is predicted by ESL plate models; therefore, it is necessary to check on an adequate number of cases as to whether DZZ are immune from such mistakes.Also, lay-ups with a quite large variation of properties of layers, and in particular soft-core sandwiches, should be investigated to enhance the existing database.Another subject that needs further studies is what happens when zig-zag functions that are different from those commonly used [22,23] are assumed.It is also necessary to better clarify what effects have an a priori assumption of certain zig-zag amplitudes at certain interfaces, or even to all of the interfaces such as in MZZ, in order to save costs.
To contribute to this matter, in this paper, zig-zag theories in displacement-based and mixed form are particularized from the ZZA 3D zig-zag theory [29], assuming different layerwise functions.Consider that ZZA, here schematically defined as the adaptive zig-zag theory, assumes a variable kinematics with fixed degrees of freedom (d.o.f.) that consist of in-plane and transverse displacements u α 0 , w 0 and transverse shear rotations Γ 0 α at a middle reference plane.Moreover, it should be noticed that ZZA is capable of great accuracy, even with strong anisotropy at a cost that is still comparable to that of ESL.
The purpose of this paper is to test the accuracy of theories derived from ZZA under simplified hypotheses on displacement, deformations, and stress fields, so as to further reduce the computational cost, in order to assess whether and when the accuracy of ZZA can be preserved.The accuracy and the efficiency of these theories, which for the most part allow a redefinition of the coefficients layer-by-layer so as to satisfy all of the physical constraints and account for the normal transverse deformability effect (although in different ways), are assessed considering the free vibration and the response behaviour of blast-loaded multilayered and sandwich plates.
A new theory referred to as ZZA* is developed from ZZA in order to prove that once the coefficients of displacements are redefined across the thickness, the same accuracy degree can be achieved irrespective of the layerwise functions that are chosen.It will also be shown that as the coefficient is redefined, the same results are obtained even when zig-zag contributions do not explicitly appear.Other theories have been previously developed by the authors [47] in displacement-based and mixed form, either with Di Sciuva's or Murakami's zig-zag-like functions, and which also have characteristics similar to other theories that are already known in the literature, are considered for comparison.
Numerical applications aim to show the superiority of ZZA* over ZZA, including its lower computational costs with the same accuracy, for the quite large range of variation of lay-ups and boundary conditions that are considered.Benchmarks that are either retaken from the literature or new are considered, the former in order to enable a comparison of available theories, and the latter in order to enhance the existing database.The results by the present theories are compared to exact solutions whenever available or to 3D finite element analysis (FEA) [48].The findings can be categorized as the new results of known benchmarks provided by the new theories and entirely new results for the new benchmarks, which could serve as test beds for future analytical and finite element models.
This study also aims to show the superiority of ZZA* over all of the other theories considered for most of the cases, since the use of simplified assumptions implies a loss of accuracy.It will be shown that the accuracy of the simplified theories, in particular those of Murakami, is strongly case-sensitive; therefore, they cannot be used interchangeably.Studies are carried out in closed form considering the free and forced vibrations of laminated and sandwich beams and plates, simply-supported and clamped edges, different lay-ups, and distinctly different material properties of constituent layers, so as to give rise to relevant 3D effects.

Preliminaries and Notations
Constituent layers are assumed to have a uniform arbitrary thickness h k and linear elastic properties.They are assumed to be perfectly bonded to each other and, as usual when the global-scale response is examined, the existing bonding resin interlayer is disregarded.For the same reason, sandwiches are described in homogenised form as multilayered structures with one or more thick soft intermediate layers as the cores, with the cell-scale effects being disregarded.
A rectangular, right-handed Cartesian coordinate reference system (x, y, z) is assumed as the reference frame, having (x, y) on the middle reference plane Ω of the laminated plate (origin in the lower left edge) and z as the thickness coordinate (z ∈ [−0.5h, 0.5h], with h being the overall thickness).Lx and Ly symbolise the plate side-length in the x and y-directions, while (k) z + and (k) z − represent the upper and lower positions of the layer interfaces, respectively.Subscripts k and superscripts k are used to indicate that a quantity belongs to the layer k, while u and l mark the upper and lower faces of the laminate, and a comma is used to indicate spatial derivatives, e.g., (.) ,x = ∂/∂x, (.) ,z = ∂/∂z.The elastic in-plane and transverse displacement components are indicated as u α and u ς , respectively.Strains are assumed to be infinitesimal and, to distinguish their origin, they are specified as σ kl , respectively when they come from kinematic [.] u or stress-strain [.] σ relations.Once assumed as primary variables, they are indicated as ε ij , ε ας , ε ςς (i, j ≡ 1, 2; α ≡ x, β ≡ y; 3 ≡ ς ≡ z; γ αβ = 2ε αβ ).Stresses from stress-strain constitutive relations are indicated as ), while when they are assumed as primary variables, they are indicated as σ ij , σ ας , σ ςς .

Recalls on Mixed Variational Theorems
Here a generalized version of Hu-Washizu theorem is used, whose primary displacement boundary condition link is weakened as S u (u i − u i )n j ∂σ ij dS = 0 and which generates the following variational statement: Here HW g represents the canonical functional, n j are the components of the external unit normal to the volume bounding surface, b i are the components of body forces (they will contain inertial forces, as explained in (4)) and dv ensure the consistency of assumed strain and stress fields with their counterparts obtained from stress-strain and strain-displacement relations.As usual, the prismatic volume V of the laminated plates is assumed to be bounded by a surface S that is split into a surface S t on which surface tractions are prescribed and a surface S u on which surface displacements are prescribed.Body forces b i on V, prescribed surface tractions t i on S t and prescribed displacements u i on S u are assumed to act.
Theories with displacements and out of plane stresses assumed separately are developed from the HR variational theorem: where δΠ HR is the first variation of the HR canonical functional, and the following assumptions are made: , being i, j = 1, 2 ≡ x, y; 3 ≡ z.Dynamic governing equations are obtained from the previous variational statements (1), (2), or from the principle of virtual work (PVW), for theories in displacement-based form, accounting for the work of inertial forces.

Construction of Analytical Solutions
Closed-form solutions to dynamic governing equations are obtained, irrespective of the theory examined, expressing functional d.o.f. as a truncated series expansion of unknown amplitudes A i ∆ and trial functions i (x, y) that individually satisfy the prescribed boundary conditions: Then, we substitute these expressions within PVW, HR, or HW functionals and operating as specified immediately below.Here, ∆ symbolises in turns u 0 , v 0 , w 0 , θ x , θ y , because middle plane displacements and rotations of the normal are assumed as the only functional d.o.f. for each theory of this paper.Mechanical boundary conditions are accounted for by determining a number of unknown amplitudes A i ∆ , in proportion to the number of boundary conditions enforced, using the Lagrange multipliers method to account for the relationships resulting from each mechanical condition.The remaining amplitudes are determined deriving the governing functional with respect to still unknown amplitudes and equating to zero, having considered the work of inertial forces: within the functionals.In this way, an algebraic system is obtained whose solution provides the numerical value of each amplitude; then, the displacement, strain, and stress fields can be computed.
Since symbolic calculus being used to construct the theories, the applied distributed loading can be conveniently defined as a continuous or discontinuous general function χ(x, y) acting on upper and/or lower faces (or just on a part of them), so that energy contributions can be constructed in exact closed form.In this way, a series expansion representation with a large number of components is not necessary to represent discontinuous or otherwise complex loading distributions; then, the construction of the structural model in the numerical applications can be simplified, and at the same time made more accurate.
The trial functions have been adopted in individual applications are explicitly defined below and in Table 2 along with the expansion order used in each case.Note that although it is different for each benchmark, the order of expansion used is the same for all theories when a specific one is examined being equal to the minimum one that makes ZZA accurate.This is done to homogenise results and in order to compare theories under the same conditions.
At the clamped edge of cantilever beams, hereafter assumed at x = 0 by way of example, the following conditions are enforced the following further conditions are enforced: to simulate that (5) holds identically across the thickness.To ensure that the transverse shear stress resultant force equals the constraint force, the following additional constraint is enforced, while at the free edge x = L such resultant force is enforced to vanish It should be noted that the boundary conditions in Equations ( 6) to (8) are enforced using Lagrange multipliers methods; no further conditions are enforced on bending moments, but if necessary, they could be coerced choosing the sufficient expansion order in Equation (3).
At the supported edge of the propped cantilever beams at x = L, the following support condition: is enforced at the lower face z = −h/2, while the condition in Equation ( 8) is reformulated as: At the simply supported edges, the following boundary conditions are enforced w 0 (0, y) = 0; w 0 (L x , y) = 0; w 0 (0, y) ,xx = 0; w 0 (L x , y) ,xx = 0 w 0 (x, 0) = 0; w 0 (x, L y ) = 0; w 0 (x, 0) ,yy = 0; w 0 (x, L y ) ,yy = 0 (11) on the reference mid-plane of plates at x = 0, x = L x , and y = 0, y = L y .Appropriate changes corresponding to the ones for simply supported beams are obtained.Table 2 provides the trial functions and expansion order that is assumed for each benchmark.  +) A uniform mesh is used; x a and y b represent the number of elements in x and y directions, respectively, z h is the number of elements across the thickness.
Any other boundary condition could be enforced in the same way, namely by choosing the trial functions that individually satisfy the prescribed boundary conditions.In the event of this not being possible and in order to satisfy the mechanical conditions, Lagrange multipliers method should be applied.In the numerical applications for cases g and h, the boundary condition to be further enforced is u ς (x, h/2) = −u ς (x, −h/2), because some of their first free vibration modes exhibit anti-symmetric displacements in the transverse direction.In this way, a closed-form solution can be obtained with a reduced computational effort, instead of resorting to FEA, as is often done for such cases.

Higher-Order Theories
The laminated plate theories that are used as structural models are thoroughly examined below.Here, they are grouped into higher-order and lower-order ones.The latter represent simplified versions of the former, having mostly features that are similar to those of the theories that have been previously proposed in the literature.In the following section, it will be specified which are new and which have been previously developed by the authors.Governing equations will not be reported in explicit form, as they can be obtained in a straightforward way using standard techniques.

Features of the ZZA Theory and the Higher Order Theories Derived from It
As it forms the basis of all theories considered in this paper, the theoretical framework of ZZA below is expounded first.The through-thickness displacement field is postulated as [29]: Three kinds of contributions that are distinctly separated into lower-[. ..] 0 , higher-[. ..] i and layerwise [. ..] c are incorporated, whose specific purpose is described below.

•
[. ..] 0 is a linear contribution and contains the five functional degrees of freedom of the theory.• [. ..] i contains higher-order terms.Any combination of independent functions could be assumed to represent [F u α (z)] i and [F ζ (z)] i ; however to include theory [28] as a particularization of ZZA [29], the following form is chosen: Higher-order contributions [(Oz 4 . ..)] i , are characteristic of ZZA, while the terms [ 3 (. ) α ] i , [ 4 (. ) ζ ] i are the same as in the previous theory [28].The closed-form expressions of coefficients C i α , D i α , b i to e i are obtained using symbolic calculus from enforcing the fulfilment of stress boundary conditions Here, p 0 (±) represents the distributed transverse loading acting on the upper p 0 (+) and lower p 0 (−) faces.Of course, also, non-homogeneous conditions σ ας ; σ βς = 0 could be enforced without any additional difficulty.For clarity, contributions [. ..] i from Equations ( 14) and (15) are rearranged in the following way: The lower-order contributions under the square brackets in Equation ( 17) are determined by enforcing the fulfilment of the boundary conditions in Equation ( 16) and the local equilibrium equations: at selected points across the thickness.Higher-order contributions A α4 z 4 + . . .+ A αn z n , A ς5 z 5 + . . .+ A ςn z n , which enable a variable-kinematics representation across the thickness, are computed for each fictitious computational layer i in which the laminate is subdivided by imposing the fulfilment of Equation (18).However, except otherwise stated, a third/fourth-order representation that embraces the whole laminate is used in the applications which is adequate to obtain accurate results.Indeed, as shown in [21], this choice is a valuable combination of accuracy and cost-saving also when extreme variations of material properties give rise to very strong layerwise effects.
A single computational layer is used for laminates, while two or three layers are used for sandwiches.Note that the in-plane position of equilibrium points can be chosen appropriately for each case.• [. ..] c represents the layerwise contributions; the expressions of zig-zag amplitudes Φ k α , Ψ k and Ω k are determined so that the continuity of out-of-plane stresses and the transverse normal stress gradient σ ςς,ς at the layer interfaces is satisfied, as prescribed by the elasticity theory through the following stress compatibility conditions: at the physical and mathematical layer interfaces.Layerwise contributions α c j u and C j ς restore the continuity of displacements at the mathematical layer interfaces.The symbols n i and n in the summation of Equations ( 12) and ( 13) are used to distinguish the number of physical interfaces from that of the mathematical layer interfaces, respectively.In detail, Φ k α enables the continuity of transverse shear stresses, while Ψ k , Ω k enable the continuity of the transverse normal stress and of its gradient at physical and mathematical layer interfaces.All together, these provide the right slope changes of displacements at the interfaces of layers with different material properties and/or orientations.
Elsewhere in this paper, the symbols − and + indicate the position just before and just after the interface, respectively.The term (z − z k )H k appearing in in Equations ( 12) and ( 13) is Di Sciuva's zig-zag function [22], while (z − z k ) 2 H k is Icardi's parabolic zig-zag function [28], with H k being the Heaviside unit step function (H k = 0 for z < z k , while H k = 1 for z ≥ z k ).
The enforcement of the equilibrium and stress compatibility conditions in Equations ( 18) and ( 19) yields to a system of algebraic equations at each interface that is solved in closed form once and for all using a symbolic calculus tool.Notice that if just the material properties and/or the orientation of layers change, but not their number, symbolic expressions representing the solution will remain the same.
It will be shown forward that (z − z k )H k and (z − z k ) 2 H k can be replaced with any other layerwise function keeping the results unchanged, provided that coefficients inside [. ..] i are re-computed for each computational layer, as outlined above.Layerwise functions can even be omitted if a sufficient number of coefficients are incorporated in [. ..] i whose expressions are determined by enforcing the fulfilment of the interfacial stress compatibility conditions in Equation (19).To show this, a new theory ZZA* is developed in Section 3.1.1,which is devoid of [. ..] c , but incorporates new unknown coefficients within [. ..] i , which will be indicated as [. ..] i+c .Since this choice speeds up the computations of the coefficients of each layer, it turns into a computational advantage that grows with the number of computational layers, so it is worth taking into consideration.
The expressions of α C j u and C j w are determined in a straightforward way by enforcing the continuity of displacements at the mathematical layer interfaces Notice that as no d.o.f.derivatives are involved, the computation of α C j u and C j w is much easier and faster than those of Φ k α , Ψ k , and Ω k , which instead involve such derivatives.However, it must be considered that the computation of all of the above indicated zig-zag terms takes only an infinitesimal fraction of the overall calculation cost, and so remains compatible with that of ESL.
The strain energy updating technique (SEUPT) [29] can be used to obtain a C 0 formulation of the ZZA theory, since derivatives of the d.o.f. are involved in the displacement field as a consequence of the enforcement of Equations ( 18) and ( 19).

ZZA* Displacement-Based Theory
This theory, which represents the new theoretical contribution brought by this paper, is developed with the intended aim to demonstrate that the choice of zig-zag functions is immaterial, provided that the coefficients of displacements [. ..] i are recomputed as indicated above.
For this purpose the displacement field of ZZA* is assumed to be similar to that of ZZA except for the layerwise functions.In numerical applications it will be shown that the same accuracy of ZZA can be achieved without explicitly incorporating contributions [. ..] c , and consequently obtaining a reduction of the computational burden.
It will be shown by the numerical results that the choice of zig-zag functions is immaterial if the coefficients [. ..] i are redefined layer-by-layer, so these functions do not even need to be explicitly incorporated in the displacement field.The same result wouldn't be achieved by keeping the coefficients fixed across the thickness, in which case the accuracy depends on the choice of zig-zag functions.
The displacement field of ZZA* is conceived in the following way: wherein the terms k B i α and k C i α serve the same purpose as Φ k α and α C j u in Equation ( 12) inside ZZA, while k b i and k c i have the same function of Ω k and Ψ k , and k d i has the function of C k ς in Equation ( 13).In the same way of ZZA, C i α , D i α , b i , c i , d i and e i allow the stress boundary conditions in Equation ( 16) and the local equilibrium equations in Equation ( 18) to be met.More specifically, b i and c i enable the fulfillment of the stress boundary conditions concerning σ ςς and σ ςς,ς over the lower bounding face, while they are cancelled in subsequent layers.Instead, C i α , D i α , d i , and e i allow satisfying the three equilibria (18) at two points for each intermediate layer.In the lower layer, such coefficients enable the two boundary conditions on σ ας to be enforced along with the three equilibrium equations at a single point.In this way, free variables still remain that allow meeting three equilibrium equations at a single point across the upper layer and the boundary conditions at the upper bounding surface.When more equilibrium points are desired, each layer can be subdivided into two or more computational layers.However, this doesn't mean an increased expansion order or number of variables, since coefficients can be recomputed using the same order of representation for each computational layer, and the d.o.f.remain fixed.As for ZZA, closed-form expressions of all of the coefficients of ZZA* within contributions [. ..] i+c are determined once and for all for a specific lay-up using a symbolic calculus tool.

HWZZ Mixed Theory
Such a theory was developed in [47] in order to reduce the computational effort of ZZA by keeping only the essential contributions to the displacement, strain, and stress fields, and hopefully preserving its accuracy.Developed within the framework of the Hu-Washizu theorem, HWZZ forms the basis of another mixed theory here used and called HWZZM, which considers a different zig-zag layerwise function.

Master Displacements
The displacements of HWZZ are derived from those of ZZA neglecting the contributions of Ω k , because they are supposed to give imperceptible slope variations compared to those by Φ k α , Ψ k .Also higher-order and adaptive contributions A α4 z 4 + . . .+ A αn z n and A ς5 z 5 + . . .+ A ςn z n are neglected and no decomposition into mathematical layers is allowed.As a consequence, contributions by α C j u and C j ς are omitted.So, the displacement field of HWZZ is written:

Master Strains
Out-of-plane strains ε zz , γ xz , γ yz are constructed on the basis of the ZZA layer-by-layer representation of displacements as: The symbol (.) states that they refer to the computational layer .On the contrary, no decomposition into mathematical layers is allowed for the in-plane strains: The expressions of membrane stresses σ xx , σ yy , and σ xy are obtained in a straightforward way from stress-strain relations and previous strains, while those of the out-of-plane counterparts are assumed as specified next.

Master Stresses
Out-of-plane stresses are obtained from membrane stresses by integrating local equilibrium equations: and then In this way, stress jumps resulting from the omission of contributions by Ω k are recovered.As a consequence of simplifying the assumptions that are made, HWZZ needs to be post-processed by ZZA in order to accurately represent displacements when very strong layerwise effects rise.

HWZZM Theory
Murakami's zig-zag contributions are assumed [51], which are variations of the canonical form where: In fact, the displacement field of HWZZM is assumed as: Differently to Murakami's theories proposed in the literature, here, multiplier coefficients A u ε k , A u ς k , and B u ς k are incorporated, which can be defined differently across the thickness to improve the accuracy.As has been done previously, their expressions are a priori determined by enforcing the continuity of the transverse shear and normal stress, and of the transverse normal stress gradient at layer interfaces, so that they are no longer uniform across the thickness.As before, C k α and C k ς restore the continuity of displacements at interfaces of mathematical layers.
HWZZM is developed starting from the displacement field in Equations ( 29) and ( 30) in the same way as HWZZ; namely, the decomposition into fictitious computational layers is not allowed for the displacement and in-plane strain fields, while it is allowed for out-of-plane strains.Contributions to in-plane displacements over the third-order and to the transverse displacement over the fourth-order are neglected.Stress-boundary conditions are enforced at the top and bottom laminate faces, while local equilibrium equations are enforced at the inner layers.Membrane stresses σ xx , σ yy , and σ xy come from stress-strain relations, while out-of-plane master stresses are derived integrating local equilibrium equations.The numerical findings will show that although different zig-zag functions contradistinguish ZZA, HWZZ, and HWZZM, their results are indistinguishable, so there will be evidence that the choice of such functions is immaterial, provided that [. ..] i and [. ..] c are recomputed at each interface.

HWZZM* New Theory and Theories Derived from It
In order to assess the effect of the choice of different zig-zag functions, a new theory called HWZZM* is derived from ZZA*, similar to how HWZZ was derived from ZZA.Then, other theories are particularized from HWZZM, assuming A u α k (z), A u ς k (z), and B u ς k (z) are equal to those at a specific interface, which differ from one another.

HWZZM* Theory
The displacement field of HWZZM* is assumed in the following form (32) whose amplitudes are computed at each interface from the enforcement of the stress compatibility conditions in Equation (19).In this case, terms [. ..] i only serve to satisfy the stress-boundary conditions and local equilibrium equations, while the satisfaction of stress compatibility conditions is demanded to [. ..] c .
It should be noted that Equations ( 31) and ( 32) imply a little reduction of the processing time, which is equal to 10% per each layer with respect to ZZA and to a somewhat lesser extent equal to 6% with respect to HWZZ.Of course, such an advantage will become more consistent as soon as the number of computational/constituent layers increase.

Theories with a Priori Chosen Zig-Zag Amplitudes
Several theories are particularized assuming a priori zig-zag amplitudes as the ones by HWZZM competing at a specific interface and then keeping them unchanged across the thickness, which is something similar to what is done using Murakami's zig-zag functions, whose amplitudes are assumed a priori.
The amplitudes of HWZZMA are assumed to be coincident with those of HWZZM at the first interface from below, while the A u ς k (z) and B u ς k (z) of HWZZMB are assumed to be the same as those of HWZZMA at the first interface from below, but instead, A u α k (z) is calculated as in HWZZM theory.In HWZZMC, only B u ς k (z) is assumed to be uniform across the thickness and coincident with that of HWZZM at the first interface from below, while the remaining are computed at each interface.B u ς k (z) are neglected in HWZZM0 theory, while A u α k (z) and A u ς k (z) are assumed in the same way as HWZZMB.
The HWZZMB2 and HWZZMC2 theories are the same as those of HWZZMB and HWZZMC, respectively, but currently, the amplitudes are assumed to be coincident with those of HWZZM at the first interface from above.Since zig-zag amplitudes are assumed, discontinuous out-of-plane stress may result; hence, the integration of local equilibrium equations is required to obviate this discontinuity, with a corresponding increase in costs by 0.9%.However, it must be considered that zig-zag amplitudes are not being computed at each interface, so a processing time saving of 10% is obtained, and in the end, a positive balance is achieved.However, applications will show a loss of accuracy, making it vain.

Lower-Order Theories
Lower-order theories are particularized from ZZA through limiting assumptions that take features that are reminiscent to those of the theories in the literature.

MHR Theory
MHR considers piecewise cubic in-plane displacements, wherein Murakami's zig-zag function in Equation ( 28) is incorporated as the layerwise function, alongside a fourth-order polynomial transverse displacement [47]: Coefficients C α , D α , a, b, c and d are still calculated by enforcing the fulfilment of the stress boundaries conditions in Equation ( 16), while u αz is calculated by enforcing the fulfilment of the first and second equilibrium in Equation ( 18) at the middle plane of the laminate.Since out-of-plane stresses may be still discontinuous at layer interfaces, and since u αz is assumed to be uniform across the thickness, the expressions of out-of-plane stresses are derived integrating local equilibrium equations within the framework of an HR variational theorem (2).A refined version that is referred as MHR±, is obtained by determining the right sign of Murakami's zig-zag function (28) on a physical basis, instead of being forced to reverse at interfaces by the coefficient (−1) k .The right slope is determined at each interface without a cost burden evaluating what sign (±1) k attain the lowest residual force norm from three local equilibrium equations.

MHR4, MHWZZA and MHWZZA4 Theories
A refined variant MHR4 of MHR is obtained assuming the in-plane displacement field u α by Equation ( 31), and a fourth-order piecewise variation of the transverse displacement [47]: Coefficients a to d are determined by enforcing the fulfillment of the stress boundary conditions in Equation ( 16), whereas w z is calculated by enforcing the fulfilment of the third local equilibrium equation at the middle plane.
MHWZZA is developed assuming the same master displacement field as that of MHR in Equation ( 33) and the same master strain and stress fields as those of the HWZZ model, in Equations ( 25) and (26), respectively.To improve the accuracy, the displacement, strain, and stress fields of MHWZZA are recovered using ZZA as the post-processor.
Similarly to MHR±, a refined theory MHR4± is obtained from MHR4 determining the sign of Murakami's zig-zag functions on a physical basis.
A further theory MHWZZA4 is derived assuming the in-plane displacement field in Equation ( 31) by MHR, the transverse displacement in Equation ( 13) by ZZA, and as the master strain and stress fields those by HWZZ in Equations ( 25) and (26), respectively.The only substantial difference of MHWZZA4 with respect to HWZZ and ZZA is a different zig-zag function; the just-mentioned theories along with ZZA* and HWZZM* make it possible to verify whether the accuracy is sensitive to the choice of zig-zag functions.

HRZZ and HRZZ4
HRZZ theory is developed postulating a uniform transverse displacement and a third-order zig-zag representation of in-plane displacements [47]: Within the framework of the HR theorem, the transverse normal stress σ 33 is assumed to be the same as that of the ZZA model, while the transverse shear stresses σ i3 are derived from the equilibrium equations assuming kinematic relations in Equation (25) to define membrane stresses.However, currently, second and higher-order derivatives of the d.o.f. are neglected, and a unique computational layer is assumed.Since a uniform transverse displacement is chosen and transformed, reduced stiffness properties are assumed.Then, (33) , and εσ 33 = 1/C 3333 (σ 33 − C 33ij ε u ij ).In order to increase the accuracy of HRZZ, the ZZA theory will be used as the post-processor and the results obtained in this way will be indicated as HRZZ PP in the figures.
HRZZ4 assumes the same in-plane representation of HRZZ, the following fourth-order polynomial approximation of the transverse displacement [47]: and the same stress fields of HRZZ.In this case, ε u 33 being no longer null, it is unnecessary to use transformed, reduced stiffness properties.Similar to in the previous theories, coefficients b to e of Equation ( 36) are determined by enforcing the stress boundary conditions at the upper + and lower − faces (16).The out-of-plane master shear stresses σ i3 are obtained through integrating local equilibrium equations, while the ∂σ 33 appearing in Equation ( 2) is assumed to be the same as that of ZZA.

Numerical Assessments and Discussion
The accuracy of previous theories is assessed, analysing free vibration modes and the transient dynamic behaviour under impulsive loading of simply supported and clamped, laminated, and sandwich beams and plates.
The lay-up, geometric, and material properties of each case are grouped in Tables 3 and 4; the trial functions and the expansion order that are used in each study are reported in Table 2; the normalization of quantities is brought in Table 5; while the results are carried in Tables 6-20 and the computational cost is reported in Tables 21 and 22.
Table 5. Normalization of displacements, stresses and frequencies.

Case Normalization
w , w static response

Propped Cantilever Sandwich Plate in Cylindrical Bending under Uniform Static Loading
Firstly, a static analysis of two propped cantilever sandwich plates in cylindrical bending under a uniform load p 0 (+) on the top face are considered, in order to preliminary assess the accuracy of the theories in predicting the displacement and stress fields.The structure is clamped at x = 0 and restrained on the lower face at x = L x ; the elastic moduli of faces and core are assumed to be E u /E l = 1.6 and E u /E c = 166.6,respectively, and with a Poisson's ratio υ = 0.3; the lower (l) and upper (u) faces are t l = c/2 and t u = c/4 thick, respectively, where c is the core thickness.The results of the theories are compared to the ones by 3D FEA [48], whose elements are formulated in order to fulfil Equations ( 16) and (18), and the displacements and out-of-plane stresses are assumed as nodal d.o.f.The length-to-thickness ratios L x /h = 5.714 (case a1) and 20 (case a2) are considered.
Regarding case a1, according to Mattei and Bardella [27] and Icardi and Sola [21], its throughthickness displacement and stress fields are very challenging, and require a very accurate description of the transverse shear force resultant at the clamped (7) and at the simply supported edges (10), and a very accurate description of the transverse normal stress; otherwise, incorrect stress predictions are obtained.Moreover, what makes this a tough case is the opposite sign that is assumed at the upper and lower faces by the transverse shear stress at the supported edge.Figure 1 shows the results for the thickest case (L x /h = 5.714).It can be seen that, contrary to what was postulated by Murakami's zig-zag function, at the supported edge, the slope of u ς never reverses, while u α doesn't reverse near the upper interface, so MHR and MHR4 obtain inaccurate results.Even MHR± and MHR4±, whose interfacial displacement slope is determined on a physical basis, are not adequate, because their kinematics are too poor.MHWZZA4 and MHWZZA, incorporating strain and stress fields from HWZZ, obtain better results than MHR and MHR4, while HRZZ calculates an incorrect null transverse displacement at the supported edge.Summarizing, all of the lower-order theories obtain inaccurate results.Since the results of HWZZM are in good agreement with those of the adaptive theories ZZA, ZZA*, HWZZ, and 3D FEA, it is proven that the choice of layerwise functions is immaterial as long as higher-order coefficients are defined as in Sections 3.1 and 3.3.HWZZMB, HWZZMB2, and HWZZM0 assuming arbitrarily zig-zag amplitudes can accurately predict the axial displacement but not the transverse one, while HWZZMC and HWZZMC2 calculate precise displacements, whereas HWZZMA always provides wrong results (so they are not reported in the figures).It should be noticed that the behaviour of theories arbitrarily assuming zig-zag amplitudes is very case-dependent, and only some are quite accurate.
The results for the slender case L x /h = 20 are reported in Figure 2; the calculations show that a still significant difference between the predictions of theories exists, and σ ας differs from that which is expected for the rather thin sandwich structures.Again, only higher-order theories obtain a good degree of accuracy, and HWZZMA is so inaccurate that results cannot be reported.The obvious conclusion is that even for a length-to-thickness ratio of 20, an accurate description of kinematics is of primary importance in the present case.
On the basis of the processing time required (see Tables 21 and 22 and Section 5.10), it can be concluded that the adaptive theories ZZA, ZZA*, HWZZ, HWZZM, and HWZZM* are the most efficient ones, as they combine accuracy and low cost.It should be remembered that the same representation with the same low expansion order (Section 2.2) dictated by industrial needs is used for all of the theories.Since the results in the literature for MZZ have shown that greatly increasing the expansion order along with that of the representation across the thickness obviously produces more accurate results, prior statements are valid only for the conditions that are examined here.Below, dynamic tests are carried out that consider progressively more challenging benchmarks that highlight the need for sophisticated theories.
Aerospace 2018, 5, x FOR PEER REVIEW 18 of 43 calculates an incorrect null transverse displacement at the supported edge.Summarizing, all of the lower-order theories obtain inaccurate results.Since the results of HWZZM are in good agreement with those of the adaptive theories ZZA, ZZA*, HWZZ, and 3D FEA, it is proven that the choice of layerwise functions is immaterial as long as higher-order coefficients are defined as in Sections 3.1 and 3.3.HWZZMB, HWZZMB2, and HWZZM0 assuming arbitrarily zig-zag amplitudes can accurately predict the axial displacement but not the transverse one, while HWZZMC and HWZZMC2 calculate precise displacements, whereas HWZZMA always provides wrong results (so they are not reported in the figures).It should be noticed that the behaviour of theories arbitrarily assuming zig-zag amplitudes is very case-dependent, and only some are quite accurate.
The results for the slender case Lx/h = 20 are reported in Figure 2; the calculations show that a still significant difference between the predictions of theories exists, and   differs from that which is expected for the rather thin sandwich structures.Again, only higher-order theories obtain a good degree of accuracy, and HWZZMA is so inaccurate that results cannot be reported.The obvious conclusion is that even for a length-to-thickness ratio of 20, an accurate description of kinematics is of primary importance in the present case.On the basis of the processing time required (see Tables 21 and 22 and Section 5.10), it can be concluded that the adaptive theories ZZA, ZZA*, HWZZ, HWZZM, and HWZZM* are the most efficient ones, as they combine accuracy and low cost.It should be remembered that the same representation with the same low expansion order (Section 2.2) dictated by industrial needs is used for all of the theories.Since the results in the literature for MZZ have shown that greatly increasing the expansion order along with that of the representation across the thickness obviously produces more accurate results, prior statements are valid only for the conditions that are examined here.Below, dynamic tests are carried out that consider progressively more challenging benchmarks that highlight the need for sophisticated theories.

Free Vibration Modes of Simply Supported Laminated and Sandwich Plates in Cylindrical Bending
First, the analysis of simply supported cross-ply [0/90/0] (case b1) and [0/90]2 (case b2) plates in cylindrical bending is presented, which primarily serves as a preliminary test of the accuracy of 3D FEA [48] in solving dynamics problems, this having yet been tested.Then, a [0/core/0] simply supported sandwich plate in cylindrical bending (cases b3) is studied, which has more marked layerwise effects.

Free Vibration Modes of Simply Supported Laminated and Sandwich Plates in Cylindrical Bending
First, the analysis of simply supported cross-ply [0/90/0] (case b1) and [0/90] 2 (case b2) plates in cylindrical bending is presented, which primarily serves as a preliminary test of the accuracy of 3D FEA [48] in solving dynamics problems, this having yet been tested.Then, a [0/core/0] simply supported sandwich plate in cylindrical bending (cases b3) is studied, which has more marked layerwise effects.
The first fundamental frequency predicted by the present theories for case b1 is given in Table 6 considering a length-to-thickness ratio of 10.Comparisons are carried out with the exact solution and the results of EFSDT and EHSDT theories by Kim [49].
It can be seen that only ZZA, HWZZ, HWZZM, ZZA*, and HWZZM* higher-order theories provide very accurate results, but a sufficient accuracy is obtained by the HRZZ, HRZZ4, MHWZZA, MHWZZA4, MHR, MHR4, MHR±, MHR4±, and HWZZMC simplified theories.Instead, HWZZMA, HWZZMB, HWZZMB2, and HWZZM0 give an inaccurate prediction of it.However, case b1 does not appear to be severe, as many theories prove adequate.Although they have dissimilar characteristics, it is noted that HSDT and FSDT (shear correction factor of 5/6 chosen to minimize the error) calculate even the fundamental frequency inaccurately.The same considerations apply for case b2, with the same length-to-thickness ratio of 10, as shown in Table 6; therefore, the previous considerations are not repeated.The fundamental frequency for a [0/core/0] simply supported sandwich plate constituting case b3 [49] is reported in Table 7 for three different length-to-thickness ratios (four, 10 and 20).Contrary to what one would expect, errors don't dramatically decrease as Lx/h increases, and the behaviour of theories remains quite diversified.In particular, MHWZZA, MHWZZA4, and HWZZMA, HSDT, and FSDT (shear correction factor is 5/6) prove to be inaccurate irrespective of the length-tothickness ratio.Modal displacements and stresses are reported In Figures 3-5, where the results are normalized as shown in Table 5, and show that the axial displacement u α and the in-plane stress σ αα are accurately calculated except by MHWZZA and MHWZZA4, while the transverse shear stress σ ας is inaccurately calculated by MHR, MHR4, MHR±, MHR4±, MHWZZA, and MHWZZA4 for the thickest case.An even bigger scattering is shown for u ς , which is accurately calculated only by ZZA, HWZZ, HWZZM, ZZA*, and HWZZM*, while on the contrary, the transverse normal stress σ ςς is accurately calculated by all of the theories.Note that the results by HWZZMA are never reported in Figures 3 and 4 as being too wrong.For the same reason, the stress and displacement distributions are not reported for FSDT and HSDT.Such big errors highlight the inapplicability of the most simplified theories even for slender cases, which is contrary to what is often claimed in the literature.
 is accurately calculated by all of the theories.Note that the results by HWZZMA are never reported in Figures 3 and 4 as being too wrong.For the same reason, the stress and displacement distributions are not reported for FSDT and HSDT.Such big errors highlight the inapplicability of the most simplified theories even for slender cases, which is contrary to what is often claimed in the literature.

Fundamental Frequency and Modal Displacements and Stresses of a Simply Supported, Anti-Symmetric Cross-Ply Plate
The fundamental frequency of [0/90/0/90], simply supported, cross-ply square plates with a length-to-thickness ratio of five [4], is examined in Table 8 for increasing the values of the degree of orthotropy E1/E2 ranging from three to 40 (case c1).
Comparisons are given with the results by a global-local higher-order theory [4] and by Matsunaga [56], which were obtained using a higher-order ESL.It can be seen that all of the theories of this paper provide fairly accurate predictions of the fundamental frequency in correspondence with the lowest orthotropy ratios, while quite dispersed results are shown as this ratio reaches the value of 20.In particular, very incorrect results are given by ESL, FSDT (shear correction factor ), and HSDT for E1/E2 values greater than 30 and by theories MHR, MHR4, MHWZZA, MHWZZA4, HWZZM0, HWZZMB2, and HWZZMC2.Adaptive theories ZZA, ZZA*, HWZZ, HWZZM, and HWZZM* instead always provide results that are in very well agreement with the exact solution, irrespective of the orthotropy ratio.

Fundamental Frequency and Modal Displacements and Stresses of a Simply Supported, Anti-Symmetric Cross-Ply Plate
The fundamental frequency of [0/90/0/90], simply supported, cross-ply square plates with a length-to-thickness ratio of five [4], is examined in Table 8 for increasing the values of the degree of orthotropy E 1 /E 2 ranging from three to 40 (case c1).Comparisons are given with the results by a global-local higher-order theory [4] and by Matsunaga [56], which were obtained using a higher-order ESL.It can be seen that all of the theories of this paper provide fairly accurate predictions of the fundamental frequency in correspondence with the lowest orthotropy ratios, while quite dispersed results are shown as this ratio reaches the value of 20.In particular, very incorrect results are given by ESL, FSDT (shear correction factor π 2 /12), and HSDT for E 1 /E 2 values greater than 30 and by theories MHR, MHR4, MHWZZA, MHWZZA4, HWZZM0, HWZZMB2, and HWZZMC2.Adaptive theories ZZA, ZZA*, HWZZ, HWZZM, and HWZZM* instead always provide results that are in very well agreement with the exact solution, irrespective of the orthotropy ratio.
The through-thickness variation of modal displacement u ς , transverse shear σ ας , and transverse normal σ ςς modal stresses are reported in Figures 6 and 7 for ratios E 1 /E 2 of 3 and 40, respectively.Note that u α is correctly captured by all of the theories only for E 1 /E 2 = 3, so just the results for this highest ratio are reported.The results by HWZZMB2 and by HWZZM0 are only reported for E 1 /E 2 = 40, as they are inaccurate for all of the other cases.All of the other theories provide a quite accurate prediction of the in-plane displacement as well as of σ ας , with the only exceptions being MHR and MHR4, while there is a bigger scattering of results regarding σ ςς .
The fundamental frequencies for case c2 with an orthotropy ratio E 1 /E 2 of 30 and a length-tothickness ratio of 10/3 are reported in Table 9. Modal displacements and modal in-plane and transverse shear and normal stresses for this case are shown in Figure 8, along with the results by 3D FEA, Zhen and Wanji [4], and Matsunaga [56].These results show that MHWZZA, MHWZZA4, MHR, and MHR4 overestimate the fundamental frequency, while MHWZZA, MHWZZA4, MHR, and MHR4 incorrectly predict the in-plane displacement and stress, while MHR± and MHR4±, whose slope is defined on a physical basis, obtain better results than their counterparts MHR and MHR4 with the slope assumed a priori.
A rather large dispersion of results is shown for u ς and for the transverse normal stress, while the transverse shear stress is erroneously provided by HRZZ, HRZZ4, MHWZZA, MHWZZA4, MHR, and MHR4 only across the first layer.Again, in this case, with the slope being defined on a physical basis, HSDT and FSDT (the latter uses π 2 /12 as the shear correction factor in order to improve the accuracy of the results) appear inadequate to perform the analysis.

Fundamental Frequency of a Cross-Ply Plate with Different Thickness Ratios and Boundary Conditions
The fundamental frequency of a [0/90/0] cross-ply square plate with a thickness ratio ranging from 4 to 100 and simply supported (case d1), all clamped, or simply supported edges and clamped on opposite sides (case d2), are reported in Tables 10 and 11, respectively.The results are compared

Fundamental Frequency of a Cross-Ply Plate with Different Thickness Ratios and Boundary Conditions
The fundamental frequency of a [0/90/0] cross-ply square plate with a thickness ratio ranging from 4 to 100 and simply supported (case d1), all clamped, or simply supported edges and clamped on opposite sides (case d2), are reported in Tables 10 and 11, respectively.The results are compared to those by Di Sciuva and Icardi [50] (RFSDT, RHSDT, theories and RHQ40 elements) and to 3D FEA [48], as well as to those by FSDT and HSDT.
The results of case d1 show a very well agreement among the theories each other, as well as with 3D FEA reference solutions when the length-to-thickness ratio increases, because the layerwise effects become less important, while remarkable differences are shown for the thickest cases.All of the theories appear to be accurate except for HWZZM0, MHWZZA, and MHWZZA4 (even for moderately thin plates), FSDT (shear correction factor 5/6), and HSDT (the latter two are accurate only for the lengths-to-thickness greater than 10).The results for the different boundary conditions considered (case d2), namely clamped and simply supported edges (CSCS) and all-clamped edges (CCCC), are reported in Table 11 for a length-to-thickness ratio Lx/h = 10, along with the results for simply supported edges (SSSS).It is noted that only ZZA, ZZA*, HWZZ, HWZZM and HWZZM* give predictions of the fundamental frequency that are always very well in agreement with 3D FEA, while all of the other theories are less accurate.
It is also noted that the errors are greater for CSCS and CCCC, which therefore prove more problematic than SSSS because the mechanical boundary constraints are more difficult to satisfy identically for clamped edges.In this case, FSDT (shear correction factor 5/6) and HSDT do not provide valid results for any boundary condition.Therefore, it is deduced that for cases d1 and d2 FSDT and HSDT are unsuitable, similarly to other theories with a fixed representation.
So, it is confirmed also in this case that only adaptive theories whose coefficients can be redefined across the thicknesses (ZZA, ZZA*, HWZZ, HWZZM, and HWZZM*) obtain always accurate results with low computational cost (see Section 5.10), and then, they should be preferred in the applications.

First and Higher-Order Free Vibration Frequencies of Potpourri Cases
Here, the vibration behaviour of plates made of different materials and with different boundary conditions is discussed.Tables 12-14 report the first five frequencies for a [0/90/0] square plate (the intermediate layer has a thickness of h/2, while the outer ones have a thickness of h/4) here referred to as case e1, which is retaken from Kapuria et al. [5] (ZIGT and TOT theories).This plate is simply supported along the edges parallel to the y-axis and free at the other two edges (SFSF).Length-to-thickness ratios of 5, 10 and 20 are considered.
Again, adaptive theories ZZA, HWZZ, ZZA*, HWZZM, ZZA*, and HWZZM* appear as the most accurate among those considered, irrespective of the mode examined.On the contrary, the other theories exhibit errors that grow with the frequency order and with the increasing thickness.Only ZZA, HWZZ, ZZA*, HWZZM, ZZA*, HWZZM*, HRZZ, HRZZ4, and MHWZZA4 obtain quite accurate results for length-to-thickness ratios of 5, while MHR, MHR4, MHR±, MHR4±, HWZZMA, and HWZZM0 can't get the fourth and the fifth frequencies for the intermediate length-to-thickness ratio of 10.Instead, for L x /h = 20, all of the theories except FSDT (shear correction factor π 2 /12) and HSDT give accurate predictions of frequencies, as the layerwise effects tend to wears off, even if contrary to what is claimed in the literature, this does not always occur in all of the examined cases.assuming a linear in-plane displacement across the thickness and a uniform transverse displacement.Errors less than 2% are shown by adaptive theories ZZA, ZZA*, HWZZ, HWZZM, and HWZZM* with respect to 3D FEA, whereas the other theories exhibit larger errors that increase with the order of frequency that is considered.However, guessing the shear correction factor, accurate results can be achieved by FSDT (in the current case with a shear correction factor of 5/6), as also shown in [56], despite the very simple kinematics.This gives a reason at least for the case that is currently examined, to those who consider ESL suitable for dynamic analysis.However, it could be argued that an appropriate shear correction factor could not easily be chosen in the industrial applications.Table 16 reports the first six frequencies of case e3, which concerns a simply supported (SSSS) [0/90/core/0/90] soft-core sandwich plate with a length-to-thickness ratio of 10.It should be noted that for this case, it is necessary to enforce two additional conditions across the upper layer in the adaptive theories, which means the first two local equilibrium equations; otherwise, the accuracy drops.Having done this, again, all of the higher-order adaptive theories obtain accurate results for all of the frequencies.
Also HRZZ, HRZZ4, HWZZMB, HWZZMC, and HWZZMC2 provide quite accurate results, while other lower-order theories are inaccurate.As in many other cases, MHR and MHR4 are not adequate, because Murakami's rule is not respected, so MHR± and MHR4±, whose slope sign is determined on a physical basis, obtain better results.Also, in this case, FSDT (shear correction factor 5/6) and HSDT (no need of shear correction factor) obtain inaccurate results.[53], who analysed it via finite elements.The interesting aspect of this case is that great cross-section deformations are already involved by the first four frequencies, which could undermine the concept of the plate on which the theories that are considered in this paper are based.Indeed, deformations from the first to the fourth frequency are respectively: a bending mode, a bending/torsional mode, a torsional mode, and finally an axial/shear mode.
The structure is 200-mm thick, while its length-to-thickness and length-to-side ratios are 10.In this case, a suited choice of even and odd trial functions should be made to obtain accurate results by the present theories, and a sufficiently high expansion order of the representation must also be considered, as indicated in Table 2.
The results for this case, which are reported in Table 17, where they are compared to the present FEA results and the ones from [53], show that adaptive ZZA, ZZA*, HWZZ, HWZZM, and HWZZM* theories only commit an error of the order of 2% or less, albeit their representation order across the thickness and their number of unknowns are lower of those of the theory that was used to construct finite elements in [53].Since there is only one interface, HWZZMA, HWZZMB, HWZZMB2, HWZZMC, and HWZZMC2 achieve the same accuracy of HWZZM, whereas all of the other lower-order theories predict wrong frequencies, which is a sign that a refined kinematics is required.In this case, the results by FSDT and HSDT aren't reported as being totally wrong.For the purpose of checking which theory is able to effectively capture strong 3D effects related to transverse normal and shear deformability, case g is examined.It concerns a simply supported cylindrically bent sandwich plate with stiff faces and a compliant core that has not yet been considered in the literature.Each face has a thickness of 0.2 h, and is made of three layers having different thickness and properties, as indicated in Table 3, while the core is 0.6-h thick.The first and the third face layers proceeding from the outside toward the inside are made of the same very stiff material (m2), while the interposed layer is made of a more compliant material in tension, compression, and shear (m1); finally, the core is made of the most compliant material (m3), as indicated in Table 4.A length-to-thickness ratio of five is considered.
Table 18 reports the first five free vibration frequencies for this case, as predicted by the theories of this paper and by 3D FEA.The results indicate that the first, second, and fifth frequencies represent the bending modes, while the third and fourth ones represent more interesting motions that occur across the thickness in a symmetrical manner with respect to the mid-plane.So, specific boundary conditions should be enforced to get rid of these modes.
Moreover, the fulfillment of local equilibrium equations should be enforced near the core interfaces in the adaptive theories, the only ones where it is possible to do it; otherwise, poor results similar to those of other theories are obtained.The results demonstrate the superior accuracy of the adaptive theories ZZA, ZZA*, HWZZA, HWZZM, and HWZZM*, which were obtained thanks to the imposition of these constraints, irrespective of whether they directly considered a piecewise transverse displacement or recovered the effects of the normal transverse deformation differently.Naturally, in this case, the gap with other theories is much more marked than in the previous cases because of the greater importance assumed by the transverse normal deformation, it being the one that was less accurately reproduced.In this case, HWZZMC and HWZZMC2 obtain inaccurate results, with their errors increasing with the frequency order, i.e., the third and fourth frequencies are progressively more inaccurate.From Figure 9, which reports the in-plane modal displacement, in-plane and transverse shear modal stresses for the first mode, and a minor subject to errors, it can be seen that all of the theories except HWZZMA, MHR, MHR4, MHR±, and MHR4± obtain quite accurate results, which is a sign that Murakami's theories are inappropriate for this case.However, it is noted that the greatest errors occur on the transverse modal displacement and the normal modal stress.Figure 10 reports the modal stresses and displacements for the third frequency by ZZA, ZZA*, HWZZA, HWZZM, HWZZM*, HWZZMC, and HWZZMC2, while the results by other theories are not reported as being too wrong.
are progressively more inaccurate.From Figure 9, which reports the in-plane modal displacement, in-plane and transverse shear modal stresses for the first mode, and a minor subject to errors, it can be seen that all of the theories except HWZZMA, MHR, MHR4, MHR±, and MHR4± obtain quite accurate results, which is a sign that Murakami's theories are inappropriate for this case.However, it is noted that the greatest errors occur on the transverse modal displacement and the normal modal stress.Figure 10 reports the modal stresses and displacements for the third frequency by ZZA, ZZA*, HWZZA, HWZZM, HWZZM*, HWZZMC, and HWZZMC2, while the results by other theories are not reported as being too wrong.The first six free vibration modes of a simply supported sandwich plate having the same thickness of faces and core of case e3 and the same orientation is studied assuming a thicker length-to-thickness ratio of five.With the same purpose, stiffer faces are considered, as indicated in    The first six free vibration modes of a simply supported sandwich plate having the same thickness of faces and core of case e3 and the same orientation is studied assuming a thicker length-to-thickness ratio of five.With the same purpose, stiffer faces are considered, as indicated in  The first six free vibration modes of a simply supported sandwich plate having the same thickness of faces and core of case e3 and the same orientation is studied assuming a thicker length-to-thickness ratio of five.With the same purpose, stiffer faces are considered, as indicated in Table 3, while the core has the same properties as case e3.These choices, which are distinctive of a benchmark that has never been investigated before, here called case h, enhance the layerwise effects, and consequently they should highlight the quite different behaviours of the theories.
Table 19 reports the first six free vibration frequencies predicted by the present theories and by 3D FEA for this new case.Similar to the previous case g, the first five frequencies are bending modes, while the sixth represents a motion that occurs in a symmetrical manner in the thickness direction.So, the same considerations about the constraints that must be imposed apply again, and consequently, it is still demonstrated that the material properties and thickness of constituent layers constitutes a strong discriminatory effect on the accuracy of the theories.
Considerations that are similar to those of the previous case g apply, because only adaptive theories are always in well agreement with 3D FEA results.All of the other theories, except HWZZMC and HWZZMC2, which obtain results accurate enough, give inaccurate predictions, especially for the sixth mode.

Blast Pulse Loading
In this section, two square plates with a different lay-up and subject to step and exponential blast pulse loadings are analysed.The first is a sandwich plate that is subject to a step pulse (case i1), which has a length Lx of 609.6 mm, its core is 25.4-mm thick, and its two faces, each one being a five-layer laminate, have a total thickness of 1.905 mm.The second is a laminated [0/90/0] square plate whose central layer is two times thicker than the outer ones, whose length is 2540 mm and whose overall thickness is 170 mm (case i2), which is subject to an exponential blast pulse.Such case studies are retaken respectively from Hause and Librescu [54], where the step blast pulse overpressure loading is described as: and from Librescu and Noisier [55], where the exponential blast pulse overpressure is considered as: where t p is 0.1 s.
The results for the first case i1 are reported in Figure 11 as the central plate deflection at the mid-plane normalized to the plate thickness.It can be seen that also in this case, only theories ZZA, ZZA*, HRZZ, HRZZ4, and HRZZ4* provide a correct time-variation of the deflection, which is in good agreement with 3D FEA.Notice that the results by [54] differ from the FEM results because of the lower-order model used therein.However, we see a similar behaviour where errors tend to disappear with respect to the present theories if a single half-wave in the x and y directions is considered instead of the expansion order that is reported in Table 2.As for the other cases considered previously, again, MHR and MHR4 appear to be inadequate because Murakami's rule is not respected, while their counterparts MHR± and MHR4±, slopes are computed on a physical basis, appear to be more accurate.
The results for the case i2 are reported in Figure 12, where again the transverse displacement is reported at z = 0 and at the center of the plate, and is still normalized to the plate thickness.Since now the layer has the same material properties and a symmetric stack-up, the layerwise effects fade, so the discrepancies between the results by theories also fade.Anyway, MHWZZA4, MHR4, MHR4±, HWZZMA, and HWZZM0 still appear inaccurate, whereas the adaptive theories once again obtain results that are in very good agreement with 3D FEA.Table 20, which reports results for specific instants (0.9 ms, 4.5 ms, 5.6 ms, 6.5 ms, and 7.4 ms for i1 and 3.9 ms, 8.0 ms, 19.96 ms, and 24.08 ms for i2) show that the erroneous predictions (both frequency and amplitude of response are wrong) given by HRZZ, HRZZ4, HWZZMB, HWZZMB2, and HWZZMC2 are erroneous, which had not been noticed for the initial instant (0.9 ms for i1 and 3.9 ms for i2, which corresponded to the occurrence of the first peak of ZZA in each of the two cases) considered in Figures 11 and 12).It is worth noting that in both cases i1 and i2, there is no detectable difference when the transverse displacement is evaluated in points across the thickness other than at the middle plane as reported in the tables and figures, because both structures are thin.However, splitting the core into two halves whose upper half is much more compliant than the lower one (case i3), similar as to when a face is damaged, and assuming Lx/h = 10 and a different orientation of layers of faces, the results of Figure 13 are obtained, which show visible differences between the results of the theories.Those assuming a uniform or a polynomial transverse displacement in this case don't account properly for this effect, so less accurate results are obtained, see e.g., MHR±, MHWZZA, MHWZZA4, HRZZ, HRZZ4, FSDT, and HSDT.
It is specified that a Newmark implicit time integration scheme was adopted for solving transient dynamic equations, since alternative explicit time integration schemes need extremely small time steps to be stable.However, for reasons of the stability of the algorithm, small time-steps are still required (30 µs) to limit the convergence and rounding errors.Since geometrical and material non-linearity are disregarded, the system to solve is a linear system, and the computational burden isn't adversely affected by such a small time step.

Computational Effort of Theories
Tables 21 and 22 reports the calculation times that are necessary to solve each of the benchmarks considered by the examined theories, which being based on the same five d.o.f., therefore have a memory storage occupation that is practically indistinguishable from one another.As closed-form solutions are considered even when other researchers recoursed to FEA due to the complexity of the solutions in the cases examined, calculation times are very short for all of the theories; that is, they remain comparable to those of FSDT and HSDT.
This testifies to the efficiency of the present adaptive and higher-order, because they require just a reduced expansion order, both with regard to the in-plane and through-thickness representation to achieve accurate results for all of the challenging benchmarks examined.So, it can be said that a level of accuracy that is comparable to that of the FEA has been obtained with a lower computational burden.However, FEA remains indispensable for solving the problems of industrial complexity, while the preliminary parametric studies can be performed as in this paper.
It appears that the MHR and MRH4 theories have the lowest processing time out of all of the theories, but this advantage is totally negated because they provide inaccurate results whenever strong layerwise effects rise, because in these cases, the slope varies differently from what is expected by Murakami's rule.Although slightly more expensive, MHR± and MHR4± obtain often rather accurate results, since their slope sign of displacements is decided on a physical basis.HRZZ and HRZZ4 Processing time reported include symbolic computation.For calculation, a laptop computer with quad-core CPU@2.60 GHz.64-bit operating system and 8.00 GB RAM was used.
Higher-order theories HWZZ, HWZZM, ZZA*, and HWZZM* provided rather more accurate results in all of the cases examined, and required a little longer processing time than HRZZ and HRZZ4.However, it is noted that HWZZM* with a priori assumed zig-zag amplitudes requires 20% less processing time than HWZZ, but they don't appear to be the most accurate theories.In particular, HWZZMA, HWZZMB, HWZZMB2, HWZZMC, and HWZZMC2 appear to be inadequate in many cases.ZZA, HWZZ, HWZZM, ZZM*, and HWZZM* don't qualify for the lowest calculation time between all of the theories, appear to be the most efficient theories, and thus are preferred in the applications, which are always very accurate and still have affordable costs.However, the best of such adaptive theories from this point of view turns out to be ZZA*, which has a slightly lower calculation cost.

Concluding Remarks
Various displacement-based and mixed zig-zag theories, which differ in the layerwise functions and in the scheme of the through-thickness representation of the displacements that are used, have been applied to investigate the free vibration behaviour and the response of blast-loaded laminated and sandwich plates with different length-to-thickness ratios, lay-ups, constituent materials, and boundary conditions.To homogenise the results, they are compared using the same type and order of representation as closed-form solutions, with the appropriate trial functions being selected for each benchmark to minimize the expansion order.The intended aim is to evaluate the merits and drawbacks of theories in order to establish which are significantly much more accurate and efficient.
The numerical applications show the importance of very accurately accounting for the transverse normal deformability whenever the layers have non-uniform mechanical properties and a different thickness.Indeed, adaptive zig-zag theories whose layerwise contributions identically satisfy interfacial stress constrains and whose displacement fields are redefined for each layer prove superiority.ZZA* theory shows that the choice of zig-zag functions is immaterial whenever the coefficients of displacements are recomputed across the computational layers.In this context, zig-zag functions can even be omitted, as the stress continuity constraints can be enforced in order to define the coefficients of displacement fields in a more computationally efficient way.
The accuracy of results is shown to be independent of the choice of zig-zag functions for ZZA*, but this result is extensible to all of the theories that in the same way provide a redefinition of the coefficients of displacements across the thickness, so as to satisfy the physical constraints.Vice versa, the theories whose coefficients of displacements are fixed fail to be accurate whenever strong layerwise effects rise or there is a strong transverse anisotropy, since finding a kind of fixed representation that is always suitable is impossible, unless a very high order of representation is used.That is the opposite of what this paper sets out, which is wanting to find accurate solutions at a low cost.Indeed, the accuracy of theories with a fixed representation appears to be largely case-dependent.Mixed theories such as MHWZZA and MHWZZA4 based on Murakami's zig-zag function (as well as all those for which zig-zag amplitudes are a priori assumed) are often proven inaccurate, although not in all cases, even though they benefit from strain and stress fields by adaptive theories.The same happens even when the slope sign of displacements at interfaces is established on a physical basis, at least for the low orders of the in-plane and through-thickness representation that are considered in this paper, which however allow the adaptive theories to be already very accurate.Anyhow, it is not easy to discern for which cases the limiting assumptions of such theories do not have weight.Therefore, it is not possible to establish a general rule, although the results undoubtedly show that the theories accounting for layerwise effects without the determination of zig-zag amplitudes on a physical basis cannot provide an adequate level of accuracy with the low expansion orders that are considered in this paper.
A simplified uniform or polynomial representation of the transverse displacement is shown to be ineffective even when the strain and stress fields are retaken from other more accurate structural models, such as for MHWZZA.In particular, FSDT and HSDT theories are proven to be inaccurate in the majority of the examined cases.
Although the adaptive theories whose coefficients of displacements are redefined across the thickness do not get the lowest processing time, they were proven to be the efficient ones, since they

Figure 1 .
Figure 1.Case a1: Propped cantilever sandwich plate in cylindrical bending (Lx/h = 5.714) under a uniform loading on the top layer.

Figure 1 .
Figure 1.Case a1: Propped cantilever sandwich plate in cylindrical bending (L x /h = 5.714) under a uniform loading on the top layer.

Figure 2 .
Figure 2. Case a2: Propped cantilever sandwich plate in cylindrical bending (Lx/h = 20) under a uniform loading on the top layer.

Figure 2 .
Figure 2. Case a2: Propped cantilever sandwich plate in cylindrical bending (Lx/h = 20) under a uniform loading on the top layer.

Figure 3 .
Figure 3. Case b3: Normalized modal displacements and stresses of a simply supported sandwich plate in cylindrical bending (Lx/h = 4).

Figure 3 .
Figure 3. Case b3: Normalized modal displacements and stresses of a simply supported sandwich plate in cylindrical bending (L x /h = 4).

Figure 4 .
Figure 4. Case b3: Normalized modal displacements and stresses of a simply supported sandwich plate in cylindrical bending (Lx/h = 10).

Figure 4 .
Figure 4. Case b3: Normalized modal displacements and stresses of a simply supported sandwich plate in cylindrical bending (Lx/h = 10).

Figure 4 .
Figure 4. Case b3: Normalized modal displacements and stresses of a simply supported sandwich plate in cylindrical bending (Lx/h = 10).

Figure 5 .
Figure 5. Case b3: Normalized modal displacements and stresses of a simply supported sandwich plate in cylindrical bending (Lx/h = 20).

Figure 5 .
Figure 5. Case b3: Normalized modal displacements and stresses of a simply supported sandwich plate in cylindrical bending (Lx/h = 20).

Figure 9 .
Figure 9. Case g: Normalized modal displacements and stresses of a simply supported sandwich plate in cylindrical bending (Lx/h = 5, first mode).

Figure 10 .
Figure 10.Case g: Normalized modal displacements and stresses of a simply supported sandwich plate in cylindrical bending (Lx/h = 5, third mode).

Figure 9 .
Figure 9. Case g: Normalized modal displacements and stresses of a simply supported sandwich plate in cylindrical bending (Lx/h = 5, first mode).

Figure 10 .
Figure 10.Case g: Normalized modal displacements and stresses of a simply supported sandwich plate in cylindrical bending (Lx/h = 5, third mode).

Figure 10 .
Figure 10.Case g: Normalized modal displacements and stresses of a simply supported sandwich plate in cylindrical bending (Lx/h = 5, third mode).

Figure 12 .
Figure 12.Case i2: Normalized transverse displacement of a simply supported laminated plate under an exponential blast pulse loading (Lx/h = 14.941).

Figure 11 .
Figure 11.Case i1: Normalized transverse displacement of a simply supported sandwich plate under a step blast pulse loading (Lx/h = 20.8696).

Figure 12 .
Figure 12.Case i2: Normalized transverse displacement of a simply supported laminated plate under an exponential blast pulse loading (Lx/h = 14.941).

Figure 12 .
Figure 12.Case i2: Normalized transverse displacement of a simply supported laminated plate under an exponential blast pulse loading (Lx/h = 14.941).

Figure 13 .
Figure 13.Case i3: Normalized transverse displacement of a simply supported sandwich plate under a step blast pulse loading (Lx/h = 10).

Figure 13 .
Figure 13.Case i3: Normalized transverse displacement of a simply supported sandwich plate under a step blast pulse loading (Lx/h = 10).

Table 3 .
Data of cases.
*** Not provided by the theory.
*** Not provided by the theory.