1. 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 C1-continuous as in homogeneous non-layered materials, but instead they should be C0-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
,
and transverse shear rotations
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.
2. Preliminaries and Notations
Constituent layers are assumed to have a uniform arbitrary thickness 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 () is assumed as the reference frame, having () on the middle reference plane of the laminated plate (origin in the lower left edge) and as the thickness coordinate (, with being the overall thickness). Lx and Ly symbolise the plate side-length in the x and y-directions, while and 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 , while u and l mark the upper and lower faces of the laminate, and a comma is used to indicate spatial derivatives, e.g., , . The elastic in-plane and transverse displacement components are indicated as and , respectively. Strains are assumed to be infinitesimal and, to distinguish their origin, they are specified as , , respectively when they come from kinematic or stress-strain relations. Once assumed as primary variables, they are indicated as , , (; , ; ; ). Stresses from stress–strain constitutive relations are indicated as , while when they are assumed as primary variables, they are indicated as , , .
2.1. Recalls on Mixed Variational Theorems
Here a generalized version of Hu-Washizu theorem is used, whose primary displacement boundary condition link is weakened as
and which generates the following variational statement:
Here represents the canonical functional, are the components of the external unit normal to the volume bounding surface, are the components of body forces (they will contain inertial forces, as explained in (4)) and , 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 of the laminated plates is assumed to be bounded by a surface that is split into a surface on which surface tractions are prescribed and a surface on which surface displacements are prescribed. Body forces on , prescribed surface tractions on and prescribed displacements on are assumed to act.
Theories with displacements and out of plane stresses assumed separately are developed from the HR variational theorem:
where
is the first variation of the HR canonical functional, and the following assumptions are made:
;
, being
;
.
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.
2.2. 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
and trial functions
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
,
,
,
,
, 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
, 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 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
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
, while the condition in Equation (8) is reformulated as:
At the simply supported edges, the following boundary conditions are enforced
on the reference mid-plane of plates at
,
, and
,
. 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.
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 , 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.
4. 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.
4.1. 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 , , , , and are still calculated by enforcing the fulfilment of the stress boundaries conditions in Equation (16), while 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 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 . The right slope is determined at each interface without a cost burden evaluating what sign attain the lowest residual force norm from three local equilibrium equations.
4.2. MHR4, MHWZZA and MHWZZA4 Theories
A refined variant MHR4 of MHR is obtained assuming the in-plane displacement field
by Equation (31), and a fourth-order piecewise variation of the transverse displacement [
47]:
Coefficients to are determined by enforcing the fulfillment of the stress boundary conditions in Equation (16), whereas 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.
4.3. 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 is assumed to be the same as that of the ZZA model, while the transverse shear stresses 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, , is the inverse of ; ; ; ; , and .
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,
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
are obtained through integrating local equilibrium equations, while the
appearing in Equation (2) is assumed to be the same as that of ZZA.
5. 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
Table 3 and
Table 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
Table 6,
Table 7,
Table 8,
Table 9,
Table 10,
Table 11,
Table 12,
Table 13,
Table 14,
Table 15,
Table 16,
Table 17,
Table 18,
Table 19 and
Table 20 and the computational cost is reported in
Table 21 and
Table 22.
5.1. 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
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 =
Lx; 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
; 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
Lx/
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 through-thickness 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 (
Lx/
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
never reverses, while
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
Section 3.1 and
Section 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
Table 21 and
Table 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.
5.2. 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-to-thickness ratio.
Modal displacements and stresses are reported In
Figure 3,
Figure 4 and
Figure 5, where the results are normalized as shown in
Table 5, and show that the axial displacement
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
, 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
Figure 3 and
Figure 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.
5.3. 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
), 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
, transverse shear
, and transverse normal
modal stresses are reported in
Figure 6 and
Figure 7 for ratios E
1/E
2 of 3 and 40, respectively. Note that
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-to-thickness 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 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.
5.4. 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
Table 10 and
Table 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.
5.5. 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.
Table 12,
Table 13 and
Table 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 Lx/h = 20, all of the theories except FSDT (shear correction factor ) 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.
The first eight frequencies for a [0/90/0] square plate with clamped edges (CCCC), which is retaken from Zhen and Wanji [
4] and Liew [
52], and is here referred to as case e2, are reported in
Table 15. In this case, in addition to 3D FEA, comparisons can be carried out with results obtained 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.
5.6. Through-Thickness Mode of a Simply Supported, Cross-Ply Plate in Cylindrical Bending
A simply supported cylindrically bent [0/90] plate (case f) is now studied, which is retaken from Pagani et al. [
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.
5.7. First Five Free Vibration Frequencies of a Soft-Core Sandwich with Strong Transverse Normal and Shear Deformability Effects
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.
5.8. Free Vibration of a Thick Simply Supported [0/90/Core/0/90] Sandwich Plate
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.
5.9. 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
tp 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±, whose 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
Figure 11 and
Figure 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 ) 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.
5.10. Computational Effort of Theories
Table 21 and
Table 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 result in slower processing times than the adaptive theories for static cases whenever stresses must be computed through time-consuming procedures, but if this is not required, that is only global quantities are required, they result in faster processing times than the adaptive theories. However, this advantage nullified the results of HRZZ and HRZZ4, which were inaccurate in almost all of the cases that were considered.
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.
6. 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 always achieve the best accuracy with a processing time that is still short, while the other theories have lower calculation times, but also a much lower accuracy.