A higher-order plate theory for the analysis of vibrations of thick orthotropic laminates

In this contribution, flexural vibrations of linear elastic laminates composed of thick orthotropic layers, such as cross-laminated timber, are addressed. For efficient computation, an equivalent single-layer plate theory with eight kinematic degrees of freedom is derived, both in terms of equations of motion at the continuum level and in terms of a finite element representation. The validity of the plate theory is demonstrated by comparing natural frequencies of a simply supported plate over a wide frequency range for which an analytical solution is available. Furthermore, the influence of individual material properties on the accuracy of the plate theory is investigated, demonstrating its broad applicability. The influence of material orientation on the accuracy of the plate theory is examined by comparative finite element simulations. It is shown that, for cross-ply laminates, the plate theory is valid for elements oriented at any angle to the material principal axes.

zig-zag theory approaches from literature: (i) The Ambartsumian Multilayered Theory [10] as extension of the classical Reissner-Mindlin theory that applies to homogeneous plates, (ii) theories based on the Reissner mixed variational theorem [11], and (iii) the Lekhnitskii Multilayered Theory [12]. All theories are based on axiomatic assumptions regarding the transverse shear stress field and/or the displacement field. In an effort to unify the plate theories in one framework, the Carrera Unified Formulation [13] was derived.
In the present paper, an approach based on the work of Lekhnitskii [12], originally developed for composite beams, is followed. Lekhnitskii's theory was extended for the analysis of plates in [14] and further refined in [15], specifically for cross-laminated timber (CLT). In this type of composite, timber panels are stacked perpendicular to each other, corresponding to a cross-ply laminate. In its linear elastic range, timber can be characterized as orthotropic material with principal axes parallel to the annual rings (longitudinal axis L) and radial (R) and tangential (T ) to the annual rings, respectively, in a cylindrical coordinate system [16]. Since the material properties in the R-and T -direction are of similar order, no distinction is made between the two in practice. Then, a Cartesian coordinate system can be used to describe the constitutive behavior, with timber represented as a transversely isotropic material, except that shear moduli and Young's moduli are decoupled. More specifically, the shear moduli G L R ≈ G LT and G RT differ by a factor of approximately 15, and the ratio of Young's modulus to shear modulus E L /G L R ≈ 17. Thus, CLT can be characterized as a highly shear-flexible laminate in which the aforementioned zig-zag deformation pattern prevails due to the abruptly changing shear moduli in the layers. The plate theory developed in [15] was extended in [17] for dynamic loading, focusing on the sound-radiating properties of CLT panels. Generally, due to the typically high stiffness-to-mass ratio of composites, vibrations and as a result, sound radiation, can become more relevant than for isotropic plates. In [17], it was shown for a simply supported plate that the proposed plate theory can more accurately capture natural frequencies in a frequency range of 20 to 550 Hz than with FSDT. However, with increasing frequency, i.e., with increasing cross-sectional curvature of the mode shapes, the discrepancies to the 3D elasticity solutions increased. Therefore, in [18], on the one hand, the plate theory was extended to capture the out-of-plane deformation of the plate, which has been neglected until then. On the other hand, insignificant terms were omitted and the kinematic assumptions were reformulated in view of a straight-forward finite element (FE) implementation, avoiding derivatives of kinematic degrees of freedom in the displacement field. The strong and weak formulations, the latter in terms of finite element matrices and load vector, were derived, and both the accuracy and efficiency of the computations were demonstrated for a realistic CLT floor construction. One limitation of the plate theory presented in [18] is that it is only applies to symmetric lay-ups. Furthermore, it was assumed that the material principal axes of the layers correspond to global coordinate axes. Although this is no limitation as long as rectangular CLT panels discretized by regular FE meshes are considered, the plate theory is extended in the present paper for more general setups, i.e., angle-ply composites with non-symmetric lay-ups. The paper is organized as follows: in Sect. 2, the strong and weak formulation of the plate theory is derived. This is followed by numerical simulations in Sect. 3 in which, first, the validity of the plate theory is proved and, second, the range of applicability in terms of material properties and layer setups is investigated. The paper ends with conclusions in Sect. 4.

Kinematics and constitutive behavior
Laminated plates consisting of an arbitrary number n of plies whose mid-planes coincide with a global x-yplane and normal direction z are considered. For an orthotropic ply k of the laminate with principal axes 1, 2, and 3 oriented at an arbitrary angle with respect to the global x yz coordinate system (see Fig. 1), the stresses σ k and strains ε k in the ply are connected via the constitutive equations where the stiffness coefficients Q k i j are related to the material properties E 1 , E 2 , E 3 (Young's moduli), G 12 , G 13 , G 23 (shear moduli) and ν 12 , ν 13 , ν 23 (Poisson ratios) [2]. When the principal material axes are aligned with the global coordinate system, Q k 16 , Q k 26 , Q k 36 and Q k 45 are zero. Following the work presented in [18], the following functions are considered for the displacements in the x-, y-an z-direction in layer k, denoted u k , v k and w k : In contrast to [18], in-plane displacements u 0 and v 0 are taken into account, so that the plate theory is also valid for non-symmetric lay-ups. The further degrees of freedom of the plate theory are the out-of-plane displacement w 0 and the correction term w to capture thickness changes via the function f (z) and the rotations ϕ x , ϕ y , ψ x and ψ y . The functions A k (z) and E k (z) capture the zig-zag deformation pattern across the plate thickness and are defined as with the functions The integration constants C k a and C k e in Eq. (3) can be determined by application of compatibility conditions at the layer interfaces h i , i = 1 . . . n + 1, and the additional constraint that at z = 0 the corresponding function (e.g., A k (z)) is zero [15]. The scaling constantQ with Q m as the mean shear stiffness of the cross section and the plate thickness h is introduced in Eq. (4) to render functions a k (z) and e k (z) in units of stress and functions A k (z) and E k (z) in units of a length. Note that the functions in Eqs. (3) and (4) are polynomials of order 3 and 2, respectively, so integration along the cross section, as introduced below, is a straightforward task.

Equations of motion
The equations of motions in terms of the eight kinematic degrees of freedom are derived by application of d'Alembert's principle [19] in the form (σ x x δε x x + σ yy δε yy + σ zz δε zz + 2σ xy δε xy + 2σ yz δε yz + 2σ xz δε xz − p z δw where the superscript k for stresses, strains and displacements has been omitted for the sake of readability. ρ denotes the mass density, S the mid-plane area and p z (x, y, t) the transverse loading of the plate. For small strains, w} and x i = {x, y, z}. Substitution of the displacement field Eq. (2) yields the strains and, taking into account Eq. (1), the stresses as a function of the kinematic degrees of freedom and, similarly virtual strains in terms of their variations. Partial integration eliminates the derivations of the variations. Performing integration over the cross section h and setting the variations to zero and one mutually, we obtain after some manipulations a system of eight coupled partial differential equations of the form where q(x, y, t) = u 0 v 0 w 0 w ϕ x ϕ y ψ x ψ y T is the vector of degrees of freedom, and derivations of (.) with respect to x, y and t are denoted by (.) ,x , (.) ,y and( .), respectively. The matrices in Eqs. (6) read . . . . . . . . . . .
while the load vector is For better illustration, zeros have been replaced by a dot (.) in the above matrices. The plate stiffness parameters apparent in Eqs. (7)-(13) are defined as and the corresponding plate mass parameters read Since, as mentioned earlier, all functions in the thickness direction z are layer-wise polynomials, the plate stiffness and mass parameters can be computed efficiently, for instance, using the polyint and polyval functions in MATLAB.

Dispersion relation
Assuming planar bending waves with wave vector k(x, y, z) = k x k y 0 T and angular frequency ω in the form (6), after some manipulations, yields the eigenvalue problem K d − ω 2 M A = 0, where M is given in Eq. (7) and elements of the complex conjugate stiffness matrix in this dispersion relation K d read The dispersion relation can be used to compute the natural frequencies of bounded plates if the wave number components k x and k y for the corresponding standing waves are known. For instance, mode shapes of a rectangular, simply supported plate with dimensions l x and l y in x-and y-direction, respectively, are characterized by k x = mπ/l x and k y = nπ/l y with integers m, n = 1, 2, . . .. Solving K d − ω 2 M = 0 for ω and considering the lowest of the eight eigenvalues, we obtain the natural frequency. Solving for the corresponding complex eigenvector A, we obtain the amplitudes for the kinematic degrees of freedom q for this mode.

Weak formulation of the governing equations
In the finite element implementation of the plate theory, the strains in the k-th layer are expressed as and With the constitutive equations σ = Qε according to Eqs. (1), the work of the internal forces in d'Alembert's principle can be written as where with the submatrices The corresponding mass property matrix M f e and load vector f f e are the same as defined in Eqs. (7) and (14). Degrees of freedom and geometry are interpolated by isoparametric shape functions [20] in the form q = N q e where q e = q (1) q (2) . . . q (r ) T are the nodal degrees of freedom of the r -noded finite elements and N is an 8 × 8r shape function matrix. Denoting the derivatives of the shape functions B = LN, inserting these approximations into d'Alembert's principle yields for integration over one finite element the elemental mass and stiffness matrix and the elemental load vector according to where J is the Jacobian relating global and local coordinates [20]. Integration is performed by Gauss quadrature in the local ξ, η coordinate system.

Computations
In the applications presented here, the main focus is on cross-laminated timber. A rectangular CLT panel with dimensions l x × l y = 4 × 5 m and an antisymmetric lay-up with layer thicknesses 2 − 4 − 2 − 4 cm is considered. Unless otherwise specified, orthotropic material properties E 1 = 12.000 N/mm 2 , E 2 = E 3 = 400 N/mm 2 , G 12 = G 13 = 800 N/mm 2 , G 23 = 50 N/mm 2 and ν 12 = ν 13 = ν 23 = 0.3 as well as a density of ρ = 450 kg/m 3 are assigned to all layers. As correction function f (z) for the out-of-plane displacement, a quadratic parabola centered in mid-plane z = 0 is used, i.e f (z) = 4/ h 2 z 2 . For simply supported laminates consisting of thick orthotropic layers, analytic solutions based on 3D elasticity theory are provided in [21]. Specifically, for free vibration the displacement field is written as Computing the relative errors f between the natural angular frequencies ω mn based on the equations in [21] and their counterparts based on the proposed plate theory, as described in Sect. 2.3, yields for the first n = 1 . . . 6 × m = 1 . . . 6 modes the result shown in Fig. 2a. While for lower modes the error is practically negligible, with increasing mode numbers m/n also f increases. However, for all 36 computed natural frequencies the relative error is less than 1 %. Alternatively, if one plots the relative errors versus the corresponding natural frequency f = ω/2π, one obtains the graph shown in Fig. 2b denoted as "8 DOF". The increase in relative error with increasing frequency can be clearly seen. In this figure, two more results are depicted: First, the relative errors for the proposed plate theory when degree of freedom w is neglected, denoted as "7 DOF". Over the entire frequency range, the relative error is approximately constant at about 1%. It can be concluded that taking w into account in the equations of motions yields more accurate results for the natural frequencies at low frequencies, but also that the employed function for f (z) seems to become more inappropriate with increasing frequency, i.e., increasing cross section rotations and curvatures. Second, in Fig. 2b results for the relative error of the natural frequencies neglecting degrees of freedom u 0 and v 0 are depicted ("6 DOF"), which follow from the plate theory proposed in [18] and apply for symmetric cross sections. In this case, the errors are an order of magnitude larger (note the logarithmic scale in Fig. 2b), which highlights the importance of the two in-plane displacements for nonsymmetric cross sections, even when the boundary conditions are free of in-plane constraints, as in the present problem. Figure 3 provides a more detailed comparison of the results based on 3D elasticity and the plate theory (denoted here as 2D). For two modes, the cross section functions (z), ψ(z) and χ(z) in Eq. (26) are shown, which govern the displacement in x-, y-and z-direction along the cross section. For the mode with m = 2, n = 3 (Fig. 3a), excellent agreement is obtained for all three functions and the natural frequency by the proposed plate theory, employing the above symmetric quadratic parabola for f (z). On the other hand, at a higher natural frequency corresponding to m = 6, n = 4, the zig-zag deformation pattern in (z) and ψ(z) is more pronounced due to the larger cross section rotation and curvature, as visible in Fig. 3b. χ(z) tends to become more asymmetric; therefore, f (z) approximates less accurately the actual out-of-plane deformation. However, the overall agreement between 3D elasticity and the proposed plate theory is still very good.
The weak formulation of the plate theory presented in Sect. 2.4 was implemented in MATLAB. Discretizing the plate with isoparametric finite elements of size of 0.2 × 0.2 m and eight nodes each, taking into account  Relative error for the natural frequencies f of 36 modes of a simply supported plate based on the FE implementation of the plate theory compared to 3D elasticity quadratic shape functions, yields an FE model with 12,728 degrees of freedom. Performing an eigenvalue analysis for the first 50 natural frequencies takes 6.7 s on a personal computer with six Intel Xeon E5-1650 v4 processors. In Fig. 4, the relative errors f in the computed natural frequencies are depicted in comparison with the analytical outcomes according to [21]. The modes from the numerical simulation were assigned to mode numbers n and m by means of the modal assurance criterion (MAC), defined for the two real mode shape vectors ϕ 1 and ϕ 2 as [22] MAC(ϕ 1 , Specifically, the component q 3 = w(x, y) from the nodal degrees of freedom was compared with the sinusoidal modes shapes given in Eq. (26) to find n and m accordingly. The deviation of natural frequencies based on the FE implementation is less than 1 % for all 36 modes, i.e., it is within the accuracy of the plate theory itself (compare with Fig. 2a). This proves the validity of the FE implementation.
To investigate the applicability of the proposed plate theory for a wider range of materials and/or properties, these parameters are varied. To limit the possible variations, the structure of the material properties, i.e., E 2 = E 3 , G 12 = G 13 is not changed in the following parametric study. As in the previous section, the results for the first n × m = 36 natural frequencies of a simply supported plate with dimensions described above are compared. The outcomes of the analytically available solution [21] are set in contrast to those of the (strong formulation of the) plate theory and shown as the relative error f in Fig. 5. Increasing Young's modulus E 2 such that it approaches E 1 has a minor effect on the accuracy of the plate theory, as shown in Fig. 5a. Up to E 2 /E 1 = 1/2, the relative errors in the considered 36 natural frequencies have almost the same appearance. When G 12 is increased with respect to E 1 (Fig. 5b), the errors in the natural frequencies increase. Since G 23 is not changed, also in this case the ratio G 23 /G 12 decreases (for G 12 /E 1 = 1/5 this would correspond to G 23 /G 12 = 1/48). This is also investigated in the third parametric study, whose results are shown in Fig. 5c. Here G 23 /G 12 is changed, while G 12 /E 1 remains the same. From the latter two parametric studies, it can be concluded that as the difference of the shear moduli G 23 to G 12 increases, i.e., as the zig-zag deformation pattern becomes increasingly pronounced, the errors in the natural frequencies increase. However, only for relatively unrealistic conditions, when G 23 /G 12 is smaller than 1/30, the relative errors in the considered frequency range exceed 1%. Therefore, it can be concluded that the plate theory is valid for a wide range of material properties.

Influence of laminate orientation
In the analytical solution [21], a layer-wise orthotropic stress-strain relation is assumed with Q i6 = 0, i = 1, 2, 3 and Q 45 = 0 in Eq. (1). This means that local axes 1 and 2 are aligned parallel to the global axes x and y (compare with Fig. 1). This is also the underlying assumption of the plate theory originally proposed in [15]. Since the plate theory presented in this paper is formulated for the general case where an angle 0 ≤ α ≤ 90 • can occur between the two coordinate systems, the influence of this orientation on the accuracy of modal properties (natural frequencies and mode shapes) is compared here. Note that the considered plate is still a cross-ply laminate, but the dimensions and boundaries are no longer oriented parallel to the principal material axes. Since no analytical solution is available for this case, the results of FE simulations conducted in ABAQUS serve as a reference. The considered four-layered laminated plate with the above dimensions and material properties is discretized with 20-noded hexahedral elements with quadratic shape functions. A very fine mesh with 24 elements over the thickness and in-plane dimensions of 0.1 × 0.1 m is used, yielding a FE model with more than 600.000 degrees of freedom. To minimize the effect of the implementation of the boundary condition, the plate is not considered simply supported in this study, but clamped at all edges, i.e., all available degrees of freedom at the edges are set equal to zero. A modal analysis is performed for the first 40 modes and compared with the outcomes of computations based on the weak formulation of the plate theory in Fig. 6. Mode shapes are compared using MAC values evaluated at 60 nodal positions of the plate on the left hand side of the figure, where "3D FEM" refers to the ABAQUS reference solution and "2D FEM" to the solution obtained by the plate theory. On the right hand side, the corresponding natural frequencies between the two simulations are compared by means of the relative errors f .
When the material axes are aligned in parallel to the global axes (α = 0 • ) as before, Fig. 6a shows a very close agreement of the modal properties: The MAC values are aligned along the diagonal with values approaching 1, which means that the two sets of modal vectors are fully correlated. The relative error of the natural frequency increases from values slightly below 0-1% in the frequency range up to 300 Hz. This behavior is consistent with the previous studies; compare, e.g., Fig. 5. On the other hand, if the laminate orientation angle α is increased to 20 and 40 • , respectively, the errors in the natural frequencies dramatically increase, as can be seen in Fig. 6b and c. The mode shapes, on the other hand, are not affected as much by the obviously incorrect plate theory, except for individual modes at higher frequencies. In Fig. 7 an explanation for this behavior is provided. The cross section functions A k (z) and E k (z) that govern the zig-zag deformation pattern (see Eq. (2)) are depicted for the considered plate for different angles α. Obviously, A k (z) for α = 0 • corresponds to E k (z) for α = 90 • and vice versa. For α between these two extreme values, A k (z) and E k (z) show a rather smooth deformation pattern without the characteristic zig-zag function. Since all the functions in Fig. 7 are commonly scaled to a minimal value of −1, the graphs illustrate that for α = 0 the magnitude of the two functions A k (z) and E k (z), i.e., its relative importance in the displacement field Eq. (2) depends strongly on α.
As a remedy, the nodal (x, y)-coordinates of the finite elements can be transformed into a coordinate system oriented in accordance with the principal material axes by means of a rotation matrix In this coordinate system, the constitutive relation for a cross-ply laminate (Eq. 1) simplifies, since Q i6 = 0, i = 1, 2, 3 and Q 45 = 0. Consequently, the entries in K f e in Eq. (20) vanish with these indices. Taking into account these local coordinate transformations in the FE simulations yields for α = 20 • and α = 40 • the results depicted in Fig. 8. In comparison with Fig. 6, the relative errors in the natural frequencies f are now of the same order of magnitude as for the case α = 0 • , and also the MAC matrix shows the targeted diagonal shape. Thus, it can be concluded for cross-ply laminates that rotating the nodal coordinates by the angle α before computing the element stiffness makes the FE implementation valid for arbitrary orientations of the element with respect to the principal material axes.

Conclusions
In this paper, an equivalent single-layer plate theory for thick laminated plates has been presented whose bending deformation is dominated by a zig-zag shape across the thickness, as is the case, for instance, for cross-laminated timber. Based on previous work by the authors, the plate theory is formulated for arbitrary cross sections (i.e., symmetric and asymmetric), considering orthotropic material behavior for the individual layers stacked at a general angle α. The deformation is governed by eight kinematic degrees of freedom, two in-plane displacements, two out-of-plane displacements, and four rotations. Equations of motion were derived, from which the dispersion relation governing the propagation of bending waves in unbounded plates was obtained. As special case, the natural frequencies and mode shapes of simply-supported rectangular plates were computed and contrasted with results from literature based on the solution of the 3D elasticity problem, revealing excellent agreement over a wide frequency range. The limiting factor in terms of accuracy was found to be the function for the out-of-plane deformation, which was assumed to be symmetric here. A parametric study with respect to material properties underscored the wide applicability of the plate theory in this sense. The weak formulation in terms of finite element elemental mass and stiffness matrices and load vector was derived, allowing the investigation of arbitrarily complex geometries. Varying the orientation of the composite plate, it was shown that the plate theory is valid for cross-ply laminates oriented at an arbitrary angle with respect to plate boundaries, if the elemental stiffness matrix is computed in a locally rotated coordinate system. On the other hand, the present study demonstrates that, as expected, the plate theory yields inaccurate results for angle-ply laminates. For these laminates, the assumed displacement field would need be adjusted. As a further extension of the plate theory, nonlinear material behavior could be implemented, e.g., by a damage formulation. Since even in symmetric composite lay-ups the failure of individual plies generally leads to asymmetric stress distributions, the proposed plate theory can be considered as important first step in this direction.
Funding Open access funding provided by University of Innsbruck and Medical University of Innsbruck.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.