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

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.


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 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 a 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.

Carrera's Unified Formulation
A plate is a flat body whose material points lie in the Cartesian closed-point subset of the three-dimensional space R 3 where: where a and b are the dimensions along the two in-plane axes, and h measures its thickness along the z-axis, where z a and b. The global reference system and plate geometry are presented in Figure 1. The displacement field is expressed as: ( The strain tensor components can be written in vector form. Two vectors are obtained, representing the in-plane and out-of-plane components: The hypothesis of small displacements helps to use a linear strain-displacement relation: where D p , D nΩ and D nz are the following differential operators: The stress vector is expressed in a similar manner: Hooke's law reads: where the terms C pp , C pn , C np and C nn are subcomponents of a material stiffness matrix C according to the stress and strain ordering in Equations (4) and (7), where the fibers lay in Ω 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 ν LT and ν TT , where subscripts L and T stand for the directions parallel and perpendicular to the fibers, respectively. For further details, see Reddy [38].

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 C reads: Superscript T stands for the transpose operator. The matrix T represents a rotation matrix that depends on an in-plane rotation angle θ. For the sake of brevity, the components of C and T are not reported here; they can be found in Reddy [38]. In a laminated VAT, the rotation angle θ is a bi-dimensional field in Ω. In this work, two different variation laws are considered for θ, a linear variation law and a parabolic one. The linear law can be expressed according to the following formula: The angle Φ describes the original direction along which θ varies, and α is a generic spatial variable defined as: x and y denote a generic in-plane reference system used for describing a fiber path, where θ 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 α-axis and the tangent to a fiber for α equal to zero and d, respectively; see Figure 2. 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 α, defined as a combination of x and y depending on the angle Φ. 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: The parameter γ 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 (±γ). 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: Φ < T 0 , T 1 >.

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 σ n as a primary unknown, whereas the PVD considers only displacements as primary variables. For the PVD case, the following variational statement applies: 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). L in is the virtual work of the inertial forces, and δ stands for a virtual variation. For the RMVT case, the variational statement is: 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: whereĈ pp ,Ĉ pn ,Ĉ np andĈ nn are (see Carrera and Demasi [35]): The superscript "−1" indicates the inverse of a matrix. The inertial work can be expressed as: where ρ is the plate material density, andü represents the acceleration vector.

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: 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 τ is an approximation function in H, and g τ is an unknown two-dimensional function in Ω. 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 τ 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 τ . In an ESL model N l is the total number of laminae, and h k is the thickness of a generic k lamina such that 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 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: 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 ∞ -continuity over 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-thethickness 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: 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: are typically used as approximation functions over 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), ξ k = 2z k h k ∈ [−1, 1] and P i = P i ξ 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 ξ k = ±1, being, by definition, F r (±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.

Acronym System
An acronym system is used in order to identify all the derived theories. Figure 3 shows this system. 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: In vector form: where F τ = z τ and u τ = u τ (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: For LMN solutions, transverse stresses are treated as primary variables. The transverse stress field can be expressed as: 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.

FE Stiffness Matrices
As far as a FEM solution is concerned, the in-plane domain is discretized into N e subdomains such as Ω = N e e=1 Ω e and Ω e ∩ Ω e = ∅ for e = e . Shape functions are then introduced for the approximation of the variation over Ω e . In the case of a bi-dimensional model, Equation (18) becomes: 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: Through the substitution of Equations (5), (8) and (27) into Equation (13), the weak PVD form can be obtained: where: Z τs wr , Z τ, z s wr , Z τs, z wr , Z τ, z s, z wr = C wr E τs , C wr E τ, z s , C wr E τs, z , C wr E τ, z s, z : w, r = p, n , (29) An axis coordinate as comma-preceded subscript stands for a derivative in that coordinate direction. In compact vector form, Equation (28) reads: where K τsij and M τsij ∈ R 3×3 are fundamental nuclei (FN) of the stiffness and mass matrices, respectively. Through the cycles on the indices τ, 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: The mass FN can be written as: It is possible to observe that M τsij is a diagonal matrix and that since the plate density is assumed to be constant, the term ρE τs can be placed outside the integral.
In the RMVT case, transverse stresses are a priori approximated: Through the substitution of Equations (5), (15), (27) and (34) into Equation (14), the RMVT weak form can be obtained: where: Ẑ τs wr ,Ẑ τ, z s wr ,Ẑ τs, z wr ,Ẑ τ, z s, z wr = Ĉ wr E τs ,Ĉ wr E τ, z s ,Ĉ wr E τs, z ,Ĉ wr E τ, z s, z : w, r = p, n . (36) In a compact form: In this case, four fundamental nuclei are obtained. The components of the FN for the RMVT case can be written as: 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.

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.  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.

Monolayer Plate
The first case corresponds to a cantilever monolayer plate with density ρ = 1540 kg/m 3 . 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. 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. 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. 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. 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 lowerorder 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.

Multilayer Plate
The second case is taken from Viglietti et al. [27] and corresponds to a multilayer clamped plate with density ρ = 1540 kg/m 3 . 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: 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. 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. 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. 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 nonnegligible 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.

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 α 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 0 and T 1 are set for each layer as follows: The stacking sequence is 0 < 0, ±30 >; see Figure 6. 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: where ω is the natural frequency, while D 0 represents a reference bending stiffness. Table 9 presents the first five non-dimensional frequencies for a/h = 100. 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. 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.    The first mode shows a simple bending of the plate on the xy plane with a single halfwave 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. 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.

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 freevibration 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.