Equivalent single-layer models in deformation analysis of laminated multilayered plates

The performance of selected equivalent single-layer (ESL) models is evaluated within several classical benchmark tests for linear static analysis of multilayered plates. The authors developed their own finite element software based on the first-order shear deformation theory (FOSD) with some modifications incorporated including a correction of the transverse shear stiffness and an application of zigzag-type functions. Seven different ESL models were considered in the study; beside the classical FOSD model, there were three FOSD models with various transverse shear corrections and three ESL models enhanced by the application of zigzag functions. In addition, particular attention was paid to investigation of differences related to the “soft” and “hard” variants of the simple support.


Introduction
A typical multilayered plate as considered in this paper can take the form of a sandwich panel or a thin composite laminate [1]. In both cases, the plate is constructed as a multilayered panel manufactured of several laminas bonded together with some adhesives. The fundamental idea of creating a multilayered construction is to combine together several different layers to obtain a new heterogeneous product with properties that are superior to those of individual laminas. A detailed representation of the multilayered plate's behavior can be quite complicated; therefore, for practical applications, one has to find a rational compromise between the requirements of precise description and realistic costs of calculations [2]. Consequently, the majority of computational models of multilayered panels is constructed within the macro-mechanical approach with the plate considered as a two-dimensional (2D) object (refer [3,4] for a comprehensive survey of computational models). Such models are usually prepared according to an appropriate lamination theory where the essential postulations describing the assumed arrangement and interaction of all components are undertaken. Typically, the laminated panel is treated as a stack of several laminas perfectly bonded together without any gaps, and without any possibility of a slip or a separation.
In the most popular equivalent single-layer (ESL) approach, the entire laminate is represented by a singlelayer panel with macro-mechanical properties estimated as a weighted average of the mechanical properties of all laminas. More complex computational models of multilayered panels are related to layer-wise (LW) theories [4,5] (known also as the discrete-layer formulations [6]) where each layer is considered separately with appropriate boundary conditions at the layer interfaces. Since each lamina is treated individually, the number of unknowns in LW theories is directly related to the number of layers, N. The most comprehensive 3D modeling requires an even larger number of unknowns; therefore, an application of such an approach is practically limited to those selected regions of the structure, where a necessity of such detailed modeling is forced by some substantial irregularities, i.e., connections with other objects.
It has been shown in many examples described in the literature (see e.g., [6,7]) that ESL models can be very effective in the investigation of a global response of multilayered panels of moderate thickness, despite substantial simplification assumptions made in that approach. Unquestionably, their main advantage lies in a restricted size of the resulting equation system, with the number of unknowns being independent of the number of layers. On the other hand, it is commonly known that the performance of most ESL models in the deformation analysis of multilayered panels significantly deteriorates for small values of the length-tothickness ratio. For that reason, many different ideas for a possible improvement in the ESL model formulation were presented in the literature during the last 50 years [3,4,8,9]. Nevertheless, to the best knowledge of the authors of the present report, the FOSD formulation is still the most common pattern of many FE models offered in commercial software for FEA of multilayered panels. One can imagine that many users of everyday use FEA systems would be extremely interested in knowing the limits and the level of error of those popular models. In search of the answer to this question, the authors of the present research examined the performance of the selected seven different ESL models within several standard benchmark tests of a linear static analysis of multilayered plates: Beside the classical FOSD model, three FOSD models with various transverse shear corrections were included in the study together with three ESL models enhanced by the application of zigzag functions. The results of the finite element analysis performed with the examined models were confronted with the 3D elasticity theory solutions.

Equivalent single-layer models
It is quite obvious that straightforward application of a regular 2D plate model in the analysis of a multilayered composite laminate or a sandwich plate is associated with a considerable approximation. As a result of significant shear deformability and mainly due to intrinsic discontinuity of material properties at the interfaces between contiguous layers, 1 the deformation profile of a multilayered plate takes the shape of a polygonal chain with rapid changes in slopes at the layer interfaces. This phenomenon is commonly classified as the zigzag effect [8]. A comparatively complex profile appears also for transverse shear stress components when the interlayer shear stress continuity conditions have to be satisfied for a multilayered structure with different material properties of individual layers. One should be aware that every attempt to model the behavior of such a complex structure with a simplified ESL approach is associated with unavoidable error of the attained results where the degree of this inaccuracy strongly depends on the particular considered case.
Most of the ESL models are based on the first-order shear deformation (FOSD) theory where a linear variation of a displacement field over the thickness of the panel is assumed. The earliest formulation of the FOSD theory for multilayered plates was published by Yang et al. [10]. The main advantage of a resulting computational model lies in a limited number of unknown parameters and their clear physical interpretation. On the other hand, due to a constant distribution of the transverse shear strains across the plate thickness, such an approach cannot provide the equilibrium of shear stress at the layer interfaces. As a consequence, computational models based on the FOSD theory require an appropriate shear correction [11].
Some improvements in ESL models can be obtained within the higher-order shear deformation (HOSD) theory by applying a power series expansion of the displacements and strains with respect to the thickness coordinate [12,13]. Such a strategy can eliminate the necessity of an implementation of a shear correction and provide a substantial improvement in results for thick panels; however, the number of variable parameters significantly increases in HOSD models when compared with FOSD ones and also some questions appear regarding a physical interpretation of certain new variables.
Another possible enhancement of the ESL models can be achieved by inserting some warping functions that can represent the zigzag feature of the deformed profile with a different slope in each layer. One of the pioneering zigzag models for multilayered plates is attributed to Ambartsumyan [14]; more examples were described in an extensive historical review on zigzag model developments published by Carrera [8].

Plate models considered in the current study
It is assumed that the considered multilayered plate can be analyzed as a layered structure built of orthotropic laminas. By disregarding the transverse normal stress (z-direction) component, the following constitutive relation holds for an individual lamina (cf. [1]): In particular, the laminate can be composed of several fiber-reinforced composite laminas with varying orientation of reinforcement in each layer. Then, a linear elastic constitutive relation for orthotropic material is given in the system of material axes (a, b, c) ( Fig. 1): where E i , ν i j and G i j stand, respectively, for appropriate Young's modulus, Poisson's ratio and shear (Kirchhoff's) modulus, along with the relation ν ab E b = ν ba E a . As was already signaled, the transverse shear stiffness in the FOSD model requires a special treatment; therefore, the corresponding elements of the constitutive matrices in (1) and (2) are marked with a tilde (˜ ). The simplest recipe that can be applied here means just a multiplication of the shear moduli G bc and G ca by the transverse shear correction factor κ: This strategy directly follows the idea of Reissner [15] who introduced κ = 5/6 in his theory for homogeneous plates. 2 A proper estimation of shear correction for laminated composites would be much more complicated than for homogeneous panels; therefore, the local transverse shear correction factor κ as introduced in (3) at the layer level should be treated as an alternative to the global (cross-sectional) transverse shear correction factors considered in the following part of the current report.
With a usual transformation rule for the stress and strain vectors between the material axes (a, b, c) and the coordinate system (x, y, z): the constitutive matrix used in (1) can be obtained as ⎡ The linear strain components are expressed as appropriate derivatives of the displacements: Within the ESL formulation, the entire multilayered laminate is represented by a single-layer plate with appropriately established macro-mechanical properties. In the classical FOSD formulation [10], the displacement distribution through the plate thickness is taken as where u, v and w are usual translational displacement components at the mid-surface, with ϕ x and ϕ y standing for rotational components. The deformation profile (8) can be furtherly enriched by taking into account the through-the-thickness warping effects represented by appropriate zigzag functions φ (k) α (z), α = x, y ( Fig. 2): with ψ x and ψ y being two supplementary displacement parameters that can be understood as the weightedover-thickness amplitudes of the zigzag effects. One of the most popular concepts of the warping function is known as the zigzag function of Murakami [17]: where h k is the thickness of the kth layer and ζ k is the local coordinate of the kth layer: Quite recently, Tessler et al. [18] proposed an alternative concept where the changes in slopes at the interfaces between two contiguous layers were related to the differences in the transverse shear stiffness of those layers: with N standing for the total number of layers and k ∈ 1, N being the number of the current layer; Q (k) αα represents the diagonal coefficients of the constitutive matrix (6) for the kth layer (Q x x = c 55 and Q yy = c 44 ). The weighted-average transverse shear stiffness coefficients G α are established as By putting (9) into the first three equations of (7), the in-plane strain components can be presented as The vector containing the mid-plane bending strain components can be related to the vector of the displacement components as With (7) and (9), the transverse shear strain components can be presented as: One can notice that only the last components on the RHS of the above equation can be related to the coordinate z (depending on the choice of the zigzag function). Moreover, the derivatives of the zigzag functions defined, either by (10) or by (12), with respect to the coordinate z would be constant over the thickness of the layer: (a) For the zigzag function of Murakami (10): (b) For the zigzag function of Tessler et al. (12): As a result of that last observation, one can conclude about the transverse shear strains keeping a constant value over the thickness of the layer: The equilibrium equations of the FE model in the most popular displacement formulation can be derived from the Principle of Virtual Displacement (PVD), stating that the virtual work of external forces, δW ext , must be balanced by the internal virtual work, δW in , with [N(x, y)] standing for shape functions and {q} being the vector of nodal displacements. Then, with the strain and stress fields established with the standard formulas: one can arrive at the classical form of the linear static FEA governing equation with the stiffness matrix and the global vector of nodal forces It should be noted that the sequence of expressions (23-28) presents just a schematic picture of the whole FE procedure; in particular, the relation between strains and nodal displacements that is implemented in computer codes is much more complex than shown in (24)(25): firstly, the bending (14) and transverse shear (20) parts of the strains are handled separately, and secondly, a distinctive rearrangement of matrices B allowed for a much effective evaluation of integrals in (27-28) by taking advantage of a through-the-thickness pre-integration. By performing the through-the-thickness integration of the appropriate stress components, one can obtain the effective stress resultants, among others the transverse shear forces: related to the transverse shear strains by the generalized constitutive relation: whereã i j represents the global (cross-sectional) transverse shear stiffness of the whole laminate 3 that can be estimated with the application of the global composite transverse shear correction factors The global composite transverse shear correction factors k 1 and k 2 used in (31) depend not only on material properties but also on such geometrical characteristics of the laminate as the stacking sequence of the layers and their ply angles. Values of the composite transverse shear correction factors can be determined by matching the transverse shear strain energy predicted by the FOSD plate model with that obtained from the three-dimensional elasticity theory (see e.g., [11,[19][20][21]). A slightly enhanced strategy was implemented by Noor and Peters [22] who calculated the transverse shear correction factors for multilayered cylindrical panels by means of the predictor-corrector approach. It is quite obvious that instead of calculating the shear correction factors one can simply apply "corrected" transverse shear stiffness in the FOSD model-this tactic was described, e.g., in [23,24]. In this context, one should also mention a recent paper of Gruttmann and Wagner [25] where the global transverse shear correction factors for layered panels were obtained as by-products in the analysis performed with a more advanced approach utilizing a coupled global-local shell model with layer-wise warping displacements.
By assuming a different weak formulation of the boundary value problem from the PVD, e.g., a mixed variational principle or hybrid stress formulation, one can introduce an independent interpolation of the transverse stress components. This strategy allows for an application of arbitrary profiles of the transverse stress components, also those that ensure the interlayer shear stress continuity. In practice, it is enough to use a mixed variational formulation with an independent interpolation of transverse stress components as proposed by Reissner [26]. This approach known in the literature as Reissner's Mixed Variational Theorem (RMVT) was used as a base to develop Murakami's theory for laminated plates [17,27], which accounts for the zigzag effect. By assuming a piecewise linear distribution of the in-plane displacements, the authors of [27] adopted in fact the FOSD hypothesis for each individual layer of the multilayered plate. However, due to implementation of the RMVT they could invoke the quadratic distribution of the transverse shear stress components in each layer and simultaneously provide continuity of inter-laminar stresses. This approach was followed i.a. by Carrera [28] (see also [5,29,30]).
According to RMVT, the variation of the internal work in the laminated plate composed of N layers can be expressed as where ε stands for strains, σ denotes stresses, k is the layer number (layers are numbered sequentially, starting at the bottom surface of the plate), h k means the thickness of the kth layer, A represents the mid-surface area; the subscripts B and S refer to bending and shear components, respectively, whereas an additional subscript (a) is used to indicate the assumed field of transverse shear stresses and the related field of transverse shear strains (via Hooke's law). The first and the second integrand in (32) represent the internal energy due to the bending and the transverse shear of the plate, respectively. They both correspond in a variational manner to the equilibrium condition, whereas the third term implies variationally the compatibility of the transverse shear tractions.
It is worth to notice that the transverse shear stress σ S(a) in (32) can adopt any possible distribution across the plate thickness. Following Murakami [17,27], a quadratic distribution of transverse shear stresses is assumed in each layer (cf. Fig. 3) according to and By invoking the stress continuity requirements at the layer interfaces together with additional assumptions of the vanishing transverse shear stress at the top and the bottom surfaces: the stress variables can be eliminated from the set of unknowns, and the final formulation can be treated as a displacement-based model with only seven kinematical unknown parameters as shown in (9) (cf. [28,31]). Finally, the set of computational models applied in the current study for the purposes of comparative analyses includes seven formulations presented in Table 1.

Numerical examples
All own results presented in the current report have been obtained by using 9-node Lagrangian isoparametric elements with selectively reduced integration (9-SRI) [32]. For each example, a convergence test has been performed to find the required density of finite elements, which provided established values of the control deflections. Only small displacement analysis has been considered in the current study. The plate deflections calculated with the considered computational models will be presented in diagrams as normalized with respect to the chosen reference solution.
To the authors' knowledge, the first 3D elasticity solutions for multilayered plates were presented in 1970 by Pagano [33], who analyzed simply supported rectangular plates under a transverse load with bi-sinusoidal distribution (Fig. 4): That pattern will be replicated in all the subsequent examples. The first example to be considered in the current study is a square (a = b) 3-layered composite cross-ply (0 • /90 • /0 • ) plate. All layers have the same thickness with the following parameters of the orthotropic material The central deflection of the plate was calculated for a varying value of the a/h ratio. Due to a double symmetry of the problem, only one quarter of the plate is analyzed numerically with the following boundary conditions: One should mention that above constraints fully match the set "SS-1" of simply support conditions defined by Reddy [2]. For the reason that Pagano [33] presented deflection solutions in the form of graphs only, all digital data necessary for the reference purposes for the considered plate were recovered with the help of [34,35]. The errors in the estimation of the central deflection with respect to the 3D elastic solution of Pagano [33] are presented in Fig. 5 in the form of graphs for varying value of the length-to-thickness ratio.
One can observe that the higher inaccuracy is related to the "pure" FOSD model (FOSD_0), and model FOSD_1 with the local shear correction factor κ = 5/6 performs only slightly better. It is almost impossible to distinguish among results obtained with both FOSD models with global shear correction (FOSD_2 and FOSD_3) which share a common graph in Fig. 5. Series FOSD_2 and FOSD_3 show a relatively small error for thin and moderately thick (a/ h ≥ 10) plates. As one could expect, the best performance is revealed by the zigzag models; however, ZZ_2 based on PVD cannot provide accuracy comparable with that offered by the both zigzag models composed with the use of RMVT (ZZ_1 and ZZ_3 series represented in Fig. 5 by a joint  line).
The boundary conditions described by (39) should be classified as hard simply support (cf. [36]) in contrast to soft simply support boundary conditions defined with: Central deflecƟon: Error with respect to 3D elasƟcity soluƟon [33] a/h The consequences of introducing such constraints release can be observed in Fig. 6 where results for boundary conditions (40) have been marked with dashed lines. Unsurprisingly, with less restrictive boundary conditions the central deflection of the plate noticeably increases and, consequently, one can observe a slight "improvement" in those models which appeared as being too stiff in the previous test. Even though the errors of models FOSD_0 and FOSD_1 slightly decreased, their performance is still far from a satisfactory one. On the other hand, the introduced softening of the computational model resulted in a spectacular "improvement" in the solution obtained with the model ZZ_2. Accordingly, the softer support conditions in the modified problem developed greater than before errors of models FOSD_2 and FOSD_3, similarly as those of ZZ_1 and ZZ_3.
In the authors' opinion, the issue of a proper description of boundary conditions for simply supported plates is worth further reflections. It can be observed that various researchers (see, e.g., [37]) accepted recommendations of Häggblad and Bathe [38] to apply soft simply support conditions in FEA of plates, as the unconstrained rotation normal to the boundary allows for formation of the edge effects that are suppressed when hard-soft support conditions are applied (cf. [36]). However, one should notice that boundary conditions specified quite precisely for the considered example by Pagano [33] match the hard simply support (39) not the soft simply support settings (40).

Simply supported cross-ply
The performance of the tested computational models is further examined for the rectangular (b = 3a) threelayer simply supported plate under bi-harmonic load (38) with the same lamination scheme and material properties as those in the previous example.
Here, we start again with the comparative study presented for the hard simply support. The error of the central deflection related to the exact 3D elastic solution of Pagano [33] (digitally verified also with [39]) is presented in Fig. 7 as a function of the length-to-thickness ratio.
At the first glimpse, one can notice a qualitative similarity between Figs. 5 and 7; however, the range of divergences with respect to the exact 3D elastic solution is bigger for the rectangular plate than was encountered in Sect. 4.1 for the square plate. It was again impossible to notice the difference between results of FOSD_2 and FOSD_3, or between those for ZZ_1 and ZZ_3.
In the next step, the less restrictive boundary conditions (40) were included in the comparative analysis (dashed lines in Fig. 8). While looking at the graphs in Fig. 8, one can be much surprised because this time the softening of constraints resulted in much less spectacular effects than in the case of a square plate. In particular, the results of the ZZ_2 did not experience such impressive "improvements" as could be observed in the previous example for a square plate with the soft simple support.  [40] who provided also the exact 3D elastic solution which will serve as the reference exact solution in our study (with numbers confirmed also by Zenkour [39]). All layers are made of the high modulus graphite/epoxy composite with the same material parameters as considered in Sect. 4.1. The thickness of each odd layer is 0.1h, whereas every even layer is 0.125h thick.
In the first computation sequence, the bending of the plate under bi-harmonic load (38) was analyzed for various values of the a/h ratio, assuming boundary conditions specified in (39). The graphs presenting the error of examined models in the central deflection estimation with respect to the reference solution [40] are presented in Fig. 9.
While analyzing the graphs in Fig. 9, one can notice that the concept of the global shear correction introduced in the models FOSD_2 and FOSD_3 allowed for achieving quite a good estimation of the central deflection with the error not exceeding 4% for the thickest considered plate (h = 0.25a). It is quite interesting that none of the other models examined in our study can provide a similar accuracy for the considered example. In particular, one could expect a better performance of the zigzag models ZZ_1, ZZ_2 and ZZ_3 (the first and the third series represented by a joint graph in Fig. 9); however, this offers a slightly better accuracy than the FOSD model with the local shear correction factor κ = 5/6 (FOSD_1). Unsurprisingly, the highest errors are related again to the "pure" FOSD model (FOSD_0).
In the next computational round, the comparative analysis was repeated for the soft simply support conditions (40). The results calculated in that second part are marked with dashed lines in Fig. 10.
Here again, similar to Sect. 4.1, one can observe that the central deflection noticeably increased for the less restrictive boundary conditions. As a consequence, the results obtained with the models FOSD_2 and FOSD_3 (sharing a united graph in Fig. 10) were shifted away from the exact solution, whereas the results calculated with all other models moved closer to the exact solution. However, there is again just an occurrence of an artificial improvement, considering that the support conditions described by Pagano and Hatfield [40] corresponded to the case of the hard simply support (39), not the soft simply support (40). This example of a simply supported square (b = a) sandwich plate under a transverse load with bi-harmonic distribution (38) was taken from the paper of Pagano [33] who provided also the 3D elasticity solution for that problem. For comparison purposes, the results acquired from the graphs in [33] were additionally verified with digital solution provided by Loredo and Castel [41]. The plate consists of a soft transversely isotropic core (thickness 0.8h) bonded between two thin face sheets (thickness 0.1h each) made of graphite/epoxy composite with the same material parameters as considered in Sect. 4.1. The parameters of both materials are given in Table 2.
The comparative study of the central deflection of the square sandwich plate under bi-harmonic transverse load was performed first for the hard simply support boundary conditions (39). The results obtained with the examined models and related to the exact 3D elastic solution are presented in Fig. 11 for varying length-tothickness ratio.
From the graphs in Fig. 11, it can be seen that the FOSD models without shear correction (FOSD_0) or even with the local shear correction factor κ = 5/6 (FOSD_1) failed to provide an acceptable estimation of the central deflection; even for the thinner plate (a = 100h) their error reached almost 1%. The FOSD models with a global shear correction performed much better, without any visible difference between the results obtained with the models FOSD_2 and FOSD_3, which slightly contradicts observations reported earlier for sandwich plates in [42].
The results of the all three zigzag models are very close to the 3D elastic solution. Just to allow for a more detailed analysis, the graphs from Fig. 11 are repeated in Fig. 12, but with the error range limited to < −1%; +3% >.
It can be observed in Fig. 12 that the model ZZ_2 offers the best approximation of the exact 3D solution, whereas both models based on RMVT (series ZZ_1 and ZZ_3 symbolized in Fig. 12 by a common representation) gave slightly worse results, but being still within a very narrow error band. On the other hand, the  observer can be somewhat confused by the fact that the ZZ_3 results show greater errors than those of ZZ_2. It is also interesting to observe that for a/ h > 10 FOSD_2 and FOSD_3 offer an accuracy which is comparable to that obtained with zigzag models. From the discussion of Sect. 4.1, it is obvious that the hard simply support conditions (39) are compatible with the boundary conditions specified by Pagano [33]; however, one could be still interested how the results of the analyzed example would change with the application of the soft support (40); that issue is illustrated in Fig. 13.
Generally, with less restrictive boundary conditions the central deflection of the plate increased; on closer inspection, one can observe that the error of the "soft support" results for FOSD_2, FOSD_3 and all the zigzag models converges to the value near 0.25% for thin plates (a = 100h). With this latter conclusion, we terminate examination of the "soft" simply support variants for all the following examples. This is the first from three benchmark examples of simply supported square sandwich plates under bi-harmonic load recommended by Carrera and Brischetto [43]. The "softness" of the core (thickness 0.8h) with respect to the stiffness of the faces (thickness 0.1h each) is even more pronounced than in the previous example-compare Table 3.
The central deflection has been normalized following [43] as The results obtained with the FOSD models for the hard simply support boundary conditions (39) are presented in Table 4 and compared with the reference 3D solution [43]. The corresponding results for the analyses performed with zigzag models are given in Table 5.
It is worth to highlight that the 3D solution provided by Carrera and Brischetto [43] took into consideration the transverse normal deformations which can be very essential for sandwich plates with soft cores. The  computational models considered in the current report neglected that effect. In this context, no one should be surprised to observe that the FOSD models can provide a reasonable estimation of the central deflection only for very thin plates. However, one can notice here a visible advantage of the FOSD models with the global shear correction (FOSD_2 and FOSD_3) giving an excellent solution also for thin plates (a = 100h). As a great surprise, one can consider a superb performance of the zigzag models with just 1% error for thick plates (a = 4h). It is also interesting to see that model ZZ_1 which does not account on the layer stiffness differences provided almost the same results as those obtained with models ZZ_2 and ZZ_3.

4.6
Simply supported 3-layer sandwich square plate with extremely soft core The second benchmark example considered by Carrera and Brischetto [43] was a simply supported square sandwich plate with artificially softened core characterized by the stiffness parameters of the Nomex core from the previous example (Table 3) but reduced by a factor of 100. With the material parameters of the faces equal to those shown in Table 3, the disproportion between Young's moduli of the faces and the core reached an extremely high level (7.3 × 10 6 ). The core thickness was 0.8h, and each of the faces had the thickness 0.1h. The results for the central deflection of the plate under bi-harmonic load (38) for the hard simply support boundary conditions (39) normalized according to (41) are compared with the reference 3D solution [43] in Table 6 (FOSD models) and Table 7 (zigzag models).
With such a high level of the stiffness discrepancy between layers as assumed in the considered example, the FOSD models are practically unable to provide a reasonable solution unless the plate is extremely thin (a = 1000h). What really strikes the eye in Table 6 are incredibly high errors of the models FOSD_2 and FOSD_3 for the thicker plates (a ≤ 10h). It seems to demonstrate the limits of the idea to determine the global shear correction by matching the transverse shear strain energy predicted by the FOSD plate model with that obtained from the 3D elasticity theory. On the other hand, as a huge surprise one can consider an excellent performance of the zigzag models as displayed in Table 7. A simply supported square sandwich plate with Nomex core (thickness 0.8h) and composite faces, each composed as two-layer cross-ply composite laminate (2 × 0.05h), was suggested by Carrera and Brischetto [43] as their third benchmark example. The material parameters used in the computations are displayed in Table 8.
The central deflection of the 5-layer [0/90/c/90/0] simply supported square plate under bi-harmonic load (38) calculated with examined computational models normalized according to (41) and compared with the reference 3D solution of Carrera and Brischetto [43] are displayed in Table 9 (FOSD models) and Table 10 (zigzag models).
As one can observe in Table 9, the FOSD models performed significantly better for this example than for the previous two. However, models FOSD_0 and FOSD_1 offered an acceptable accuracy only for an extremely thin plate (a = 1000h), but models FOSD_2 and FOSD_3 behaved much better providing 2% error for a moderately thick plate (a = 10h) and 13% error for the thickest of the analyzed plates (a = 4h).
At the first glimpse, the performance of the ZZ_1 model (cf . Table 10) can be considered as a profound disappointment as its accuracy is not better than that of the model FOSD_0, whereas the models ZZ_2 and ZZ_3 provide excellent results even for the thick plate (error 0.8% for a = 4h). After careful consideration, however, one can discover that this evident advantage of the models ZZ_2 and ZZ_3 results directly from the application of the zigzag function proposed by Tessler (12) with variable slope dependent on the variation of the layers' stiffness. Apparently, the essential difference between the slope change at the first and at the second interlayer boundaries available in the models ZZ_2 and ZZ_3 allows for a much more accurate representation of the through-the-thickness deformation profile. On the other hand, the zigzag function of Murakami (10) does not depend on the material parameters; therefore, the zigzag part shown in (9) practically is not active for the particular considered example of [0/90/c/90/0] sandwich plate which results in the performance corresponding to that of a pure FOSD model.

Simply supported rectangular composite [90/0] plate
The next example originates from a popular benchmark problem of a square composite [90/0] plate introduced in [33]; however, Pagano [33] presented only graphs without any numerical results for his example; on the other  [39] hand, the numerical results available for that problem in the literature (see e.g., [44]) are slightly questionable. Seeing that Zenkour [39] presented tabular results for an analogous two-layered plate with the only difference in the rectangular shape, the latter case has been selected to be analyzed as the current example. The comparative study of the central deflection of a rectangular (b = 3a) two-layered [90/0] plate under bi-harmonic transverse load (38) was performed for the hard simply support boundary conditions (39). The analyzed plate consisted of two layers (thickness 0.5h each) made of the same high modulus graphite/epoxy composite material as considered in Sects. 4.1-4.3. The results obtained with the examined models and related to the reference 3D solution of Zenkour [39] are presented in Fig. 14 for a varying value of the length-to-thickness ratio.
In contrast to what was observed in the previous examples, the best performance of the all FOSD-based formulations is shown by the model FOSD_1, while the deflections obtained with the models FOSD_2 and FOSD_3 are noticeably bigger. It can also be interesting that the model ZZ_2 performed just slightly better than the model FOSD_0 which possesses no shear correction. The best estimation of the central deflection of the plate is provided by the zigzag models based on the RMVT series ZZ_1 and ZZ_3 which share the same line in Fig. 14.

Simply supported antisymmetric composite [−15/+15] square plate
This example proposed by Savoia and Reddy [45] was considered afterward by many researches including Kulikov and Plotinkova [46] or Loredo and Castel [41]. A two-layered simply supported square plate is constructed of two composite layers with the same thickness but with a different orientation of the fiber reinforcement. With the antisymmetric lamination scheme [−15/+15], the considered problem would not satisfy the stipulations of a double symmetry, despite the double symmetry of the geometry, the loading and the boundary conditions. Consequently, in contrast to Sects. 4.1-4.6, the whole plate must be discretized in the FEA. Another crucial difference with respect to all the previous examples consists in a distinctive specification of the simply support applied by Savoia and Reddy [44]. As a result, the boundary conditions used in the current FEA calculations were given by the following specifications: that match neither formula (39) nor (40), but fully harmonize with the set "SS-2" of simply support conditions defined in [2]. It has to be stressed that the zigzag function proposed by Tessler (12) for the antisymmetric lamination schemes [−θ/+θ ], [−θ/+θ/−θ/+θ ] and other of that kind characterized by the same absolute value of the fiber orientation angle in subsequent layers tends to zero which finally results in a singularity of the stiffness matrix. For that reason, a special approach should be implemented to ensure a solution for the computations performed with the models ZZ_2 and ZZ_3. In fact, the authors of the current report tested two different tactics. In the first approach (with corresponding results denoted as ZZ_2A and ZZ_3A), the lamination scheme [−15/+15] was replaced simply with the sequence [−14.95/ + 15.05], which corresponds to an introduction of a small imperfection in trajectories of the fiber reinforcement. The second tested remedy consisted in an appropriate modification of the zigzag function, which for the antisymmetric lamination schemes [−θ/+θ ] reflected as a matter of fact the shape of the zigzag profile of Murakami. The results calculated with that second tactic will be presented as ZZ_2B and ZZ_3B.
The comparative study performed for the central deflection of the square composite plate with antisymmetric layup [−15/+15] under bi-harmonic transverse load is illustrated in Fig. 15. The results obtained with the examined models and related to the exact 3D elastic solution of Loredo and Castel [41] are presented in the form of graphs for a varying value of the length-to-thickness ratio.
It is quite interesting that the global transverse shear correction factors calculated for the FOSD_2 model according to the procedure suggested by Whitney [19] are equal to 5/6, which exactly matches the value of the classical shear correction factor κ = 5/6 introduced by Reissner [15]. As a consequence, the FOSD_2 model yields exactly the same response as the FOSD_1 model. Better results may be obtained with the model FOSD_3 where the procedure for the evaluation of the global composite transverse shear stiffness following Rolfes and Rohwer [23] involves an independent calculation of the off-diagonal termsã 45 in (30), which seems to be crucial for angle-ply laminates. In fact, this is the first example in our report where differences between the responses obtained with the models FOSD_2 and FOSD_3 are so visible. A "pure" FOSD formulation applied in the FOSD_0 model allowed for obtaining a reasonable solution only for a very thin plate (a = 100h).  Quite surprising might be the somewhat disappointing performance of the zigzag models. Although the altered lamination scheme [−14.95/+15.05] in the models ZZ_2A and ZZ_3A allowed for avoiding the singularity of the stiffness matrix, the obtained results did not converge to the correct solution even for very thin plates. As can be seen from the graph, the second alternative implemented in the models ZZ_2B and ZZ_3B appears to be more effective. The series ZZ_2B moved closer to the reference solution, while the line representing ZZ_3B approached the exact results of the model ZZ_1. However, it must be considered as a big surprise that the model ZZ_1 (and ZZ_3B) could not provide a better accuracy for the current problem than the model FOSD_1. A more detailed examination would show that the series ZZ_1, ZZ_2B and ZZ_3B converge to the error 0.3% (for a = 100h), whereas the error of the series FOSD_1, FOSD_2 and FOSD_3 approaches a level < 0.01%. This reflection seem to be too important to be left without at least an attempt of an explanation. To explore that issue more deeply, additional series of computations have been performed with a little modification of the boundary conditions where the set of constraints (42) was supplemented with the following conditions: The results of computations for those new boundary conditions are marked with a single star (*) in the graph presented in Fig. 16.
As one can observe in Fig. 16, the analyzed model with modified boundary conditions became more stiff, which resulted in a spectacular "improvement" in the accuracy of the models ZZ_1*, ZZ_3A*, ZZ_3B* and also ZZ_2A* as compared to the results obtained earlier with the models ZZ_1, ZZ_3A, ZZ_3B and ZZ_2A. On the first glimpse, the correctness of the model ZZ_2B* was noticeably aggravated with respect to that of the model ZZ_2B. However, after a more detailed inspection, one can find that the errors of the series ZZ_2B*and ZZ_2A* are represented now by collinear graphs which converge approximately to zero (actually less than 0.04%-as can be seen in Table 11).
It is also quite symptomatic that there is almost no difference between solutions obtained with the models ZZ_1* (= ZZ_3B*) and ZZ_3A* (cf. Table 12).
At a first glance, this coherence is somewhat surprising; however, it seems rather not to be just a pure coincidence. More inquiry would be definitely required, although a collection of 3D elastic solution results available in the open literature for composite plates with antisymmetric angle-ply layup is extremely limited. One should also be aware that the matter of a proper description of the boundary conditions for a composite  plate with antisymmetric layup is not a trivial issue. The 3D elastic solutions presented in [41] for composite plates with antisymmetric layups were obtained with a special set of boundary conditions (for the details see [47]) that actually did not match the set "SS-2" of simply support conditions defined by Reddy [2]. On the other hand, most practitioners would be more concerned with the representation of the real-world physical conditions in their computational models rather than trying to fulfill somewhat elaborated constraints resulted from the limits of available analytical solutions.

Conclusions
The performance of various ESL models in a linear static finite element analysis of multilayered plates has been investigated in the paper. The authors have considered altogether seven computational models gathered in two collections: the FOSD Set consisting of four variants of the classical first-order shear deformation formulation diversified by an application of different transverse shear correction, and the ZZ Set containing three refined models based on the application of the zigzag functions. The comparative examination of the selected models has been performed for the assortment of well-established benchmark examples with available exact 3D elastic solutions for static bending problems of simply supported plates under a transverse load with bisinusoidal distribution. A special attention has been paid to more detailed analyses concerned with examination of differences related to the "soft" and "hard" variants of the simply support. As one could expect, the computational models enhanced by the application of the zigzag functions (gathered in the ZZ Set) generally performed better than the FOSD formulations. However, they seem to be too stiff for the analysis of laminates with larger number of layers. Such multilayered laminates are better approximated by FOSD models with global shear stiffness correction. From all the analyzed ZZ models, a better performance was observed in the case of approaches constructed based on the RMVT (models ZZ_1 and ZZ_3) compared to the model ZZ_2 constructed with PVD.
Generally, the models based on a pure FOSD formulation (FOSD_0) or eventually modified by the application of a local transverse shear correction (FOSD_1) can provide reasonable results only for thin, or even, extremely thin plates. Strategies for an appropriate correction of the global transverse shear stiffness incorporated in models FOSD_2 and FOSD_3 shared a common background consisting in a matching of the transverse shear strain energy predicted by the plate model with that established by the three-dimensional elasticity theory. As a consequence, the results provided by both those models were almost identical for most of the considered examples. The only differences were encountered for the plate with the antisymmetric lamination scheme [−15/ + 15], where the strategy proposed by Rolfes and Rohwer [23] could fully reveal its superiority resulting from a more comprehensive approach. On the basis of obtained results, one can conclude that the proper estimation of a global (cross-sectional) transverse shear correction is essential for achieving an acceptable estimation of the deflection of moderately thick multilayered plates with FOSD formulations. One can expect that a proper value of the shear correction factor would be even more significant in a nonlinear analysis as indicated quite recently by Chróścielewski et al. [48] for the six-parameter finite rotation analysis of elastic shells (see also [49] for examples of nonlinear analysis of composite shells under thermal loads).
At the very final stage of preparing of this report, the authors came across two interesting papers [50,51] related to the application of zigzag models in the analysis of multilayered plates. In the first one [50], we can find an alternative RMVT-based implementation of Tessler's zigzag functions. The second of those publications [51] contains a stimulating assessment review of various zigzag formulations. The focus of both of these publications is aimed on the stress recovery problems being, however, beyond the scope of the present report, but planned to be reported in the subsequent publication [52].