Geometrically nonlinear analysis and vibration of in-plane-loaded variable angle tow composite plates and shells

This work intends to present a novel numerical approach for studying the vibration behaviours of variable angle tow (VAT) composite structures in their quasi-static nonlinear equilibrium states. This methodology is able to predict the buckling load, to investigate the natural frequencies variation for progressively higher loads, and to provide a means for verifying experimental Vibration Correlation Technique results. The use of VAT composites, in which the fibre orientations are allowed to vary along with a curvilinear pattern within each lamina, dramatically increases the design space and provides a significant improvement in buckling performance and benefits in the postbuckling regime. This study has been performed using an innovative methodology based on the well-established Carrera Unified Formulation able to describe several kinematic models for two-dimensional structures. In detail, layerwise theories are employed to characterize the complex phenomena that may appear in VAT composite structures. All Green-Lagrange strain components are employed because far nonlinear regimes are investigated. Furthermore, the geometrical nonlinear equations are written in a total Lagrangian framework and solved with an opportune Newton–Raphson method along with a path-following approach based on the arc-length constraint. Different VAT composite structures have been analyzed to validate the proposed approach and provide some benchmark solutions. The computed equilibrium paths are compared with results obtained using the commercial code ABAQUS. The results document the good accuracy and reliability of the presented methodology and show this numerical tool’s potentialities.

advanced composite structures, called variable angle tow (VAT), have attracted a lot of interest because they allow the manufacturing of Variable Stiffness Composite Laminates (VSCL) without discontinuity in the material while maximizing the ratio of stiffness to mass.
The concept of VAT composites is not new, but it was proposed several decades ago. However, these new classes of materials have recently gained renewed interest thanks to significant technological improvements in the automatic manufacturing process [5]. In particular, the advanced automated fibre placements (AFP) techniques, the Continuous Tow Shearing (CTS) process, and the Automated Tape Laying (ATL), which allow the fibre orientation angle of a layer to vary as a function of one or more spatial directions, have led to developing accurate and computationally light mathematical models of VAT structures. For this reason, the study of VAT materials is still one of the most crucial aspects for scientists and researchers. In fact, the literature concerning the VAT theories is vast. For example, Hyer and Lee [6] adopted curvilinear fibres to change the stress concentration around a hole. Lopes et al. [7] conducted several progressive damage and failure analyses on classical and VAT structures with and without cutouts, proving the potentialities of the variable stiffness panels to redirect the load fluxes to the stiffer edge area to improve the structural performance of the structure. VAT techniques were used and compared with the traditional composite structures by Stodieck et al. [8] to investigate the possibility to improve the aeroelastic tailoring of a rectangular unswept laminate wing. Several works based on numerical and experimental analysis have demonstrated the advantages given by the application of VAT structure [9][10][11].
Among the numerous benefits introduced by the use of VAT structures, elastic stability is undoubtedly one of the crucial aspects in the design of plate and shell structures commonly used in aerospace applications. Indeed, the usage of VAT composite structures has significantly increased the buckling performance, which is often the primary design criterion of structures. Numerous researches have amply demonstrated the remarkable improvement in the buckling carrying capacity given to the use of VAT placement techniques [12][13][14]. These works generally relied on the classical weak-form finite element method (FEM) for analysis and design of VAT structures [10,15,16]. In addition, the pseudo-spectral differential quadrature method (DQM) [17][18][19] has also been widely adopted for buckling and postbuckling analyses of VAT plates and shells. Numerous articles were written about the buckling and the postbuckling of VAT two-dimensional (2D) structures, and some of these are presented here. Madeo et al. [20] performed a postbuckling analysis of VAT composite 2D structures by employing the Koiter's approach and the FEM. Zucco et al. [21] formulated a new mixed shell element to carry out static and buckling analysis of VAT plate structures. An efficient postbuckling analysis of VAT plates using a perturbation approach was performed by Vescovini et al. [22]. For brevity, the reader is referred to [23][24][25] for other important works on buckling and postbuckling of VATs.
The stability of structures subjected to compressive loads, essential for the design and validation of the safety of new structures, continues to be extensively verified through experimental tests. However, one of the researchers' goals over the years has been to reduce both the time and cost of operations of these complex analyses to determine the buckling load. One of the effective ways to do this is to adopt nondestructive experimental tests to evaluate the critical load of structures. The most used nondestructive method is represented by the Vibration Correlation Technique (VCT) [26,27]. This method calculates the buckling load and the equivalent boundary conditions by interpolating the natural frequencies of the structures for progressively increasing applied loads without reaching the point of instability. The first experimental VCT investigations were data from the 1950s, with the studies conducted by Lurie [28], Meier [29] and Chu [30]. In some cases, the relationship between squared frequencies and compression loads is totally linear, whereas, in others, there are significant deviations from linearity, particularly in the presence of initial geometrical imperfections. Therefore, in these cases, nonlinear analysis is essential to perform accurate analysis and verification of the results. Different approaches were adopted to obtain reliable results of buckling prediction and natural frequency variation. For instance, Souza et al. [31] presented a new relationship between the natural frequencies and the applied compressive loadings. A semi-analytical VCT formulation to compute the critical load of shells was proposed by Jansen et al. [32]. Arbelo et al. [33,34] implemented a modified-VCT based on the considerations made by Souza. An important book that reports a detailed and complete description about experimental VCT setup and results was written by Abramovich [35]. The literature on VCT analyses of isotropic and classical composite 2D structures is vast, while the one on VAT plates and shells is limited. Readers are referred to [36][37][38] for some interesting studies on the dynamic investigation of VAT structures.
The present paper is not primarily intended to show the benefits of variable angle tow designs. This is a topic that past research efforts have already proved, but rather at illustrating an efficient numerical model for studying the eigenfrequencies and eigenmodes of VAT structures in their quasi-static nonlinear (and eventually postcritical) equilibrium states. As a matter of fact, the main objectives are to provide an efficient methodology based on high-accuracy but efficient layerwise (LW) models to investigate the dynamic characteristic of VAT plates and shells under extreme loadings and a means for verifying experimental VCT results.
In this context, the VAT 2D structures are formulated in the Carrera Unified Formulation (CUF) [39] framework in order to obtain accurate results. CUF gives the possibility to develop refined structural models quickly and easily with no ad-hoc definitions of the finite element (FE) matrices. According to CUF, any theory can be degenerated into a generalized kinematics assuming an arbitrary expansion of the generalized variables. Thus, by adopting this procedure, the nonlinear governing equations and the relative FE arrays of the 2D theories are written in terms of Fundamental Nuclei (FNs). The latter are the main building blocks of the presented theory. The plate and shell formulation was introduced by Carrera [40], and it has recently extended to investigate VAT composite structures [41,42]. Multiscale defects in tow-steered components are studied adopting a stochastic response analysis in [43][44][45]. In this research, both the equivalent single layer (ESL) based on the Taylor-like expansions (TE model) and the LW adopting the Lagrange expansion (LE model) are considered.
The content of this paper is arranged as follows: (i) first, preliminary information about the 2D CUF is presented in Sect. 2, presenting the Green-Lagrange nonlinear geometrical relations and techniques adopted for modelling the variation law of the fibre angle; (ii) then, the virtual VCT methodology used in this research is briefly reported in Sect. 3; (iii) next, in Sect. 4 some numerical examples are discussed; (iv) finally, Sect. 5 summarizes the main conclusions.

Preliminary considerations
In this paper, the VAT composite structures are modelled employing refined 2D CUF models. Consider a 2D model that is described by using a Cartesian system (x, y, z) for plate structures and an orthogonal curvilinear system (α,β,z) for shells, in which the in-plane surface is described by x-y and α-β, respectively, and z along the thickness direction. The three-dimensional (3D) displacement, strain and stress vectors of a generic point within the structural domain of a plate can be defined as: where the superscript k indicates the kth layer of the laminated structure. Although derivation is carried out for plates in the following sections, it should be underlined that similar relations hold for shells. The presented methodology adopts a total Lagrangian formulation employing the Green-Lagrange strain tensor components and the second Piola-Kirchhoff stresses σ . Considering a geometrical nonlinear behaviour, the strain-displacement relation is formulated as follows: where b l and b nl stand for the 6 × 3 linear and nonlinear differential operators. These are reported in [46,47] for both plates and shells. Accordingly, the stress vector is computed from the following constitutive equation: where C k is the material linear elastic matrix for orthotropic materials and is defined in [48,49]. By considering a VAT structure, the fibre has a general orientation function of the space coordinates, i.e. θ(x, y) or θ(α, β). Consequently, we obtain where in which T is the rotation matrix [50] and superscript T denotes transposition.
To characterize a VAT structure, the fibre orientation angle is allowed to change continuously along a curvilinear path in each ply. In this way, the laminate has a different stiffness value at each position. In this article, VAT plates/shells with a linear fibre angle variation over the lamina are considered, see Fig. 1a. The fibre path can be described by adopting the well-known expression provided by Gürdal [51]; it reads: in which the fibre path has a rotation of an angle with respect to a certain reference direction placed at an arbitrary point O, see Fig. 1b. The fibre orientation angle at this point is T 0 and varies along a direction x oriented by angle from the original coordinate axis x. The fibre orientation assumes the value T 1 at a characteristic distance d from the reference point. By considering this rotation angle, the fibre orientation path θ(x, y) can be formulated as θ(x ), where x = x cos + y sin . The parameter d is commonly equal to a/2 or b/2 when = 0 • or = 90 • , whereas a and b are the width and length of the 2D structure, respectively. It should be underlined that, as a finite element approximation is adopted in this paper, fibre angle is allowed to vary at each Gauss point in the plate/shell element. In this manner, the FE arrays can be obtained by integrating the material elastic coefficients coherently with the shape function approximation order. This is not always the case of FE models of VAT structures obtained with commercial software tools, where the lamination angle is considered generally constant in the whole element domain. The presented method guarantees a more accurate and efficient analysis of VAT structures. More details can be found in [42].

Displacement field, nonlinear governing equations and FE approximation
By adopting CUF, the 3D displacement field can be defined as an arbitrary through-the-thickness expansion of the in-plane variables: where F τ is a set of thickness expansion functions, u τ represents the generalized displacement vector depending on the in-plane coordinates x and y, M denotes the order of expansion in the thickness direction, t stands for time, and the repeated index τ indicates summation. Depending on the choice of F τ , LW and ESL high-order models can be developed straightforwardly in CUF. Note that in this paper, particular emphasis is given to LW theories based on LE expansion. Offering an independent kinematic description for each layer of the laminated structure, LW models are demonstrated to provide more accurate results and will be preferred over ESL.
Interested readers are referred to [52] for an exhaustive description of LW and ESL. Note that the acronym LEN (LW) and TEN (ESL), considered in this paper, indicate, respectively, LE and TE of order N assumed in the z direction. The FEM is used to approximate the in-plane generalized displacement vector through the shape function in which N i are the ith shape function, q τ i denotes the unknown nodal variables, N n represents the number of nodes per element, and i stands for summation. In this manuscript, the classical 2D nine-node quadratic (Q9) FEs are used for the shape functions in the x−y plane. The nonlinear governing equations of the static problem are formulated adopting the principle of virtual work. Namely, where the virtual variation of the internal strain energy is expressed as: . . . ; M) and N i shape functions (i; j = 1; . . . ; N n ) to obtain the final stiffness matrix. The reader is referred to [39] for the CUF assembly technique. The virtual variation of the external work can be formulated as where p s j stands for the nodal loading vector. Readers are referred to [53] for its complete derivation. After some mathematical operations, the nonlinear equilibrium equations can be expressed in a compact form as Following expansion of the FNs and assembling, one can obtain the final nonlinear algebraic system that can be solved with an incremental linearized scheme based, for example, on the path-following Newton-Raphson method (or tangent method) [54,55].

Vibration around nonlinear equilibrium states
Vibration analysis is briefly introduced in this section. For this purpose, the virtual variation of the inertial loadings is derived in the following: where M i jτ s represents the FN of the mass matrix andq τ i indicates the nodal acceleration vector; the dot stands for time derivative. The derivation of FN of the mass matrix is provided in [39] . Because the modal behaviour of a structure is not a property of the geometrical and mechanical characteristics, but it is a property of the state of equilibrium, eigenfrequencies and eigenmodes may suffer abrupt aberrations in deep nonlinear regimes. To investigate this aspect, the nonlinear governing equations are linearized around a nontrivial equilibrium state as follows: σ nl is the so-called geometric stiffness. The derivation of the tangent stiffness matrix is not given here for the sake of brevity, but it can be found in [46,47].
Displacement variations in Eq. (14) are small so that harmonic vibration can be assumed and the system solved as a linear eigenvalue problem. Summarizing, the presented methodology to investigate the vibration around nonlinear equilibrium states can be described in the following steps: • First, the static geometrical nonlinear problem is solved employing the Newton-Raphson method based on the arc-length approach. • Once the nonlinear equilibrium curve is calculated, the tangent stiffness matrix K T is obtained in each states of interest, see Fig. 2. • Then, by considering the incremental linearized equilibrium condition of Eq. (14) and assuming harmonic displacements around nontrivial equilibrium states, the equations of motion are simplified into a linear eigenvalue problem from which it is possible to evaluate natural frequencies and mode shapes: ω represents the natural frequency and δ q τ i is the eigenvector. • For the sake of clarity, it is important to underline how the nonlinear vibrations exhibit low amplitudes; consequently, it is legitimate to use a linearization around the state of equilibrium for the resolution of the problem.
If the linearization is performed around the trivial equilibrium state (u = 0), the tangent stiffness matrix can be written as where λ is a load factor.

Numerical examples and discussions
First, to explain the adoption of the LE theory, a comparison between the layerwise and equivalent single-layer approach is reported, showing the need to adopt layerwise when VAT structures with complex laminations are investigated. Then, the postbuckling behaviours and the capabilities of the presented approach to provide an efficient methodology based on high-accuracy but efficient LW models to investigate the dynamic characteristic of VAT plates and shells under extreme loadings and a means for verifying experimental VCT results are discussed in this section.

Layerwise vs. equivalent single-layer approach
This section aims to compare the LE and the TE theories in the case of VAT composite structures. In detail, LE is able to provide a layerwise formulation since each layer has a proper description, whereas TE models introduce a global expansion that leads to an equivalent single-layer approach.

Stress analysis
The investigated VAT structure is a square plate characterized by two layers of equal thickness with the following lamination: [0 + < 90/45 > / 0 + < 0/45 >]. This plate, illustrated in Fig. 3 and analyzed in [41], has the following geometrical data: a = b = 1 m and the thickness (h) equal to 0.1 m. All four edges are clamped. The global coordinate system (x, y, z) is located at the centre of the plate. The material is graphite epoxy with the following properties, E 1 = 137.9 GPa, E 2 = E 3 = 8.96 GPa, G 12 = G 13 = 7.1 GPa G 23 = 6.21 GPa, ν 12 = ν 13 = 0.3 and ν 23 = 0.49. The present structure is loaded with a normal pressure P z = −10 kPa applied on the top surface. The convergent model is obtained modelling the structure with at least 10 × 10Q9 for the in-plane mesh approximation and two LD3 in each layer in the z-direction.
The 3D stress distribution, calculated at the point Q = (−a/4, −b/4), is provided in Fig. 4. Both the LW and ESL approaches are adopted.
Clearly, the LW model kinematics should be exploited to evaluate the stress values accurately. The solutions prove that the ESL models are sufficient to only calculate the normal stress, while they are inadequate to predict the transverse shear stress component accurately. Since the next works will concern accurate stress analyses in the nonlinear regime and failure analyses, LE models are already adopted from this work. For clarity, this analysis is presented only for validation and assessment of the presented methodology and to support the choice of the LE model in the next studies.

Nonlinear vibration-buckling investigation
In this section, representative benchmark problems are examined, showing the potentialities of the proposed methodology to carry out virtual VCT in order to compute the critical load and characterize the natural   The results in the following analyses is presented in the nondimensional form using the equation

Square plate with simply supported edge conditions
A 2-layer [0+ < 0/15 > /90+ < 0/15 >] square VAT composite plate is considered as the first example. This structure has the following geometrical data: a = b = 1 m and h = 2 mm. The material properties of the lamina are tabulated in Table 1.
As regards the boundary conditions, simply supported conditions are applied to all the edges, i.e. x = 0, a has v = w = 0, while y = 0, b satisfies u = w = 0. The square plate is subjected to an uniform axial compression in x-direction as illustrated in Fig. 5.
First of all, a convergence analysis is needed for carrying out an accurate investigation. For this square plate model, the convergence is achieved by employing at least 18 × 18Q9 for the in-plane mesh approximation and only one LE1 in each layer in the z-direction. The equilibrium path for the 2-layer [0+ < 0/15 > /90+ < 0/15 >] square VAT composite plate is displayed in Fig. 6, which plots the normalized values of the displacement at the centre of the plate versus the applied compressive load.
For completeness, solutions obtained using the commercial code ABAQUS (ABQ) are also provided for this first case for comparison and validation reasons. In detail, both the ABQ 2D shell model (18 × 18 S8R) and the ABQ 3D solid model (54 × 54 × 2 C3D20R) are used. The results demonstrate that the 2D CUF nonlinear (NL) model agrees well with the ABQ solutions. In addition, also the linear 2D CUF load-deflection curve and the linear buckling load value, representing by the horizontal line, are depicted in Fig. 6. In particular, the linear critical buckling load (P cr ) value is equal to 546.53 N/m. For the sake of clarity, it is important to underline that the considered structure does not show any bifurcation because of the mechanical coupling due to material nonsymmetry. Nevertheless, it is quite common in the literature to identify the linearized buckling loads with the first eigenvalues of the stiffness matrix. Figure 7 illustrates the deformed configuration with the displacements contours based on the 2D CUF NL model and ABQ 3D solid model at the fixed load N x = 2130.71 N/m for the square VAT plate under compressive loads.
Furthermore, the values of the first ten free vibration modes obtained using the CUF and ABAQUS are tabulated in Table 2, whereas the most relevant mode shapes are depicted in Figs. 8, 9, 10, 11, 12, 13 and 14.      Furthermore, the comparison between the variation of the natural frequencies obtained via the trivial linearized (L) and via full nonlinear (NL) approach is reported in Fig. 17.
This figure illustrates the trend of the linearized and nonlinear frequencies associated with the first eight modes. In the case of the linear analysis, it can be pointed out that the first natural frequency, mode (1, 1), disappears when the buckling phenomenon occurs within the structure. Thus, the first mode has a null frequency in correspondence with the buckling load as the tangent stiffness is singular. On the other hand, in the case of nonlinear analysis, different behaviour is observed. In fact, no vibration mode reaches zero frequency. In particular, the mode (1, 1) reaches a minimum and then it increases in the postbuckling regime. This change in slope shows an instability behaviour, but the buckling phenomena are irrelevant in this case. The results  suggest that a nonlinear approach should be adopted to carry out accurate virtual VCT analyses with the aim of predicting the buckling and characterizing the natural frequencies variation. For completeness, the first eight vibration mode shapes of the square VAT composite plate structure at N x = 3.47 N/m are depicted in Fig. 18. Readers are invited to compare these modes with those shown in Figs. 8, 9, 10, 11, 12, 13 and 14 (trivial state, i.e. N x = 0); mode aberration is evident.

Rectangular plate with simply supported edge conditions
A 4-layer [0± < 0/15 >] 2 rectangular VAT laminate plate is investigated as a second assessment. The plate has the following geometrical characteristics: length (a) equal to 0.2 m, the width (b) is 0.05 m and a total thickness (h) of 0.002 m (ply thickness is 0.0005 m). This VAT structure is subjected to uniform axial compressive load in the x-direction, N x as depicted in Fig. 19. All the plate edges are simply supported. In detail, the opposite edges along width x = 0, a satisfy v = w = 0, whereas the other opposite edges along the length y = 0, b have w = 0 at z = 0. In addition, a further constraint condition, u = v = 0 at the centre point of the model is employed to avoid the rigid body motion of the structure. The lamina properties are given in Table 3.
As in the previous case, a convergence analysis on the discretization of the model was performed. As a result, the convergent model for this rectangular VAT plate structures is reached using at least 20 × 5Q9 for Fig. 19 Geometry of the rectangular VAT plate under uniform axial compressive loads Table 3 Material properties of the rectangular VAT plate under uniform axial compression  Figure 20 depicts the equilibrium curves of the rectangular VAT plate computed by the 2D CUF linear model, 2D CUF nonlinear model, and, for the sake of comparison, using the ABQ 3D solid nonlinear model. Furthermore, in this figure (Fig. 20), some linear buckling load values, representing by the horizontal lines, calculated with the CUF method are also reported. In particular, the first linear critical buckling load value (P cr ) is equal to 2.31 × 10 5 N/m. Although bifurcations are not there in the considered problem, buckling shall be identified as the point in which the slope change in the equilibrium curve is not acceptable any more for the given problem. In this case, the linearized buckling has a mere mathematical and not physical meaning. As displayed in Fig. 20, the equilibrium curve obtained by using the proposed approach has a very good correlation compared to the solution computed with the commercial code ABAQUS. The small discrepancies present in the equilibrium curves are mainly due to the different applications of the variable angle of the laminations. In fact, in ABAQUS, the application of VAT is applied to each element, while in CUF, it is performed for each gauss point, giving the possibility of having greater accuracy. For the sake of completeness, the convergent ABQ model used in this analysis is discretized with a fine mesh (60 × 15 × 4) by adopting C3D20R elements  to overcome the mesh instability problem due to the hourglassing. A particular deformed configuration with the displacements contours based on the 2D CUF NL model and ABQ 3D solid model at the fixed load N x = 2.1 × 10 6 N/m for the rectangular VAT plate is illustrated in Fig. 21. Figure 22 depicts the natural frequencies variation versus compressive loading via trivial linearized solution. Instead, the natural frequencies variation versus compressive loading via full nonlinear solution is reported in Fig. 23.
Furthermore, the comparison between the variation of the nondimensional natural frequencies for progressively increasing loads via trivial linearized and full nonlinear approach is depicted in Fig. 24.
The results displayed in this figure show that even for small load values an important difference is presented between the two trends of the frequencies. The results suggest that a nonlinear approach is needed to perform accurate vibration analysis and, in particular, to have a reliable buckling prediction. In fact, it can be observed from Fig. 24 that using the linear approach the natural frequency of the mode (1, 1) tends to zero at the linear buckling load. On the other hand, the nonlinear solutions exhibit a different behaviour. In detail, the mode (2, 1) reaches a minimum at about 3.7 × 10 5 N/m, after which the frequency increases in the postbuckling regime.  This particular behaviour, already highlighted in the literature in numerous nonlinear studies on different panels [56][57][58][59], shows that the natural frequencies increase after the structure had buckled due to the nonlinear effect of postbuckling distortions. This definite change in the slope of the frequency represents a criterion for the nonlinear buckling prediction. Moreover, it is clear that crossing and veering phenomena are completely unforeseen by the trivial linearized approach. These results obtained using the proposed nonlinear virtual VCT approach prove the potentialities of this methodology to compute the buckling load of the structure and to evaluate the natural frequency variation in a large nonlinear regime with high reliability. For the sake of completeness, the first six vibration mode shapes of the rectangular VAT plate structure at the trivial solution N x = 0 N/m and in far nonlinear regime N x = 1.83 × 10 6 N/m are provided in Figs. 25 and 26, respectively. Readers are invited to compare these modes reported in these figures; various crossing phenomena occur and mode aberration is evident.

Hinged VAT composite shell undergoing snap-through
As a final example, a hinged VAT composite shell under uniform axial compression and a transverse load is considered. This nonlinear benchmark problem is very popular due to the snapping behaviour. This structure consists of three layers with the following stacking sequence: [90+ < 0/45 > /0+ < 0/45 > /90+ < 0/45 >]. This shell model has: L = 508 mm, R α = 2540 mm, θ = 0.1 rad and the thickness (h) equal to 12.7 mm. The material data of the present shell structure are tabulated in Table 4. As displayed in Fig. 27, all nodal displacements are restrained along the hinged edges.  The present shell is modelled adopting 10 × 10Q9 for the in-plane mesh approximation and one LE2 in each layer in the thickness direction.
The transverse deflection, evaluated in the middle of the shell, versus the compressive load P β , considering the ratio P β 0 /P 0 = 2.2, is presented in Fig. 28. Moreover, also the linear buckling load value (P cr = 2.07 N/mm) is shown in this figure.
In addition, some representative points with the relative angle at which the local tangent stiffness is calculated are shown in this figure. At each step of the analysis, free vibrations are computed, as presented in Sect. 3, by employing the tangent stiffness matrix obtained in the previous nonlinear static analysis.   Figure 29 illustrates the natural frequencies variation versus compressive loading via trivial linearized solution, while the variation of the natural frequencies obtained via the full nonlinear approach is provided in Fig. 30.
In addition, for the sake of clarity, the comparison between the variation of the first two vibration modes using the two approaches is shown in Fig. 31. In this figure, the red dots and the relative numbers refer to those displayed in Fig. 28.
The results prove that for this type of structure to evaluate the trend of the natural frequencies accurately is needed to perform a nonlinear analysis. In particular, the results provide the possibility to predict both the linear and nonlinear critical buckling load. In detail, the two peaks of the blue line, Fig. 31, indicate Characteristics first ten free vibration mode shapes for the hinged VAT composite shell under compressive and transverse loads the linear (1.25 N/mm) and nonlinear (5.30 N/mm) buckling load, respectively. Finally, for completeness, the first ten vibration mode shapes of the hinged VAT composite shell are reported in Fig. 32. Furthermore, a Modal Assurance Criterion (MAC) graphical representation is illustrated in Fig. 33. MAC is a scalar indicating correspondence between two sets of mode shapes, see [60,61]. In particular, this figure compares the first 10 modes at the state 3 and 5, see Fig. 31. As reported by this graphical representation, the natural modes at state 3 are identical to those relating to state 5, with the exception of modes 7 and 8 which show a crossing phenomena.

Concluding remarks
A novel numerical methodology to study the eigenfrequencies and eigenmodes of variable angle tow (VAT) structures in their quasi-static nonlinear equilibrium states has been proposed in the present manuscript. This method allows to evaluate the critical buckling load of different structures, to characterize the natural frequencies variation for progressively increasing loadings, and to provide a means for verifying experimental Vibration Correlation Technique (VCT) results. To do this, the two-dimensional (2D) Carrera Unified Formulation (CUF) and the layerwise (LW) approach were adopted to investigate different VAT composite plate and shell structures. The results demonstrated the potential of the presented approach to predict the buckling load and to investigate the natural frequencies variation adopting both an approach based on trivial linearized solution and full nonlinear solution. In detail, the need to adopt a nonlinear approach to perform accurate assess- ments was underlined, particularly for the hinged VAT composite shell structure undergoing snap-through. It was shown that eigenfrequencies and eigenmodes can suffer abrupt aberrations in deep nonlinear regimes. Moreover, mode aberration is evident compared to the modes calculated using the full nonlinear approach with those obtained in the trivial state. In addition, for problems in which bifurcations do not exist or whenever the pre-buckling state is unstable, full geometrical nonlinear analysis shall be preferred to trivial linearized solutions. Finally, this method provides reasonable confidence for future applications in this topic. In particular, future works will concern virtual VCT analyses of unstiffened and stiffened VAT cylindrical shell structures.