# A FEM Free Vibration Analysis of Variable Stiffness Composite Plates through Hierarchical Modeling

Gaetano Giunta, Domenico Andrea Iannotta, Marco Montemurro

<figure><img src="image_1.png" /></figure>

## To cite this version:

Gaetano Giunta, Domenico Andrea Iannotta, Marco Montemurro. A FEM Free Vibration Analysis of Variable Stiffness Composite Plates through Hierarchical Modeling. Materials, 2023, 16 (13), pp.4643. 10.3390/ma16134643. hal-04214571

HAL Id: hal-04214571

https://hal.inrae.fr/hal-04214571

Submitted on 22 Sep 2023

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Article

# A FEM Free Vibration Analysis of Variable Stiffness Composite Plates through Hierarchical Modeling

Gaetano Giunta <fnref n="1" />,* Domenico Andrea Iannotta <fnref n="1" />,<fnref n="2" /> D and Marco Montemurro <fnref n="3" /> D

<figure><img src="image_3.png" /></figure>

check for
updates

**Citation:** Giunta, G.; Iannotta, D.A.; Montemurro, M. A FEM Free Vibration Analysis of Variable Stiffness Composite Plates through Hierarchical Modeling. *Materials* **2023**, *16*, 4643. https://doi.org/10.3390/ma16134643

Academic Editors: Angelo Marcello Tarantino, Carmelo Majorana, Raimondo Luciano and Michele Bacciocchi

Received: 17 May 2023

Revised: 14 June 2023

Accepted: 17 June 2023

Published: 27 June 2023

<figure><img src="image_4.png" /></figure>

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

**Abstract:** Variable Angle Tow (VAT) laminates offer a promising alternative to classical straight-fiber composites in terms of design and performance. However, analyzing these structures can be more complex due to the introduction of new design variables. Carrera's unified formulation (CUF) has been successful in previous works for buckling, vibrational, and stress analysis of VAT plates. Typically, one-dimensional (1D) and two-dimensional (2D) CUF models are used, with a linear law describing the fiber orientation variation in the main plane of the structure. The objective of this article is to expand the CUF 2D plate finite elements family to perform free vibration analysis of composite laminated plate structures with curvilinear fibers. The primary contribution is the application of Reissner's mixed variational theorem (RMVT) to a CUF finite element model. The principle of virtual displacements (PVD) and RMVT are both used as variational statements for the study of monolayer and multilayer VAT plate dynamic behavior. The proposed approach is compared to Abaqus three-dimensional (3D) reference solutions, classical theories and literature results to investigate the effectiveness of the developed models. The results demonstrate that mixed theories provide the best approximation of the reference solution in all cases.

**Keywords:** free vibration analysis; finite element method; variable angle tow plates; Carrera's unified formulation; Reissner's mixed variational theorem

## 1. Introduction

Over the last decades, composite structures have gained significant attention across diverse application fields, including aerospace, automotive and construction, due to their unique properties. Due to their high stiffness-to-weight ratio, composites help to build light structures with interesting mechanical properties. Despite this, a common thought is that the potential of fiber-reinforced structures could be better exploited by improving the directional properties through the variation of the fiber angle along the in-plane directions. The choice to keep the fiber orientation constant in each layer is particularly restrictive for geometries that present geometrical discontinuities such as cut-outs. VAT plates are characterized by an in-plane variation of fiber angle, helping to expand the design space of a specific structure. This is particularly useful for optimization problems, where a wider design space can positively affect the search of an optimal solution. For example, in the context of vibrational analyses, the maximization of fundamental frequencies can be improved by using curvilinear fibers. VATs were originally obtained through automated tape placement (ATP) and automated fiber placement (AFP). ATP helps the automated placement of composite material tapes with a specific angle in order to reproduce a desired path. AFP is similar to ATP, since the main difference is related to the width of the material that is laid down: while ATP handles a tape with a width between 75 and

<footnote n="1">¹ Luxembourg Institute of Science and Technology, 5 Avenue des Hauts-Fourneaux, L-4362 Esch-sur-Alzette, Luxembourg; domenico-andrea.iannotta@list.lu</footnote>

<footnote n="2">² Doctoral School in Science and Engineering (DSSE), University of Luxembourg, 2 Avenue de l'Université, L-4365 Esch-sur-Alzette, Luxembourg</footnote>

<footnote n="3">³ Arts et Métiers Institute of Technology, Université de Bordeaux, CNRS, INRA, Bordeaux INP, HESAM Université, I2M UMR 5295, F-33405 Talence, France; marco.montemurro@ensam.eu</footnote>

<footnote>* Correspondence: gaetano.giunta@list.lu</footnote>
300 mm, AFP involves the placement of the material with a typical width between 3.1 and 12.7 mm. By consequence, AFP allows for better control of fiber angles, achieving a wider design flexibility; see Dirk et al. [1]. However, automated processes show some limitations related to manufacturing defects, such as gaps and overlaps, or constraints such as the minimum steering radius. These aspects can be partially overcome due to new technologies such as additive manufacturing (AM), also known as 3D printing. AM involves the layer-by-layer deposition of materials to create a three-dimensional object. In the case of variable angle tow composites, AM techniques are used to deposit and cure layers of composite materials with varying fiber orientations; see Zhuo et al. [2]. VAT composites have diverse applications ranging from aerospace engineering and wind energy to automotive and construction contexts, offering enhanced structural performance, weight reduction, and tailored properties for improved efficiency and functionality in a wide range of industries. For example, these materials can be employed in the optimization of aircraft wings to enhance structural weight and fuel consumption, as presented in Brooks et al. [3]. In the space context, VATs can be used for the design of liquid oxygen in order to reduce the mass and increase the payload of space launchers, as discussed by Gren et al. [4]. Despite the significant advantages associated with curvilinear fiber composites, these materials have limitations from both manufacturing and design perspectives. The production of VAT composites can be more complex in comparison with traditional laminates, since the material behavior is strongly affected by process-induced defects. Moreover, not all fiber patterns can be realized, because of the technological limitations that characterize their production. The complexity of analysis is one of the main disadvantages of VATs, because a greater number of unknowns must be taken into account and unfeasible fiber patterns could be obtained during the optimization process.

Several methods for the study of VAT mechanical responses are available in the literature. In the following text, a brief review of these approaches is presented, with a particular focus on free vibration analyses. To the best of the authors' knowledge, the first works that have been presented on the topic are based on the assumption of a constant fiber angle within each element in a finite element method (FEM) solution. Therefore, the continuous variation of fiber direction was approximated in a step-wise discrete way. This approach can be used in commercial FEM software tools that, at the moment, cannot handle continuous fiber variation. Hyer and Charette [5] and Hyer and Lee [6] used this method to improve the VAT tensile strength and buckling response, respectively. One of the main disadvantages of this step-wise approach is that, as the true variation is continuous, the discrete representation of fiber angle variation imposes a further approximation. A p-version FEM based on the third-order shear deformation theory (TSDT) was applied by Akhavan and Ribeiro [7] to preform vibrational analyses. The results showed that fiber variation helps to increase (or decrease) natural frequencies and that thin plates are more affected by this phenomenon if compared with thick ones. Ribeiro and Akhavan [8] used the p-version FEM approach with elements based on the first-order shear deformation theory (FSDT) to perform non-linear vibration analyses. The advantage of the p-version of the FEM is that the accuracy of the approximation is improved by increasing the order of shape functions over the elements. Vibration analyses were performed on VAT plates with a central circular cut-out considering parabolic fibers by Hachemi et al. [9]. Zhao and Kapania [10] investigated the free vibration of prestressed VAT stiffened plates, where plates and stiffeners were modeled separately through Mindlin plate theory and Timoshenko beam theory, respectively. The compatibility conditions at the interface between the plate and stiffeners were satisfied by using a transformation matrix. Honda and Narita [11] used the classical plate theory within the Ritz method in order to evaluate the natural frequencies and vibrational modes. An experimental approach was used in Rodrigues et al. [12] for the free vibration analysis of a plate with free boundary conditions that was subjected to random excitation via an electromagnetic shaker. Subsequently, the results were compared to the ones obtained through FEM, where a four-node isoparametric element based on the Reissner–Mindlin theory was used. Stodieck et al. [13] showed that curvilinear fibers
can be useful for improving the aeroelastic response of composite wings. The Rayleigh–Ritz method and classical lamination theory were used to develop a 1D beam model, considering the assumption of null chamber deformation of the wing chord-wise section. The aeroelastic response was computed by introducing quasi-static aerodynamic forces in a model developed for the plate structural analysis. A parametric study showed that by using VATs, it is possible to influence wing response both positively and negatively.

Curvilinear fibers can improve the modal response, as shown in several works. Abdalla et al. [14] used the classical lamination theory in combination with a successive approximation method in order to solve an optimization problem. The results showed that curvilinear fibers increased the optimal fundamental frequency in comparison with straight ones. A similar approach was presented in Blom et al. [15], where the maximization of the first natural frequency considering manufacturing constraints was obtained for VAT conical shells. In Carvalho et al. [16], a genetic algorithm and shell elements based on FSDT were used for maximization of the fundamental frequency. The multi-scale two-level (MS2L) approach helps to split the optimization problem in two parts. The composite is modeled as an equivalent homogeneous anisotropic plate in the first step, which aims to find the ideal distribution of the polar parameters that represent the mechanical design variables. The main goal of a second step is to establish the best stacking sequence in relation to the mechanical property distribution that has been obtained in the first step. The MS2L method was applied by Montemurro and Catapano [17] to VAT plates in order to optimize the buckling response. In order to evaluate the polar parameters, B-spline surfaces were introduced, while manufacturing constraints were considered during the second step. More details about the MS2L approach can be found in Catapano et al. [18], Montemurro and Catapano [19] and Fiordilino et al. [20], where both stiffness and buckling optimization problems were solved.

VAT structures have also been studied by using Carrera's unified formulation. CUF is a mathematical framework that helps the derivation of different theories, such as classical lamination plate theory, higher-order shear deformation theories, or LW approaches, within a unique formulation; see Carrera [21,22]. The a priori approximation over the thickness (typical of plates' structural modeling) can also be freely assumed as a generic combination of functions whose number is a free parameter of the formulation. When polynomial functions are used, as in this article, the expansion order along the thickness of the plate is arbitrary in the formulation, and it can be set when performing a specific analysis. This flexibility is beneficial because it helps to tailor the accuracy and computational efficiency of the analysis to the specific requirements of the problem at hand. Carrera et al. [23] used CUF in order to develop a Navier closed-form solution for the static analysis of isotropic plates under several loading conditions. The same approach was used in Carrera and Giunta [24] in order to perform failure analyses on isotropic plates. A further extension of this method was shown in Giunta et al. [25], where an indentation failure analysis of composite sandwich plates was performed. Giunta et al. [26] performed free vibration analyses of composite beams. In Viglietti et al. [27] and Fallahi et al. [28], free vibration and buckling analyses of VATs were performed through the use of a 1D CUF model. Within this framework, shell models were developed as well for VAT cases in order to perform stress analyses; see Sánchez-Majano et al. [29]. In Pagani and Sánchez-Majano [30,31] and Sánchez-Majano et al. [32], manufacturing defects were taken into account by using stochastic techniques. Vescovini and Dozio [33] used the Ritz method within CUF for vibrational and buckling analyses. A generalization of CUF was developed in order to allow for the use of different expansions for every component of the displacement vector. Demasi et al. [34] applied this approach to the study of VAT plates with an ESL model. A further advantage of CUF is that it can be used in combination with different variational statements. An alternative to the classic PVD is represented by the RMVT, where both displacements and transverse out-of-plane stresses are considered as primary variables. RMVT has been widely used within CUF for the study of straight-fiber composite structures.
For example, Carrera and Demasi [35,36] developed RMVT-based CUF models to perform the static analysis of straight fiber plates.

The free vibration analysis is an important problem in engineering (see Babaei et al. [37]), and within this context, CUF has been applied to the study of VATs considering PVD as the main variational statement. For this reason, this work aims to extend this framework with the RMVT formulation in order to develop a family of hierarchical plate finite elements. This will help to better predict the natural frequencies of composite plates characterized by curvilinear fibers. Section 2 shows the theoretical derivation for free vibration problems. Section 3 presents the numerical results where three cases are investigated. Analyses are performed that consider a varying side-to-thickness ratio in order to investigate thin and thick plates, and the differences between models are discussed regarding PVD or RMVT statements. The results are compared to reference solutions for validation. Concluding observations and remarks are presented in Section 4.

## 2. Carrera’s Unified Formulation

A plate is a flat body whose material points lie in the Cartesian closed-point subset

$$ \mathcal{P} = \Omega \times \mathcal{H} \qquad (1) $$

of the three-dimensional space $\mathbb{R}^3$ where:

$$ \begin{aligned} \Omega &= \{(x, y) : \frac{x}{a}, \frac{y}{b} \in [0, 1]\} \subset \mathbb{R}^2, \\ \mathcal{H} &= \{z : \frac{2z}{h} \in [-1, 1]\}, \end{aligned} \qquad (2) $$

where $a$ and $b$ are the dimensions along the two in-plane axes, and $h$ measures its thickness along the z-axis, where $z \ll a$ and $b$. The global reference system and plate geometry are presented in Figure 1.

<figure><img src="image_5.png" /><figcaption>Figure 1. Plate geometry and reference system.</figcaption></figure>

The displacement field is expressed as:

$$ \mathbf{u} = \begin{cases} u_x \\ u_y \\ u_z \end{cases}. \qquad (3) $$

The strain tensor components can be written in vector form. Two vectors are obtained, representing the in-plane and out-of-plane components:

$$ \epsilon_p = \begin{cases} \epsilon_{xx} \\ \epsilon_{yy} \\ \epsilon_{xy} \end{cases}, \quad \epsilon_n = \begin{cases} \gamma_{xz} \\ \gamma_{yz} \\ \gamma_{zz} \end{cases}. \qquad (4) $$
The hypothesis of small displacements helps to use a linear strain-displacement relation:

$$ \epsilon_p = \mathbf{D}_p \mathbf{u}, \\ \epsilon_n = (\mathbf{D}_{n\Omega} + \mathbf{D}_{nz}) \mathbf{u}, \qquad (5) $$

where $\mathbf{D}_p$, $\mathbf{D}_{n\Omega}$ and $\mathbf{D}_{nz}$ are the following differential operators:

$$ \mathbf{D}_p = \begin{bmatrix} \frac{\partial}{\partial x} & 0 & 0 \\ 0 & \frac{\partial}{\partial y} & 0 \\ \frac{\partial}{\partial y} & \frac{\partial}{\partial x} & 0 \end{bmatrix}, \quad \mathbf{D}_{n\Omega} = \begin{bmatrix} 0 & 0 & \frac{\partial}{\partial x} \\ 0 & 0 & \frac{\partial}{\partial y} \\ 0 & 0 & 0 \end{bmatrix}, \quad \mathbf{D}_{nz} = \begin{bmatrix} \frac{\partial}{\partial z} & 0 & 0 \\ 0 & \frac{\partial}{\partial z} & 0 \\ 0 & 0 & \frac{\partial}{\partial z} \end{bmatrix}. \quad (6) $$

The stress vector is expressed in a similar manner:

$$ \sigma_p = \begin{cases} \sigma_{xx} \\ \sigma_{yy} \\ \sigma_{xy} \end{cases}, \quad \sigma_n = \begin{cases} \sigma_{xz} \\ \sigma_{yz} \\ \sigma_{zz} \end{cases}. \qquad (7) $$

Hooke's law reads:

$$ \sigma_p = \tilde{C}_{pp} \epsilon_p + \tilde{C}_{pn} \epsilon_n, \\ \sigma_n = \tilde{C}_{np} \epsilon_p + \tilde{C}_{nn} \epsilon_n, \qquad (8) $$

where the terms $\tilde{C}_{pp}$, $\tilde{C}_{pn}$, $\tilde{C}_{np}$ and $\tilde{C}_{nn}$ are subcomponents of a material stiffness matrix $\tilde{C}$ according to the stress and strain ordering in Equations (4) and (7), where the fibers lay in $\Omega$ and where they are not, in general, aligned with the x-axis. C stands for the stiffness matrix in the global reference system, and its components can be written in terms of the Young's moduli $E_L$ and $E_T$, shear moduli $G_{LT}$ and $G_{TT}$ and Poisson's ratios $\nu_{LT}$ and $\nu_{TT}$, where subscripts L and T stand for the directions parallel and perpendicular to the fibers, respectively. For further details, see Reddy [38].

## 2.1. Variable Stiffness Composite Plates

Laminated VAT structures are considered in this work. For this reason, the material stiffness coefficients can change layer-wise along the thickness and pointwise along the in-plane directions. The mapping of C into $\tilde{C}$ reads:

$$ \tilde{C} = \mathbf{TCT}^T. \qquad (9) $$

Superscript T stands for the transpose operator. The matrix T represents a rotation matrix that depends on an in-plane rotation angle $\theta$. For the sake of brevity, the components of $\tilde{C}$ and T are not reported here; they can be found in Reddy [38]. In a laminated VAT, the rotation angle $\theta$ is a bi-dimensional field in $\Omega$. In this work, two different variation laws are considered for $\theta$, a linear variation law and a parabolic one. The linear law can be expressed according to the following formula:

$$ \theta(\alpha) = \Phi + T_0 + \frac{T_1 - T_0}{d} |\alpha|. \qquad (10) $$

The angle $\Phi$ describes the original direction along which $\theta$ varies, and $\alpha$ is a generic spatial variable defined as:

$$ \alpha = x' \cos(\Phi) + y' \sin(\Phi). \qquad (11) $$

$x'$ and $y'$ denote a generic in-plane reference system used for describing a fiber path, where $\theta$ is measured. The introduction of a new reference system is useful in order to represent the local fiber orientation independently from the global reference system identified by axes x and y. $T_0$ and $T_1$ are the angles between the a-axis and the tangent to a fiber for $\alpha$ equal to zero and d, respectively; see Figure 2.
<figure><img src="image_6.png" /><figcaption>Figure 2. Example of in-plane fiber orientation.</figcaption></figure>

As shown in the Figure 2, the fiber angle is always measured with respect to the x'-axis, and it can change along a generic direction $\alpha$, defined as a combination of x' and y' depending on the angle $\Phi$. Further details about the fiber linear variation law can be found in Gürdal et al. [39]. The parabolic law can be expressed according to the following equation:

$$ \theta(\alpha) = \Phi + T_0 + \tan^{-1}\left(\gamma \frac{\alpha}{d}\right). \qquad (12) $$

The parameter $\gamma$ is used to control the shape of the parabola, and it is related to the final fiber angle $T_1$ as $T_1 = \tan^{-1}(\pm\gamma)$. More details about the parabolic fiber path can be found in Hachemi et al. [9] and Honda et al. [40]. The following notation, based upon the above introduced parameters, is used in order to describe the in-plane linear and parabolic fiber behavior: $\Phi < T_0, T_1 >$.

## 2.2. Variational Statements

PVD and RMVT variational statements are considered to derive the governing equations for the free vibration problem for a laminated VAT plate. The fundamental distinction is that the RMVT considers the vector of the out-of-plane stresses $\sigma_n$ as a primary unknown, whereas the PVD considers only displacements as primary variables. For the PVD case, the following variational statement applies:

$$ \int_{\Omega} \int_{H} \left( \delta \epsilon_{pG}^{T} \sigma_{pH} + \delta \epsilon_{nG}^{T} \sigma_{nH} \right) dz d\Omega + \delta \mathcal{L}_{in} = 0, \qquad (13) $$

where the subscript G refers to the components obtained from the geometrical relations in Equation (5), and subscript H refers to the components obtained from Hooke's law in Equation (8). $\mathcal{L}_{in}$ is the virtual work of the inertial forces, and $\delta$ stands for a virtual variation. For the RMVT case, the variational statement is:

$$ \int_{\Omega} \int_{H} \left[ \delta \epsilon_{pG}^{T} \sigma_{pH} + \delta \epsilon_{nG}^{T} \sigma_{nM} + \delta \sigma_{nM}^{T} (\epsilon_{nG} - \epsilon_{nH}) \right] dz d\Omega + \delta \mathcal{L}_{in} = 0. \qquad (14) $$

The M subscript refers to the transverse stress components considered as primary unknowns in the mixed formulation. For the RMVT formulation, Hooke's law is rewritten as follows:

$$ \begin{aligned} \sigma_{pH} &= \hat{C}_{pp} \epsilon_{pG} + \hat{C}_{pn} \sigma_{nM}, \\ \epsilon_{nH} &= \hat{C}_{np} \epsilon_{pG} + \hat{C}_{nn} \sigma_{nM}, \end{aligned} \qquad (15) $$

where $\hat{C}_{pp}$, $\hat{C}_{pn}$, $\hat{C}_{np}$ and $\hat{C}_{nn}$ are (see Carrera and Demasi [35]):

$$ \begin{aligned} \hat{C}_{pp} &= \tilde{C}_{pp} - \tilde{C}_{pn}\tilde{C}_{nn}^{-1}\tilde{C}_{np}, \\ \hat{C}_{pn} &= \tilde{C}_{pn}\tilde{C}_{nn}^{-1}, \\ \hat{C}_{np} &= -\tilde{C}_{nn}^{-1}\tilde{C}_{np}, \\ \hat{C}_{nn} &= \tilde{C}_{nn}^{-1}. \end{aligned} \qquad (16) $$
The superscript “-1” indicates the inverse of a matrix. The inertial work can be expressed as:

$$ \delta\mathcal{L}_{in} = \int_{\Omega} \int_{\mathcal{H}} \delta\mathbf{u}^T \rho \ddot{\mathbf{u}} d\Omega dz, \qquad (17) $$

where $\rho$ is the plate material density, and $\ddot{\mathbf{u}}$ represents the acceleration vector.

### 2.3. Kinematic Assumptions

CUF uses an axiomatic approach along the through-the-thickness direction to represent the primary unknowns; see Carrera [22]. The generic unknown component $f = f(x, y, z)$ is approximated as:

$$ f(x, y, z) = F_{\tau}(z)g_{\tau}(x, y), \quad \tau = 0, 1, \dots, N, \qquad (18) $$

where $f$ is a displacement component in a formulation derived from the PVD, but it can also be an out-of-plane stress component when a RMVT formulation is considered. $F_{\tau}$ is an approximation function in $\mathcal{H}$, and $g_{\tau}$ is an unknown two-dimensional function in $\Omega$. According to Einstein's notation, a twice-repeated index implies a sum over the index range. Finally, $N$ is the approximation order. Both $N$ and $F_{\tau}$ are a priori defined. This feature of CUF helps to obtain multiple theories in the same formulation. Within CUF, ESL or LW models can also be obtained depending on the support of $F_{\tau}$. In an ESL model

$$ F_{\tau}: \mathcal{H} \mapsto \mathbb{R}, \text{ whereas for a LW model } F_{\tau}: \mathcal{H}^k \mapsto \mathbb{R} \text{ where } \mathcal{H}^k = \left\{ z^k : \frac{2z^k}{h^k} \in [-1, 1] \right\} $$

such that $\mathcal{H} = \bigcup_{k=1}^{N_l} \mathcal{H}^k$ and $\mathcal{H}^k \cap \mathcal{H}^{k'} = \emptyset$ for $k \neq k'$ with $k, k' = 1, 2, \dots, N_l$, where $N_l$ is the total number of laminae, and $h^k$ is the thickness of a generic $k$ lamina such that

$k = \sum_{k=1}^{N_l} h^k$. The number of unknowns in the ESL case is independent of the number of layers in the lamination since the approximation is imposed globally over $\mathcal{H}$. The total stiffness contributions can be seen as a weighted average of each layer stiffness along the thickness. Maclaurin's series approximation is considered for the ESL models as a linear combination

of the power functions:

$$ F_{\tau}(z) = z^{\tau}, \quad \tau = 0, 1, \dots, N, \qquad (19) $$

where $N$ is the expansion order. The computational cost of ESL models depends on $N$ only, and for a given $N$, it is lower than a LW model since this latter model depends on the total number of layers in the lamination. ESL models are suitable for relatively thick laminates. However, they are unable to accurately predict the behavior of thick plates with a high degree of anisotropy. ESL models have $C^{\infty}$-continuity over $\mathcal{H}$ because of the used approximation functions, whereas laminated composites present a $C^0$-continuity since the interface between the two consecutive layers of the different materials introduces a change in the slope of the displacements (also known as zig-zag displacement through-the-thickness variation). This behavior can be accommodated within an ESL theory by means of Murakami's function. This approach is not considered here; for more details, refer to Carrera [41]. In an LW model, the kinematics of each layer are formulated independently:

$$ f^k(x, y, z) = F_b(z)g_b^k(x, y) + F_r(z)g_r^k(x, y) + F_t(z)g_t^k(x, y), \quad r = 2, \dots, N, \qquad (20) $$

where subscripts $b$ and $t$ stand for the bottom and top layers, respectively. Congruence at the interface is retrieved via a through-the-thickness assembly procedure similar to that used in the finite element method. For this reason, Lagrange polynomials (which ensure partition of unity), or the following linear combination of Legendre polynomials, which are represented as:

$$ F_t(z(\xi^k)) = \frac{P_0 + P_1}{2}, \quad F_b(z(\xi^k)) = \frac{P_0 - P_1}{2}, \quad F_r(z(\xi^k)) = P_r - P_{r-2}, \quad r = 2, \dots, N \qquad (21) $$
are typically used as approximation functions over $\mathcal{H}^k$. The use of Lagrange or Legendre polynomials along the thickness changes according to the used model, and this is specified at the end of the next subsection. In Equation (21), $\xi^k = \frac{2z^k}{\hbar^k} \in [-1, 1]$ and $P_i = P_i(\xi^k)$ are an $i$-order Legendre polynomial. Equation (21) creates a base where $F_t$ and $F_b$ are the two linear Lagrange polynomials, and $F_r$ is a kind of $p$-version-enriching function since it does not contribute to a base linear combination for $\xi^k = \pm 1$, being, by definition, $F_r(\pm 1) = 0$. Since LW base functions have local support, inter-layer $C^0$-continuity for layers made of different materials is ensured, but the computational costs are higher than for ESL models.

## 2.4. Acronym System

An acronym system is used in order to identify all the derived theories. Figure 3 shows this system.

<figure><img src="image_7.png" /><figcaption>Figure 3. Acronym system.</figcaption></figure>

The first letter addresses the approximation level that is applied: ‘E’ denotes the ESL models, whereas ‘L’ denotes the LW models. The second letter denotes the variational statement: PVD or RMVT are denoted by ‘D’ or ‘M’, respectively. The last number is the order of expansion along the plate thickness. A number at the beginning of the acronym, when present, indicates how many virtual layers have been used to approximate each physical layer in an LW model to improve the results for a given approximation order. If this number is not present, only one virtual layer has been used to represent each physical layer.

As an example, in EDN models, the displacement field can be expressed as:

$$ 
\begin{aligned}
u_x &= u_{x0} + u_{x1}z + u_{x2}z^2 + \cdots + u_{xN}z^N, \\
u_y &= u_{y0} + u_{y1}z + u_{y2}z^2 + \cdots + u_{yN}z^N, \\
u_z &= u_{z0} + u_{z1}z + u_{z2}z^2 + \cdots + u_{zN}z^N.
\end{aligned} 
\quad (22) $$

In vector form:

$$ \mathbf{u} = F_0\mathbf{u}_0 + F_1\mathbf{u}_1 + \cdots + F_N\mathbf{u}_N = F_{\tau}\mathbf{u}_{\tau}, \quad \tau = 0, 1, \dots, N, \quad (23) $$

where $F_{\tau} = z^{\tau}$ and $\mathbf{u}_{\tau} = \mathbf{u}_{\tau}(x, y)$. Additionally, classical theories can be taken into account. Classical lamination theory (CLT) and first-order shear deformation theory are obtained as a particular case of a first-order ESL theory. FSDT is obtained through the penalization of the $uz_1$ term, while for CLT, transverse shear stresses are disregarded by using a fictitiously
high value of the material shear moduli. The material stiffness matrix needs to be reduced in a plane-stress sense to overcome thickness locking.

For LDN solutions, only displacements are considered as the primary variables:

$$ \mathbf{u}^k = F_0 \mathbf{u}_0^k + F_1 \mathbf{u}_1^k + \dots + F_N \mathbf{u}_N^k = F_\tau \mathbf{u}_\tau^k, \quad \tau = 0, 1, \dots, N, \quad k = 1, 2, \dots, N_l. \quad (24) $$

For LMN solutions, transverse stresses are treated as primary variables. The transverse stress field can be expressed as:

$$ \sigma_n^k = F_0 \sigma_0^k + F_1 \sigma_1^k + \dots + F_N \sigma_N^k = F_\tau \sigma_\tau^k, \quad \tau = 0, 1, \dots, N, \quad k = 1, 2, \dots, N_l. \quad (25) $$

It can be observed that ESL theories can be considered as a particular case for LW theories. While in the first case the integration along the thickness is performed in order to represent composite properties through a unitary layer, for the second case, the integration is computed layer by layer. This helps to represent the kinematics of each layer separately for LW models. LDN solutions are obtained with Lagrange polynomials with equally spaced nodes, whereas LMN ones are obtained with Legendre polynomials.

## 2.5. FE Stiffness Matrices

As far as a FEM solution is concerned, the in-plane domain is discretized into $N_e$ subdomains such as $\Omega = \bigcup_{e=1}^{N_e} \Omega_e$ and $\Omega_e \cap \Omega' = \emptyset$ for $e \neq e'$. Shape functions are then introduced for the approximation of the variation over $\Omega_e$. In the case of a bi-dimensional model, Equation (18) becomes:

$$ f(x, y, z) = F_{\tau}(z) N_i(x, y) g_{\tau i}, \quad \tau = 0, 1, \dots, N, \quad i = 1, \dots, N_n, \quad (26) $$

where $N_i$ stands for the shape functions, and $N_n$ is the number of nodes in the used finite element. Classical Lagrange shape functions are used. They are not presented here for the sake of brevity, but they can be found in Bathe [42]. FE stiffness matrices are obtained by the weak form of the variational principles. In the PVD case, considering Equation (26), the displacement field can be written as:

$$ \mathbf{u} = F_{\tau} N_i \begin{Bmatrix} q_{x\tau i} \\ q_{y\tau i} \\ q_{z\tau i} \end{Bmatrix} = F_{\tau} N_i \mathbf{q}_{\tau i}. \quad (27) $$

Through the substitution of Equations (5), (8) and (27) into Equation (13), the weak PVD form can be obtained:

$$ \begin{aligned} & \int_{\Omega_e} \delta q_{\tau i}^T [\mathbf{D}_p^T(N_i \mathbf{I}) \tilde{\mathbf{Z}}_{pp}^{\tau s} \mathbf{D}_p(N_j \mathbf{I}) + \mathbf{D}_p^T(N_i \mathbf{I}) \tilde{\mathbf{Z}}_{pn}^{\tau s} \mathbf{D}_{n\Omega}(N_j \mathbf{I}) + \mathbf{D}_p^T(N_i \mathbf{I}) \tilde{\mathbf{Z}}_{pn}^{\tau s,z}(N_j \mathbf{I}) \\ & + \mathbf{D}_{n\Omega}^T(N_i \mathbf{I}) \tilde{\mathbf{Z}}_{np}^{\tau s} \mathbf{D}_p(N_j \mathbf{I}) + \mathbf{D}_{n\Omega}^T(N_i \mathbf{I}) \tilde{\mathbf{Z}}_{nn}^{\tau s} \mathbf{D}_{n\Omega}(N_j \mathbf{I}) + \mathbf{D}_{n\Omega}^T(N_i \mathbf{I}) \tilde{\mathbf{Z}}_{nn}^{\tau s,z}(N_j \mathbf{I}) \\ & + (N_i \mathbf{I}) \tilde{\mathbf{Z}}_{np}^{\tau,z s} \mathbf{D}_p(N_j \mathbf{I}) + (N_i \mathbf{I}) \tilde{\mathbf{Z}}_{nn}^{\tau,z s} \mathbf{D}_{n\Omega}(N_j \mathbf{I}) + (N_i \mathbf{I}) \tilde{\mathbf{Z}}_{nn}^{\tau,z s,z}(N_j \mathbf{I})] \\ & q_{sj} d\Omega = - \int_{\Omega_e} \delta q_{\tau i}^T (N_i \mathbf{I}) \rho E_{\tau s} (N_j \mathbf{I}) \ddot{q}_{sj} d\Omega, \end{aligned} \quad (28) $$

where:

$$ (\tilde{\mathbf{Z}}_{wr}^{\tau s}, \tilde{\mathbf{Z}}_{wr}^{\tau,z s}, \tilde{\mathbf{Z}}_{wr}^{\tau s,z}, \tilde{\mathbf{Z}}_{wr}^{\tau,z s,z}) = (\tilde{\mathbf{C}}_{wr} E_{\tau s}, \tilde{\mathbf{C}}_{wr} E_{\tau,z s}, \tilde{\mathbf{C}}_{wr} E_{\tau,s,z}, \tilde{\mathbf{C}}_{wr} E_{\tau,z s,z}): w, r = p, n, \quad (29) $$

$$ (E_{\tau s}, E_{\tau,z s}, E_{\tau,s,z}, E_{\tau,z s,z}) = \int_{\mathcal{H}} (F_{\tau} F_s, F_{\tau,z} F_s, F_{\tau} F_{s,z}, F_{\tau,z} F_{s,z}) dz. \quad (30) $$

An axis coordinate as comma-preceded subscript stands for a derivative in that coordinate direction. In compact vector form, Equation (28) reads:

$$ \delta q_{\tau i}^T \mathbf{K}^{\tau sij} q_{sj} = -\delta q_{\tau i}^T \mathbf{M}^{\tau sij} \ddot{q}_{sj}, \quad (31) $$
where $K^{\tau sij}$ and $M^{\tau sij} \in \mathbb{R}^{3 \times 3}$ are fundamental nuclei (FN) of the stiffness and mass matrices, respectively. Through the cycles on the indices $\tau, s, i$ and $j$, it is possible to build the stiffness and mass matrices of a finite element. The components of the stiffness FN for the PVD case can be written as:

$$K_{xx}^{\tau sij} = \int_{\Omega_e} (\tilde{Z}_{pp11}^{\tau s} N_{j_{ix}} N_{i_{ix}} + \tilde{Z}_{pp16}^{\tau s} N_{j_{iy}} N_{i_{iy}} + \tilde{Z}_{pp16}^{\tau s} N_{j_{ix}} N_{i_{ix}} + \tilde{Z}_{pp66}^{\tau s} N_{j_{iy}} N_{i_{iy}} + \tilde{Z}_{nn44}^{\tau s, z_{z}} N_j N_i) d\Omega,$$

$$K_{xy}^{\tau sij} = \int_{\Omega_e} (\tilde{Z}_{pp12}^{\tau s} N_{j_{iy}} N_{i_{ix}} + \tilde{Z}_{pp16}^{\tau s} N_{j_{ix}} N_{i_{ix}} + \tilde{Z}_{pp26}^{\tau s} N_{j_{iy}} N_{i_{iy}} + \tilde{Z}_{pp66}^{\tau s} N_{j_{ix}} N_{i_{iy}} + \tilde{Z}_{nn45}^{\tau s, z_{z}} N_j N_i) d\Omega,$$

$$K_{xz}^{\tau sij} = \int_{\Omega_e} (\tilde{Z}_{pn13}^{\tau s, z} N_j N_{i_{ix}} + \tilde{Z}_{pn36}^{\tau s, z} N_j N_{i_{iy}} + \tilde{Z}_{nn44}^{\tau s, z} N_{j_{ix}} N_i + \tilde{Z}_{nn45}^{\tau s, z} N_{j_{iy}} N_i) d\Omega,$$

$$K_{yx}^{\tau sij} = \int_{\Omega_e} (\tilde{Z}_{pp12}^{\tau s} N_{j_{ix}} N_{i_{iy}} + \tilde{Z}_{pp26}^{\tau s} N_{j_{iy}} N_{i_{iy}} + \tilde{Z}_{pp16}^{\tau s} N_{j_{ix}} N_{i_{ix}} + \tilde{Z}_{pp66}^{\tau s} N_{j_{iy}} N_{i_{ix}} + \tilde{Z}_{nn45}^{\tau s, z} N_j N_i) d\Omega,$$

$$K_{yy}^{\tau sij} = \int_{\Omega_e} (\tilde{Z}_{pp22}^{\tau s} N_{j_{iy}} N_{i_{iy}} + \tilde{Z}_{pp26}^{\tau s} N_{j_{iy}} N_{i_{ix}} + \tilde{Z}_{pp66}^{\tau s} N_{j_{ix}} N_{i_{ix}} + \tilde{Z}_{nn55}^{\tau s, z} N_j N_i) d\Omega, \quad (32)$$

$$K_{yz}^{\tau sij} = \int_{\Omega_e} (\tilde{Z}_{pn23}^{\tau s, z} N_j N_{i_{iy}} + \tilde{Z}_{pn36}^{\tau s, z} N_j N_{i_{ix}} + \tilde{Z}_{nn45}^{\tau s, z} N_{j_{ix}} N_i + \tilde{Z}_{nn55}^{\tau s, z} N_{j_{iy}} N_i) d\Omega,$$

$$K_{zx}^{\tau sij} = \int_{\Omega_e} (\tilde{Z}_{nn44}^{\tau s, z} N_j N_{i_{ix}} + \tilde{Z}_{nn45}^{\tau s, z} N_j N_{i_{iy}} + \tilde{Z}_{np13}^{\tau s, z} N_{j_{ix}} N_i + \tilde{Z}_{np36}^{\tau s, z} N_{j_{iy}} N_i) d\Omega,$$

$$K_{zy}^{\tau sij} = \int_{\Omega_e} (\tilde{Z}_{nn45}^{\tau s, z} N_j N_{i_{ix}} + \tilde{Z}_{nn55}^{\tau s, z} N_j N_{i_{iy}} + \tilde{Z}_{np23}^{\tau s, z} N_{j_{iy}} N_i + \tilde{Z}_{np36}^{\tau s, z} N_{j_{ix}} N_i) d\Omega,$$

$$K_{zz}^{\tau sij} = \int_{\Omega_e} (\tilde{Z}_{nn44}^{\tau s} N_{j_{ix}} N_{i_{ix}} + \tilde{Z}_{nn45}^{\tau s} N_{j_{iy}} N_{i_{ix}} + \tilde{Z}_{nn45}^{\tau s} N_{j_{ix}} N_{i_{iy}} + \tilde{Z}_{nn55}^{\tau s, z} N_j N_i + \tilde{Z}_{nn33}^{\tau s, z} N_j N_i) d\Omega.$$

The mass FN can be written as:

$$M^{\tau sij} = \int_{\Omega_e} (N_i \mathbf{I}) \rho E_{\tau s} (N_j \mathbf{I}) d\Omega. \quad (33)$$

It is possible to observe that $M^{\tau sij}$ is a diagonal matrix and that since the plate density is assumed to be constant, the term $\rho E_{\tau s}$ can be placed outside the integral.

In the RMVT case, transverse stresses are a priori approximated:

$$\sigma_n = F_{\tau} N_i \begin{Bmatrix} g_{xz\tau i} \\ g_{yz\tau i} \\ g_{zz\tau i} \end{Bmatrix} = F_{\tau} N_i g_{\tau i}. \quad (34)$$

Through the substitution of Equations (5), (15), (27) and (34) into Equation (14), the RMVT weak form can be obtained:

$$\int_{\Omega_e} \delta \boldsymbol{q}_{\tau i}^T [\mathbf{D}_p^T (N_i \mathbf{I}) \hat{\boldsymbol{Z}}_{pp}^{\tau s} \mathbf{D}_p (N_j \mathbf{I})] \boldsymbol{q}_{sj} + \delta \boldsymbol{q}_{\tau i}^T [\mathbf{D}_p^T (N_i \mathbf{I}) \hat{\boldsymbol{Z}}_{pn}^{\tau s} (N_j \mathbf{I}) + \mathbf{D}_{n\Omega}^T (N_i \mathbf{I}) (E_{\tau s} \mathbf{I}) (N_j \mathbf{I}) \\ + (N_i \mathbf{I}) (E_{\tau z s} \mathbf{I}) (N_j \mathbf{I})] \boldsymbol{g}_{sj} + \delta \boldsymbol{g}_{\tau i}^T [(N_i \mathbf{I}) (E_{\tau s} \mathbf{I}) \mathbf{D}_{n\Omega} (N_j \mathbf{I}) + (N_i \mathbf{I}) (E_{\tau s, z} \mathbf{I}) (N_j \mathbf{I}) \\ -(N_i \mathbf{I}) \hat{\boldsymbol{Z}}_{np}^{\tau s} \mathbf{D}_p (N_j \mathbf{I})] \boldsymbol{q}_{sj} - \delta \boldsymbol{g}_{\tau i}^T (N_i \mathbf{I}) \hat{\boldsymbol{Z}}_{nn}^{\tau s} (N_j \mathbf{I}) \boldsymbol{g}_{sj} d\Omega = - \int_{\Omega_e} \delta \boldsymbol{q}_{\tau i}^T (N_i \mathbf{I}) \rho E_{\tau s} (N_j \mathbf{I}) \ddot{\boldsymbol{q}}_{sj} d\Omega, \quad (35)$$

where:

$$(\hat{\boldsymbol{Z}}_{wr}^{\tau s}, \hat{\boldsymbol{Z}}_{wr}^{\tau z s}, \hat{\boldsymbol{Z}}_{wr}^{\tau s, z}, \hat{\boldsymbol{Z}}_{wr}^{\tau z, s z}) = (\hat{\mathbf{C}}_{wr} E_{\tau s}, \hat{\mathbf{C}}_{wr} E_{\tau z s}, \hat{\mathbf{C}}_{wr} E_{\tau s, z}, \hat{\mathbf{C}}_{wr} E_{\tau z, s z}): w, r = p, n. \quad (36)$$

In a compact form:

$$\begin{aligned} \delta \boldsymbol{q}_{\tau i}^T \mathbf{K}_{uu}^{\tau sij} \boldsymbol{q}_{sj} + \delta \boldsymbol{q}_{\tau i}^T \mathbf{K}_{u\sigma}^{\tau sij} \boldsymbol{g}_{sj} &= -\delta \boldsymbol{q}_{\tau i}^T \mathbf{M}^{\tau sij} \ddot{\boldsymbol{q}}_{sj}, \\ \delta \boldsymbol{g}_{\tau i}^T \mathbf{K}_{\sigma u}^{\tau sij} \boldsymbol{q}_{sj} + \delta \boldsymbol{g}_{\tau i}^T \mathbf{K}_{\sigma \sigma}^{\tau sij} \boldsymbol{g}_{sj} &= 0. \end{aligned} \quad (37)$$
In this case, four fundamental nuclei are obtained. The components of the FN for the RMVT case can be written as:

$$K_{uuxx}^{\tau sij} = \int_{\Omega_e} (\hat{Z}_{pp11}^{\tau s} N_{j,x} N_{i,x} + \hat{Z}_{pp31}^{\tau s} N_{j,x} N_{i,y} + \hat{Z}_{pp13}^{\tau s} N_{j,y} N_{i,x} + \hat{Z}_{pp33}^{\tau s} N_{j,y} N_{i,y}) d\Omega,$$

$$K_{uuxy}^{\tau sij} = \int_{\Omega_e} (\hat{Z}_{pp12}^{\tau s} N_{j,y} N_{i,x} + \hat{Z}_{pp32}^{\tau s} N_{j,y} N_{i,y} + \hat{Z}_{pp13}^{\tau s} N_{j,x} N_{i,x} + \hat{Z}_{pp33}^{\tau s} N_{j,x} N_{i,y}) d\Omega,$$

$$K_{uuyx}^{\tau sij} = \int_{\Omega_e} (\hat{Z}_{pp21}^{\tau s} N_{j,x} N_{i,y} + \hat{Z}_{pp31}^{\tau s} N_{j,x} N_{i,x} + \hat{Z}_{pp23}^{\tau s} N_{j,y} N_{i,y} + \hat{Z}_{pp33}^{\tau s} N_{j,y} N_{i,x}) d\Omega,$$

$$K_{uuyy}^{\tau sij} = \int_{\Omega_e} (\hat{Z}_{pp22}^{\tau s} N_{j,y} N_{i,y} + \hat{Z}_{pp32}^{\tau s} N_{j,y} N_{i,x} + \hat{Z}_{pp23}^{\tau s} N_{j,x} N_{i,y} + \hat{Z}_{pp33}^{\tau s} N_{j,x} N_{i,x}) d\Omega,$$

$$K_{uuxz}^{\tau sij} = 0, \quad K_{uuyz}^{\tau sij} = 0, \quad K_{uuxz}^{\tau sij} = 0, \quad K_{uuzy}^{\tau sij} = 0, \quad K_{uuzz}^{\tau sij} = 0,$$

$$K_{u\sigma xx}^{\tau sij} = \int_{\Omega_e} (E_{\tau,zs} N_j N_i) d\Omega, \quad K_{u\sigma xz}^{\tau sij} = \int_{\Omega_e} (\hat{Z}_{pm13}^{\tau s} N_j N_{i,x} + \hat{Z}_{pn33}^{\tau s} N_j N_{i,y}) d\Omega,$$

$$K_{u\sigma yy}^{\tau sij} = \int_{\Omega_e} (E_{\tau,zs} N_j N_i) d\Omega, \quad K_{u\sigma yz}^{\tau sij} = \int_{\Omega_e} (\hat{Z}_{pn23}^{\tau s} N_j N_{i,y} + \hat{Z}_{pn33}^{\tau s} N_j N_{i,x}) d\Omega,$$

$$K_{u\sigma zx}^{\tau sij} = \int_{\Omega_e} (E_{\tau s} N_j N_{i,x}) d\Omega, \quad K_{u\sigma zy}^{\tau sij} = \int_{\Omega_e} (E_{\tau s} N_j N_{i,y}) d\Omega, \quad K_{u\sigma zz}^{\tau sij} = \int_{\Omega_e} (E_{\tau,zs} N_j N_i) d\Omega, \quad (38)$$

$$K_{u\sigma xy}^{\tau sij} = 0, \quad K_{u\sigma yx}^{\tau sij} = 0,$$

$$K_{\sigma uxx}^{\tau sij} = \int_{\Omega_e} (E_{\tau s, \tau z} N_j N_i) d\Omega, \quad K_{\sigma uyz}^{\tau sij} = \int_{\Omega_e} (E_{\tau s} N_{j,x} N_i) d\Omega, \quad K_{\sigma uuy}^{\tau sij} = \int_{\Omega_e} (E_{\tau s, \tau y} N_j N_i) d\Omega,$$

$$K_{\sigma uyz}^{\tau sij} = \int_{\Omega_e} (E_{\tau s} N_{j,y} N_i) d\Omega, \quad K_{\sigma uzx}^{\tau sij} = - \int_{\Omega_e} (\hat{Z}_{np31}^{\tau s} N_{j,x} N_i - \hat{Z}_{np33}^{\tau s} N_{j,y} N_i) d\Omega,$$

$$K_{\sigma uzy}^{\tau sij} = - \int_{\Omega_e} (\hat{Z}_{np32}^{\tau s} N_{j,y} N_i - \hat{Z}_{np33}^{\tau s} N_{j,x} N_i) d\Omega, \quad K_{\sigma uzz}^{\tau sij} = \int_{\Omega_e} (E_{\tau s, \tau z} N_j N_i) d\Omega,$$

$$K_{\sigma uxy}^{\tau sij} = 0, \quad K_{\sigma uyx}^{\tau sij} = 0,$$

$$K_{\sigma\sigma xx}^{\tau sij} = - \int_{\Omega_e} (\hat{Z}_{nn11}^{\tau s} N_j N_i) d\Omega, \quad K_{\sigma\sigma xy}^{\tau sij} = - \int_{\Omega_e} (\hat{Z}_{nn12}^{\tau s} N_j N_i) d\Omega,$$

$$K_{\sigma\sigma yx}^{\tau sij} = - \int_{\Omega_e} (\hat{Z}_{nn21}^{\tau s} N_j N_i) d\Omega, \quad K_{\sigma\sigma xx}^{\tau sij} = - \int_{\Omega_e} (\hat{Z}_{nn22}^{\tau s} N_j N_i) d\Omega,$$

$$K_{\sigma\sigma xz}^{\tau sij} = 0, \quad K_{\sigma\sigma yz}^{\tau sij} = 0, \quad K_{\sigma\sigma zx}^{\tau sij} = 0, \quad K_{\sigma\sigma zy}^{\tau sij} = 0, \quad K_{\sigma\sigma zz}^{\tau sij} = 0.$$

The mass FN is the same as the PVD case; see Equation (33). Since the in-plane integrals are calculated via Gauss quadrature, it is crucial to consider an appropriate number of Gauss points in accordance with the variational rule of the fiber angle.

### 3. Results and Discussion

Three cases are analyzed in this work: a cantilever monolayer plate, a clamped multilayer plate and a clamped multilayer plate with a central circular cut-out. For each case, a square geometry is considered ($a = b = 1$ m). Parametric studies are performed considering different side-to-thickness ratios ($a/h = 100, 10, 5$). Material properties are represented in Table 1 for all the considered analyzed cases.

<table><caption>Table 1. Material properties.</caption><thead><tr><th>Case</th><th>E<sub>L</sub> (GPa)</th><th>E<sub>T</sub> (GPa)</th><th>G<sub>LT</sub> = G<sub>TT</sub> (GPa)</th><th>v<sub>LT</sub> = v<sub>TT</sub></th></tr></thead><tbody><tr><td>1</td><td>50.0</td><td>10.0</td><td>5.0</td><td>0.25</td></tr><tr><td>2</td><td>173.0</td><td>7.2</td><td>3.8</td><td>0.29</td></tr><tr><td>3</td><td>138.0</td><td>9.0</td><td>7.1</td><td>0.30</td></tr></tbody></table>

Reference solutions are represented by an Abaqus 3D model where quadratic elements (C3D20R) were used. For CUF solutions, nine-node square elements were used. For each case study, a preliminary convergence analysis was carried out to identify the appropriate mesh for both CUF and Abaqus solutions.
### 3.1. Monolayer Plate

The first case corresponds to a cantilever monolayer plate with density $ρ = 1540$ kg/m³. For this problem, axes $x'$ and $y'$ of the angle reference system are coincident with axes $x$ and $y$ of the plate. This means that the origin of the angle reference system is the same as the global one and that $x'$ and $y'$ are parallel to $x$ and $y$, respectively. It is assumed that the fiber angle is a linear function of $y'$; see Equation (10). The length parameter corresponds to $d = b$, while the direction of fiber variation $α$ corresponds to $y'$, which means that $Φ = 90°$. In this case, $T_0 = 0°$ and $T_1 = 90°$. The fiber orientation changes only along $y'$ from a value of $θ(0) = Φ + T_0 = 90°$ to $θ(b) = Φ + T_1 = 180°$. The angle variational law in this case can be expressed as $90° < 0, 90°>$, and it is presented in Figure 4.

<figure><img src="image_8.png" /><figcaption>Figure 4. Stacking sequence; case 1.</figcaption></figure>

This law has been taken from Viglietti et al. [27]. The reference solution contains 80 elements along each in-plane side and 12 elements along the thickness. The only clamped side of the plate is the one that lies on the $xz$ plane, corresponding to $y' = 0$. For the CUF results, a 10 × 10 mesh is considered. Table 2 shows the degrees of freedom (DOF) for some considered solutions.

<table><caption>Table 2. Degrees of freedom; case 1.</caption><thead><tr><td>Model</td><td>DOF</td></tr></thead><tbody><tr><td>Abaqus 3D</td><td>997,515</td></tr><tr><td>3LM4</td><td>34,398</td></tr><tr><td>2LM2</td><td>13,230</td></tr><tr><td>3LD4</td><td>17,199</td></tr><tr><td>2LD2</td><td>6615</td></tr><tr><td>ED4</td><td>6615</td></tr><tr><td>ED2</td><td>3969</td></tr><tr><td>FSDT</td><td>2646</td></tr><tr><td>CLT</td><td>2646</td></tr></tbody></table>

It is possible to observe that higher-order CUF models allow for a DOF reduction of one order of magnitude in comparison with the Abaqus 3D reference solution. Table 3 shows the first five natural frequencies for $a/h = 100$.
<table><caption>**Table 3.** Natural frequencies (Hz), *a/h* = 100; case 1.</caption><thead><tr><th></th><th>1</th><th>2</th><th>Mode<br>3</th><th>4</th><th>5</th></tr></thead><tbody><tr><td>Abaqus 3D</td><td>7.397</td><td>16.354</td><td>37.158</td><td>48.025</td><td>63.349</td></tr><tr><td>3LM4</td><td>7.399</td><td>16.334</td><td>37.164</td><td>47.988</td><td>63.310</td></tr><tr><td>2LM2</td><td>7.398</td><td>16.333</td><td>37.162</td><td>47.986</td><td>63.309</td></tr><tr><td>3LD4</td><td>7.400</td><td>16.362</td><td>37.179</td><td>48.053</td><td>63.378</td></tr><tr><td>2LD2</td><td>7.400</td><td>16.362</td><td>37.179</td><td>48.054</td><td>63.379</td></tr><tr><td>ED4</td><td>7.400</td><td>16.362</td><td>37.179</td><td>48.053</td><td>63.378</td></tr><tr><td>ED2</td><td>7.401</td><td>16.368</td><td>37.186</td><td>48.069</td><td>63.399</td></tr><tr><td>FSDT</td><td>7.398</td><td>16.363</td><td>37.171</td><td>48.054</td><td>63.388</td></tr><tr><td>CLT</td><td>7.403</td><td>16.414</td><td>37.213</td><td>48.175</td><td>63.537</td></tr></tbody></table>

For this case, classic and higher-order theories show very good approximations of the reference solution, where the maximum difference from the reference solution is 0.4% for the second natural frequency computed via CLT. Table 4 shows the first five natural frequencies for *a/h* = 10.

<table><caption>**Table 4.** Natural frequencies (Hz), *a/h* = 10; case 1.</caption><thead><tr><th></th><th>1</th><th>2</th><th>Mode<br>3</th><th>4</th><th>5</th></tr></thead><tbody><tr><td>Abaqus 3D</td><td>72.229</td><td>151.762</td><td>338.517</td><td>389.336</td><td>431.011</td></tr><tr><td>3LM4</td><td>72.244</td><td>151.751</td><td>338.577</td><td>389.554</td><td>431.004</td></tr><tr><td>2LM2</td><td>72.233</td><td>151.705</td><td>338.432</td><td>389.546</td><td>430.824</td></tr><tr><td>3LD4</td><td>72.250</td><td>151.796</td><td>338.625</td><td>389.587</td><td>431.151</td></tr><tr><td>2LD2</td><td>72.269</td><td>151.906</td><td>338.939</td><td>389.589</td><td>431.577</td></tr><tr><td>ED4</td><td>72.253</td><td>151.810</td><td>338.669</td><td>389.588</td><td>431.207</td></tr><tr><td>ED2</td><td>72.466</td><td>153.069</td><td>342.179</td><td>389.592</td><td>435.990</td></tr><tr><td>FSDT</td><td>72.437</td><td>153.021</td><td>342.036</td><td>389.510</td><td>435.853</td></tr><tr><td>CLT</td><td>73.825</td><td>163.064</td><td>365.813</td><td>389.510</td><td>472.565</td></tr></tbody></table>

It is possible to observe that classical and lower-order ESL theories are now less accurate, especially for the prediction of higher frequencies. For example, CLT, FSDT and ED2 models, corresponding to the third natural frequency, present a percentage error equal to 8.1%, 1.0% and 1.1%, respectively. This can be explained by considering that the side-to-thickness ratio *a/h* = 10 corresponds to a thick plate. In this case, higher-order theories are needed to obtain an accurate approximation. Since a moderately thick plate is considered, transverse shear stresses affect the solution. This is the reason that CLT, which neglects those stresses, is less close to the reference solution. The best approximations of plate natural frequencies are given by 2LM2 and 3LM4 mixed theories, which show a maximum percentage error of 0.1% each for the fourth natural frequency. In particular, it is possible to observe that the 2LM2 solution is globally closer to Abaqus in comparison with the 3LD4 solution, even if the last one is characterized by a higher number of degrees of freedom. Table 5 shows the first five natural frequencies for *a/h* = 5.

Because of the low side-to-thickness ratio, a very thick plate is considered, and lower-order theories do not provide a correct prediction of the natural frequencies. For CLT, the sixth mode is the same as the fifth mode of the reference solution, that is, the order of appearance is swapped. In this regard, mode tracking was performed by visually comparing the modes of each proposed solution with those of the reference solution obtained in Abaqus. The corresponding percentage error is as high as 27.1%. On the other hand, a 3LM4 model matches the Abaqus reference results.
<table><caption>**Table 5.** Natural frequencies (Hz), *a/h* = 5; case 1.</caption><thead><tr><th></th><th>1</th><th>2</th><th>Mode<br>3</th><th>4</th><th>5</th></tr></thead><tbody><tr><td>Abaqus 3D</td><td>136.723</td><td>264.080</td><td>389.391</td><td>556.394</td><td>704.284</td></tr><tr><td>3LM4</td><td>136.742</td><td>264.077</td><td>389.557</td><td>556.404</td><td>704.295</td></tr><tr><td>2LM2</td><td>136.667</td><td>263.747</td><td>389.550</td><td>555.332</td><td>703.121</td></tr><tr><td>3LD4</td><td>136.755</td><td>264.119</td><td>389.638</td><td>556.511</td><td>704.442</td></tr><tr><td>2LD2</td><td>136.875</td><td>264.684</td><td>389.643</td><td>558.145</td><td>706.381</td></tr><tr><td>ED4</td><td>136.774</td><td>264.224</td><td>389.640</td><td>556.855</td><td>704.832</td></tr><tr><td>ED2</td><td>138.015</td><td>269.553</td><td>389.651</td><td>570.563</td><td>721.354</td></tr><tr><td>FSDT</td><td>137.947</td><td>269.463</td><td>389.510</td><td>570.329</td><td>721.159</td></tr><tr><td>CLT</td><td>146.479</td><td>319.929</td><td>389.510</td><td>696.687</td><td>895.089</td></tr></tbody></table>

## 3.2. Multilayer Plate

The second case is taken from Viglietti et al. [27] and corresponds to a multilayer clamped plate with density $ρ = 1540$ kg/m³. The plate is composed of three layers with the same thicknesses. It is assumed that fiber angle is a function of $y'$ only, which means that $α$ is parallel to $y'$ ($Φ = 90°$). In this case, a linear law is considered for the fiber path, according to Equation (10). For this problem, axes $x'$ and $y'$ of the angle reference system are aligned with axes $x$ and $y$ of the plate, but their origin is placed on the center of the plate ($a/2, b/2$). In this case, $d = b/2$ is considered as the length parameter in Equation (10). $T_0$ and $T_1$ are set for each layer as follows: $T_0^{layer1} = T_0^{layer3} = 0°$, $T_0^{layer2} = -45°$, $T_1^{layer1} = T_1^{layer3} = 45°$, $T_1^{layer2} = -60°$. The lamination of the plate is $90° < 0, 45° >$ for layer 1, $90° < -45, -60° >$ for layer 2 and $90° < 0, 45° >$ for layer 3. The stacking sequence is presented in Figure 5.

<figure><img src="image_9.png" /><figcaption>**Figure 5.** Stacking sequence; case 2.</figcaption></figure>

As for the previous case, the Abaqus reference solution contains 80 elements along each side and 12 elements along the thickness. For the CUF results, a 10 × 10 mesh is considered. Table 6 shows the first five natural frequencies for $a/h = 100$, together with the results presented in Viglietti et al. [27].

In this case, the best approximation is given by the LM2 and LM4 theories. The LM2 and LM4 models both have a maximum percentage error as high as 0.4% corresponding to the third frequency. In addition, classical and low-order theories provide good results since a thin plate is considered. For this reason, transverse stresses do not play an important role. For example, the maximum error given by CLT is 2.1% for the fifth frequency. The case for $a/h = 10$ is presented in Table 7.
<table><caption>**Table 6.** Natural frequencies (Hz), *a/h* = 100; case 2.</caption><thead><tr><th></th><th>1</th><th>2</th><th>Mode<br>3</th><th>4</th><th>5</th></tr></thead><tbody><tr><td>Abaqus 3D</td><td>92.18</td><td>130.68</td><td>194.96</td><td>237.56</td><td>274.60</td></tr><tr><td>Ref. [27]</td><td>92.90</td><td>132.28</td><td>198.97</td><td>240.46</td><td>278.75</td></tr><tr><td>LM4</td><td>92.35</td><td>131.01</td><td>195.77</td><td>238.25</td><td>275.60</td></tr><tr><td>LM2</td><td>92.34</td><td>130.99</td><td>195.74</td><td>238.23</td><td>275.58</td></tr><tr><td>LD4</td><td>92.36</td><td>131.03</td><td>195.81</td><td>238.30</td><td>275.67</td></tr><tr><td>LD2</td><td>92.36</td><td>131.04</td><td>195.84</td><td>238.31</td><td>275.69</td></tr><tr><td>ED4</td><td>92.37</td><td>131.06</td><td>195.88</td><td>238.32</td><td>275.72</td></tr><tr><td>ED2</td><td>92.49</td><td>131.23</td><td>196.16</td><td>238.97</td><td>276.48</td></tr><tr><td>FSDT</td><td>92.38</td><td>131.01</td><td>195.75</td><td>238.74</td><td>276.20</td></tr><tr><td>CLT</td><td>93.04</td><td>131.85</td><td>197.00</td><td>242.48</td><td>280.40</td></tr></tbody></table>

<table><caption>**Table 7.** Natural frequencies (Hz), *a/h* = 10; case 2.</caption><thead><tr><th></th><th>1</th><th>2</th><th>Mode<br>3</th><th>4</th><th>5</th></tr></thead><tbody><tr><td>Abaqus 3D</td><td>606.67</td><td>896.70</td><td>1208.24</td><td>1313.26</td><td>1458.25</td></tr><tr><td>Ref. [27]</td><td>609.79</td><td>903.63</td><td>1216.00</td><td>1328.41</td><td>1469.33</td></tr><tr><td>LM4</td><td>606.90</td><td>897.26</td><td>1208.80</td><td>1314.85</td><td>1459.23</td></tr><tr><td>LM2</td><td>606.33</td><td>896.52</td><td>1206.86</td><td>1313.56</td><td>1457.30</td></tr><tr><td>LD4</td><td>607.22</td><td>897.73</td><td>1209.64</td><td>1315.80</td><td>1460.16</td></tr><tr><td>LD2</td><td>608.65</td><td>901.20</td><td>1213.06</td><td>1322.93</td><td>1465.20</td></tr><tr><td>ED4</td><td>609.84</td><td>905.18</td><td>1214.60</td><td>1331.82</td><td>1469.17</td></tr><tr><td>ED2</td><td>633.68</td><td>941.96</td><td>1272.39</td><td>1396.16</td><td>1540.10</td></tr><tr><td>FSDT</td><td>632.82</td><td>940.46</td><td>1271.42</td><td>1393.96</td><td>1538.74</td></tr><tr><td>CLT</td><td>921.28</td><td>1287.71</td><td>2368.22</td><td>1885.61</td><td>2699.22</td></tr></tbody></table>

Here, the CLT model shows that the inversion of the third and fourth modes can be observed by the corresponding values of the frequencies that are not in ascending order as the mode number increases. In comparison with the monolayer plate, in this case, the mode inversions of the CLT model can be seen for higher side-to-thickness ratios and lower frequencies. For the third mode, CLT shows a percentage error of 96.0%, while the best approximation is given by LM4, which has a percentage error of 0.17% for the same mode. Table 8 shows the first five frequencies for *a/h* = 5.

<table><caption>**Table 8.** Natural frequencies (Hz), *a/h* = 5; case 2.</caption><thead><tr><th></th><th>1</th><th>2</th><th>Mode<br>3</th><th>4</th><th>5</th></tr></thead><tbody><tr><td>Abaqus 3D</td><td>794.730</td><td>1201.916</td><td>1439.956</td><td>1701.328</td><td>1810.250</td></tr><tr><td>LM4</td><td>794.760</td><td>1202.101</td><td>1440.092</td><td>1701.788</td><td>1811.113</td></tr><tr><td>LM2</td><td>792.734</td><td>1199.331</td><td>1433.897</td><td>1696.266</td><td>1805.942</td></tr><tr><td>LD4</td><td>795.213</td><td>1202.777</td><td>1441.080</td><td>1702.986</td><td>1812.317</td></tr><tr><td>LD2</td><td>799.063</td><td>1209.706</td><td>1448.714</td><td>1713.982</td><td>1820.716</td></tr><tr><td>ED4</td><td>802.019</td><td>1216.744</td><td>1450.930</td><td>1723.900</td><td>1825.405</td></tr><tr><td>ED2</td><td>845.154</td><td>1294.481</td><td>1523.246</td><td>1847.193</td><td>1930.364</td></tr><tr><td>FSDT</td><td>844.048</td><td>1292.846</td><td>1522.478</td><td>1845.945</td><td>1928.631</td></tr><tr><td>CLT</td><td>1790.121</td><td>2411.198</td><td>-</td><td>-</td><td>-</td></tr></tbody></table>

In this case, lower-order theories have an evident loss of accuracy. The CLT model can predict only the first two modes. In addition, the FSDT and ED2 models show non-negligible errors, which become bigger with the increase in frequency. On the other hand, mixed models are able to correctly predict the dynamic behavior of the plate for both low and high frequencies.
### 3.3. Multilayer Plate with Central Hole

Case 3 is taken from Hachemi et al. [9] and corresponds to a multilayer clamped plate that presents a circular cut-out. The center of the cut-out is placed at the plate center (a/2, b/2), and its radius is r = 0.2 m. It is assumed that the fiber angle is a parabolic function of x', which means that a is parallel to x' (Φ = 0°). As in the previous case, the x' and y' axes are parallel, respectively, to x and y, and their origin is placed at the center of the plate. The angle variational law is defined in Equation (12), considering d = a/2. The plate is composed of two layers that have the same thicknesses. The values of T₀ and T₁ are set for each layer as follows: T₀<sup>layer1</sup> = T₀<sup>layer2</sup> = 0°, T₁<sup>layer1</sup> = 30°, T₁<sup>layer2</sup> = -30°. The stacking sequence is 0 < 0, ±30 >; see Figure 6.

<figure><img src="image_10.png" /><figcaption>Figure 6. Stacking sequence; case 3.</figcaption></figure>

In this case, the Abaqus reference solution is made of 73728 elements: 4608 elements are defined on the plate plane, and 16 elements are defined along the thickness. For the CUF results, 128 elements are used on the plate plane. The natural frequencies are expressed in the following dimensionless form:

$$ \bar{\omega} = (\omega a^2) \sqrt{\rho h / D_0}, \qquad (39) $$

$$ D_0 = E_2 h^3 / 12 (1 - v_{LT} v_{TL}), \qquad (40) $$

where ω is the natural frequency, while D₀ represents a reference bending stiffness. Table 9 presents the first five non-dimensional frequencies for a/h = 100.

<table><caption>Table 9. Non-dimensional frequencies $\bar{\omega}$, a/h = 100; case 3.</caption><thead><tr><th></th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th></tr></thead><tbody><tr><td>Abaqus 3D</td><td>87.079</td><td>106.407</td><td>147.559</td><td>184.034</td><td>197.096</td></tr><tr><td>LM4</td><td>87.281</td><td>106.622</td><td>147.070</td><td>184.554</td><td>197.522</td></tr><tr><td>LM2</td><td>87.259</td><td>106.593</td><td>147.045</td><td>184.500</td><td>197.489</td></tr><tr><td>LD4</td><td>87.327</td><td>106.704</td><td>147.911</td><td>184.789</td><td>197.969</td></tr><tr><td>LD2</td><td>87.336</td><td>106.719</td><td>147.952</td><td>184.821</td><td>198.022</td></tr><tr><td>ED4</td><td>87.331</td><td>106.708</td><td>147.921</td><td>184.798</td><td>197.984</td></tr><tr><td>ED2</td><td>87.364</td><td>106.768</td><td>148.169</td><td>184.931</td><td>198.228</td></tr><tr><td>FSDT</td><td>87.184</td><td>106.538</td><td>148.047</td><td>184.525</td><td>198.029</td></tr><tr><td>CLT</td><td>87.387</td><td>106.942</td><td>150.080</td><td>185.420</td><td>199.725</td></tr></tbody></table>

It is possible to observe that the theories show a good approximation of the reference results. In addition, the percentage errors of FSDT and CLT are less than 2%. Mixed theories match the Abaqus results. Table 10 shows the results for a/h = 10 in order to compare the
Abaqus reference solution with the one presented in Hachemi et al. [9] and the solutions obtained with CUF.

<table><caption>**Table 10.** Non-dimensional frequencies $\bar{\omega}$, $a/h = 10$; case 3.</caption><thead><tr><th rowspan="2"></th><th colspan="5">Mode</th></tr><tr><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th></tr></thead><tbody><tr><td>Abaqus 3D</td><td>72.645</td><td>86.745</td><td>104.279</td><td>136.366</td><td>140.278</td></tr><tr><td>Ref. [9]</td><td>72.432</td><td>86.626</td><td>103.910</td><td>135.828</td><td>139.747</td></tr><tr><td>LM4</td><td>72.699</td><td>86.830</td><td>104.307</td><td>136.467</td><td>140.408</td></tr><tr><td>LM2</td><td>72.573</td><td>86.700</td><td>104.051</td><td>136.137</td><td>140.143</td></tr><tr><td>LD4</td><td>72.744</td><td>86.888</td><td>104.376</td><td>136.558</td><td>140.516</td></tr><tr><td>LD2</td><td>73.107</td><td>87.263</td><td>105.144</td><td>137.567</td><td>141.231</td></tr><tr><td>ED4</td><td>72.868</td><td>86.990</td><td>104.630</td><td>136.851</td><td>140.725</td></tr><tr><td>ED2</td><td>73.977</td><td>88.609</td><td>107.143</td><td>140.522</td><td>143.556</td></tr><tr><td>FSDT</td><td>74.075</td><td>88.782</td><td>107.645</td><td>141.221</td><td>143.885</td></tr><tr><td>CLT</td><td>84.751</td><td>104.166</td><td>143.133</td><td>190.321</td><td>174.656</td></tr></tbody></table>

As already observed in previous cases, classical theories and, in general, low-order ones are not able to provide an accurate approximation of natural frequencies, because of the low side-to-thickness ratio value. It is also possible that this generates the inversion of modes four and five for the CLT case. On the other hand, the best approximation is given by mixed theories, which are closer to the Abaqus solution for high frequencies. The shapes of the modes are presented in Figures 7–11 for $a/h = 10$. The modal shapes obtained with the LM4 model are compared with those of Abaqus 3D.

<figure><img src="image_11.png" /><figcaption>**Figure 7.** Mode 1, $a/h = 10$; case 3.</figcaption></figure>

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 8.** Mode 2, $a/h = 10$; case 3.</figcaption></figure>

<figure><img src="image_14.png" /></figure>
<figure><img src="image_15.png" /><figcaption>Figure 9. Mode 3, *a/h* = 10; case 3.</figcaption></figure>

<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /><figcaption>Figure 10. Mode 4, *a/h* = 10; case 3.</figcaption></figure>

<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /><figcaption>Figure 11. Mode 5, *a/h* = 10; case 3.</figcaption></figure>

<figure><img src="image_20.png" /></figure>

The first mode shows a simple bending of the plate on the *xy* plane with a single half-wave along each in-plane direction. The second and the third modes show two half-waves in the *y* and *x* directions, respectively. Mode number four shows three half-waves along the plate diagonally between the *x* and *y* axes. The fifth mode shows three half-waves along the *y* direction. Finally, Table 11 shows the frequencies for *a/h* = 5.
<table><caption>**Table 11.** Non-dimensional frequencies $\bar{\omega}$, a/h = 5; case 3.</caption><thead><tr><th rowspan="2"></th><th colspan="5">Mode</th></tr><tr><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th></tr></thead><tbody><tr><td>Abaqus 3D</td><td>54.333</td><td>64.456</td><td>70.572</td><td>90.875</td><td>98.086</td></tr><tr><td>LM4</td><td>54.326</td><td>64.456</td><td>70.541</td><td>90.866</td><td>98.098</td></tr><tr><td>LM2</td><td>54.038</td><td>64.201</td><td>70.036</td><td>90.292</td><td>97.612</td></tr><tr><td>LD4</td><td>54.388</td><td>64.514</td><td>70.619</td><td>90.956</td><td>98.191</td></tr><tr><td>LD2</td><td>54.875</td><td>64.963</td><td>71.421</td><td>91.868</td><td>98.955</td></tr><tr><td>ED4</td><td>54.554</td><td>64.623</td><td>70.913</td><td>91.224</td><td>98.408</td></tr><tr><td>ED2</td><td>56.062</td><td>66.756</td><td>73.181</td><td>94.442</td><td>101.535</td></tr><tr><td>FSDT</td><td>56.253</td><td>66.985</td><td>73.702</td><td>95.219</td><td>102.017</td></tr><tr><td>CLT</td><td>76.928</td><td>95.975</td><td>119.513</td><td>-</td><td>-</td></tr></tbody></table>

Since a thick plate is considered, the effect of transverse stresses is not negligible, which causes the classical and lower-order theories to be inaccurate. This can be observed for CLT, which is not able to predict the fourth and fifth modes and has an error as high as 69.4% for the third mode. Considering the FSDT, ED4 and LD4 models, this error can be reduced to 4.4%, 0.5% and 0.1%, respectively.

## 4. Conclusions

In this paper, a new framework for the dynamic analysis of VAT structures is presented. RMVT is developed within CUF in order to obtain a new family of 2D models for the free-vibration analysis of VAT plates. The results are obtained via either RMVT or PVD and are compared in order to show the effective capabilities of the proposed method in the prediction of VAT plates' natural frequencies. The Abaqus 3D reference solutions and results from Refs. [9,27] are also included to further validate the models proposed in this article. Linear and parabolic laws are both considered in order to describe the in-plane path of fiber variation. The possibility to use a polynomial order defined a priori through CUF and the introduction of the transverse stress field as a primary variable of the problem through RMVT both help to obtain a valid approach for the prediction of VAT dynamic behavior. After the results analysis, the following remarks can be made:

* Classical theories (FSDT and CLT) provide the best trade-off between accuracy and computational costs for thin plates ($a/h = 100$), whereas they are not able to correctly predict the behavior of thicker plates ($a/h = 10$ and 5), specially at high frequencies. The loss of accuracy is more evident for CLT results, since this theory does not consider transverse shear stresses, which become important in thick plates. This error is particularly evident in the second- and third-order theories, where the inversion of modes can be observed.

* The PVD results show monotonic convergence to the reference solution: the lower the DOF number, the higher the frequency value. For a given mode, frequency values decrease when higher-order models are employed, and they move closer to the reference solution.

* In all the cases, layer-wise mixed theories yield the best match of the reference 3D solution, independently from the plate geometry or fiber variational law. This is justified by the fact that RMVT considers both displacements and transverse stresses as primary variables, assuring a better approximation of the transverse stresses field into the problem domain, and improving the overall solution accuracy.

* For a given expansion order, models based on RMVT are more computationally expensive than PVD models. For this reason, the use of LW mixed models is advantageous in the cases where a more precise representation of the through-the-thickness behavior is needed, as in the case of higher frequencies or for thick plates, whereas low-order ESL and classical models are accurate for lower frequencies and thin plates.

In conclusion, the application of RMVT within CUF has demonstrated significant potential for improving the accuracy and efficiency of modeling VAT plates for free-vibration
analyses. The promising results suggest, as future perspectives, the extension to buckling and failure analyses where an accurate and efficient modeling of VAT structures under various loading and operational conditions is required.

**Author Contributions:** Methodology, G.G., D.A.I. and M.M.; Software, D.A.I.; Validation, D.A.I.; Writing—original draft, D.A.I.; Writing—review & editing, G.G. and M.M.; Supervision, G.G. and M.M.; Funding acquisition, G.G. and M.M. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Luxembourg National Research Fund (FNR), grant reference INTER/ANR/21/16215936 GLAMOUR-VSC. For the purpose of open access, and in fulfillment of the obligations arising from the grant agreement, the author has applied a Creative Commons Attribution 4.0 International (CC BY 4.0) license to any Author Accepted Manuscript version arising from this submission. M. Montemurro is grateful to French National Research Agency for supporting this work through the research project GLAMOUR-VSC (Global–LocAl two-level Multi-scale optimisation strategy accOUnting for pRocess-induced singularities to design Variable Stiffness Composites) ANR-21-CE10-0014.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

## References

1. Dirk, H.J.L.; Ward, C.; Potter, K.D. The engineering aspects of automated prepreg layup: History, present and future. *Compos. Part B Eng.* **2012**, *43*, 997–1009.

2. Zhuo, P.; Li, S.; Ashcroft, I.A.; Jones, A.I. Material extrusion additive manufacturing of continuous fiber reinforced polymer matrix composites: A review and outlook. *Compos. Part B Eng.* **2021**, *224*, 109143. [CrossRef]

3. Brooks, T.R.; Martins, J.R.; Kennedy, G.J. High-fidelity aerostructural optimization of tow-steered composite wings. *J. Fluids Struct.* **2019**, *88*, 122–147. [CrossRef]

4. Grenoble, R.W.; Nguyen, T.; McKenney, M.J.; Przekop, A.; Juarez, P.D.; Gregory, E.D.; Jegley, D.C. Fabrication of a composite tow-steered structure for air-launch vehicle applications. In Proceedings of the AIAA/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, Kissimmee, FL, USA, 8–12 January 2018.

5. Hyer, M.W.; Charette, R.F. *Innovative Design of Composite Structures: The Use of Curvilinear Fiber Format in Composite Structure Design*; NASA Langley Research Center: Hampton, VA, USA, 1990.

6. Hyer, M.W.; Lee, H.H. Innovative design of composite structures: The use of curvilinear fiber format to improve buckling resistance of composite plates with central circular holes. *Compos. Struct.* **1991**, *18*, 239–261. [CrossRef]

7. Akhavan, H.; Ribeiro, P. Natural modes of vibration of variable stiffness composite laminates with curvilinear fibers. *Compos. Struct.* **2011**, *93*, 3040–3047. [CrossRef]

8. Ribeiro, P.; Akhavan, H. Non-linear vibrations of variable stiffness composite laminated plates. *Compos. Struct.* **2012**, *94*, 2424–2432. [CrossRef]

9. Hachemi, M.; Hamza-Cherif, S.M.; Houmat, A. Free vibration analysis of variable stiffness composite laminate plate with circular cutout. *Aust. J. Mech. Eng.* **2020**, *18*, 63–79. [CrossRef]

10. Zhao, W.; Kapania, R.K. Prestressed vibration of stiffened variable-angle tow laminated plates. *AIAA J.* **2019**, *57*, 2575–2593. [CrossRef]

11. Honda, S.; Narita, Y. Natural frequencies and vibration modes of laminated composite plates reinforced with arbitrary curvilinear fiber shape paths. *J. Sound Vib.* **2012**, *331*, 180–191. [CrossRef]

12. Rodrigues, J.D.; Ribeiro, P.; Akhavan, H. Experimental and finite element modal analysis of variable stiffness composite laminated plates. In Proceedings of the 11th Biennial International Conference on Vibration Problems (ICOVP-2013), Lisbon, Portugal, 9–12 September 2013; Volume 30.

13. Stodieck, O.; Cooper, J.E.; Weaver, P.M.; Kealy, P. Improved aeroelastic tailoring using tow-steered composites. *Compos. Struct.* **2013**, *106*, 703–715. [CrossRef]

14. Abdalla, M.M.; Setoodeh, S.; Gürdal, Z. Design of variable stiffness composite panels for maximum fundamental frequency using lamination parameters. *Compos. Struct.* **2007**, *81*, 283–291. [CrossRef]

15. Blom, A.W.; Setoodeh, S.; Hol, J.; Gürdal, Z. Design of variable-stiffness conical shells for maximum fundamental eigenfrequency. *Comput. Struct.* **2008**, *86*, 870–878. [CrossRef]

16. Carvalho, J.; Sohouli, A.; Suleman, A. Fundamental Frequency Optimization of Variable Angle Tow Laminates with Embedded Gap Defects. *J. Compos. Sci.* **2022**, *6*, 64. [CrossRef]
17. Montemurro, M.; Catapano, A. On the effective integration of manufacturability constraints within the multi-scale methodology for designing variable angle-tow laminates. *Compos. Struct.* **2017**, *161*, 145–159. [CrossRef]

18. Catapano, A.; Montemurro, M.; Balcou, J.-A.; Panettieri, E. Rapid prototyping of variable angle-tow composites. *Aerotec. Missili Spaz.* **2019**, *98*, 257–271. [CrossRef]

19. Montemurro, M.; Catapano, A. A general B-Spline surfaces theoretical framework for optimisation of variable angle-tow laminates. *Compos. Struct.* **2019**, *209*, 561–578. [CrossRef]

20. Fiordilino, G.A.; Izzi, M.I.; Montemurro, M. A general isogeometric polar approach for the optimisation of variable stiffness composites: Application to eigenvalue buckling problems. *Mech. Mater.* **2021**, *153*, 103574. [CrossRef]

21. Carrera, E. Theories and finite elements for multilayered, anisotropic, composite plates and shells. *Arch. Comput. Meth. Eng.* **2002**, *9*, 87–140. [CrossRef]

22. Carrera, E. Theories and finite elements for multilayered plates and shells: A unified compact formulation with numerical assessment and benchmarking. *Arch. Comput. Meth. Eng.* **2003**, *10*, 215–296. [CrossRef]

23. Carrera, E.; Giunta, G.; Brischetto, S. Hierarchical closed form solutions for plates bent by localized transverse loadings. *J. Zhejiang Univ. Sci. A* **2007**, *8*, 1026–1037. [CrossRef]

24. Carrera, E.; Giunta, G. Hierarchical models for failure analysis of plates bent by distributed and localized transverse loadings. *J. Zhejiang Univ. Sci. A* **2008**, *9*, 600–613. [CrossRef]

25. Giunta, G.; Catapano, A.; Belouettar, S. Failure indentation analysis of composite sandwich plates via hierarchical models. *J. Sandw. Struct. Mater.* **2013**, *15*, 45–70. [CrossRef]

26. Giunta, G.; Biscani, F.; Belouettar, S.; Ferreira, A.J.M.; Carrera, E. Free vibration analysis of composite beams via refined theories. *Compos. Part B Eng.* **2013**, *44*, 540–552. [CrossRef]

27. Viglietti, A.; Zappino, E.; Carrera, E. Analysis of variable angle tow composites structures using variable kinematic models. *Compos. Part B Eng.* **2019**, *171*, 272–283. [CrossRef]

28. Fallahi, N.; Viglietti, A.; Carrera, E.; Pagani, A.; Zappino, E. Effect of fiber orientation path on the buckling, free vibration, and static analyses of variable angle tow panels. *Facta Univ. Ser. Mech. Eng.* **2020**, *18*, 165–188. [CrossRef]

29. Sánchez-Majano, A.R.; Azzara, R.; Pagani, A.; Carrera, E. Accurate Stress Analysis of Variable Angle Tow Shells by High-Order Equivalent-Single-Layer and Layer-Wise Finite Element Models. *Materials* **2021**, *14*, 6486. [CrossRef]

30. Pagani, A.; Sánchez-Majano, A.R. Influence of fiber misalignments on buckling performance of variable stiffness composites using layerwise models and random fields. *Mech. Adv. Mater. Struct.* **2022**, *29*, 384–399. [CrossRef]

31. Pagani, A.; Sánchez-Majano, A.R. Stochastic stress analysis and failure onset of variable angle tow laminates affected by spatial fiber variations. *Compos. Part C Open Access* **2021**, *4*, 100091. [CrossRef]

32. Sánchez-Majano, A.R.; Pagani, A.; Petrolo, M.; Zhang, C. Buckling sensitivity of tow-steered plates subjected to multiscale defects by high-order finite elements and polynomial chaos expansion. *Materials* **2021**, *14*, 2706. [CrossRef]

33. Vescovini, R.; Dozio, L. A variable-kinematic model for variable stiffness plates: Vibration and buckling analysis. *Compos. Struct.* **2016**, *142*, 15–26. [CrossRef]

34. Demasi, L.; Biagini, G.; Vannucci, F.; Santarpia, E.; Cavallaro, R. Equivalent Single Layer, Zig-Zag, and Layer Wise theories for variable angle tow composites based on the Generalized Unified Formulation. *Compos. Struct.* **2017**, *177*, 54–79. [CrossRef]

35. Carrera, E.; Demasi, L. Classical and advanced multilayered plate elements based upon PVD and RMVT. Part 1: Derivation of finite element matrices. *Int. J. Numer. Methods Eng.* **2002**, *55*, 191–231. [CrossRef]

36. Carrera, E.; Demasi, L. Classical and advanced multilayered plate elements based upon PVD and RMVT. Part 2: Numerical implementations. *Int. J. Numer. Methods Eng.* **2002**, *55*, 253–291. [CrossRef]

37. Babaei, M.; Kiarasi, F.; Tehrani, M.S.; Hamzei, A.; Mohtarami, E.; Asemi, K. Three dimensional free vibration analysis of functionally graded graphene reinforced composite laminated cylindrical panel. *Proc. Inst. Mech. Eng. Part L J. Mater. Des. Appl.* **2022**, *236*, 1501–1514. [CrossRef]

38. Reddy, J.N. *Mechanics of Laminated Composite Plates and Shells: Theory and Analysis*; CRC Press: Boca Raton, FL, USA, 2003.

39. Gürdal, Z.; Tatting, B.F.; Wu, C.K. Variable stiffness composite panels: Effects of stiffness variation on the in-plane and buckling response. *Compos. Part A Appl. Sci. Manuf.* **2008**, *39*, 911–922. [CrossRef]

40. Honda, S.; Oonishi, Y.; Narita, Y.; Sasaki, K. Vibration analysis of composite rectangular plates reinforced along curved lines. *J. Syst. Des. Dyn.* **2008**, *2*, 76–86. [CrossRef]

41. Carrera, E. On the use of the Murakami's zig-zag function in the modeling of layered plates and shells. *Comput. Struct.* **2004**, *82*, 541–554. [CrossRef]

42. Bathe, K.-J. *Finite Element Procedures*; Prentice Hall: Upper Saddle River, NJ, USA, 2006.

**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.