# Hierarchical Free Vibration Analysis of Variable-Angle Tow Shells Using Unified Formulation

Domenico Andrea Iannotta, Gaetano Giunta, Levent Kirkayak, Marco Montemurro

Montemurro

► To cite this version:

Domenico Andrea Iannotta, Gaetano Giunta, Levent Kirkayak, Marco Montemurro. Hierarchical Free Vibration Analysis of Variable-Angle Tow Shells Using Unified Formulation. Journal of Composites Science, 2025, 9 (2), pp.55. 10.3390/jcs9020055. hal-04998234

HAL Id: hal-04998234

https://hal.inrae.fr/hal-04998234v1

Submitted on 20 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Article

# Hierarchical Free Vibration Analysis of Variable-Angle Tow Shells Using Unified Formulation

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

**Abstract:** This paper investigates the dynamic behavior of shell structures presenting variable-angle tow laminations. The choice of placing fibers along curvilinear patterns allows for a broader structural design space, which is advantageous in several engineering contexts, provided that more complex numerical analyses are managed. In this regard, Carrera's unified formulation has been widely used for studying variable-angle tow plates and shells. This article aims to expand this formulation through the derivation of the complete formulation for a generic shell reference surface. The principle of virtual displacements is used as a variational statement for obtaining, in a weak sense, the stiffness and mass matrices within the finite element solution method. The free vibration problem of singly and doubly curved variable-angle tow shells is then addressed. The proposed approach is compared to Abaqus three-dimensional reference solutions and classical theories to investigate the effectiveness of the developed models in predicting the vibrational frequencies and modes. The results demonstrate a good agreement between the proposed approach and reference solutions.

**Keywords:** modal analyses; finite element method; variable-angle tow composites; Carrera's unified formulation; shells of revolution

## 1. Introduction

The study of shell structures has significantly expanded in recent decades. Composite shells, in particular, have been utilized across various engineering contexts due to their excellent mechanical properties and low weight. The versatility and efficiency of shell structures make them an ideal choice for a wide range of applications, including aerospace and civil engineering and the automotive industry, as well as several other industrial fields. On the other hand, shell theory is more complicated than classical plate theories because the middle surface curvature generates more complex equations of motion. Thanks to the development of new manufacturing technologies, variable-angle tow (VAT) composites can be arranged in shell configurations in order to further improve the mechanical properties of these structures. Among the various processes used for manufacturing VAT composites, automated fiber placement (AFP) and modern 3D printing (3DP) techniques play a pivotal role. AFP allows one to place prepreg composite materials onto a mold or mandrel. Usually, this technique uses a robotic arm that moves the placement head through a predefined path. When the tow lays down along a straight path, conventional composites can be

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

check for
updates

Academic Editor: Francesco
Tornabene

Received: 19 December 2024

Revised: 17 January 2025

Accepted: 21 January 2025

Published: 24 January 2025

**Citation:** Iannotta, D.A.; Giunta, G.; Kirkayak, L.; Montemurro, M. Hierarchical Free Vibration Analysis of Variable-Angle Tow Shells Using Unified Formulation. *J. Compos. Sci.* **2025**, *9*, 55. https://doi.org/10.3390/jcs9020055

**Copyright:** © 2025 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/).

<footnote n="1">1 Luxembourg Institute of Science and Technology, 5 Avenue des Hauts-Fourneaux, L-4362 Esch-sur-Alzette, Luxembourg; gaetano.giunta@list.lu (G.G.); levent.kirkayak@list.lu (L.K.)</footnote>

<footnote n="2">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">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: domenico-andrea.iannotta@list.lu</footnote>
obtained. During the process, a heating system allows one to obtain a tacky resin, which is compressed on the surface to avoid gaps. The choice of the material is crucial and can heavily influence the final overall quality. Three kinds of materials are typically used for the AFP process: thermosets, thermoplastics and dry fibers. The term "prepreg" refers to thermosets and thermoplastics, since the fibers are already pre-impregnated with resin while dry fibers are not. More details about AFP techniques can be found in the work by Brasington et al. [1]. The manufacturing technique strongly influences the actual placement of fibers and the presence of gaps/overlaps. For example, the standard AFP head operates by rotating to ensure that the compaction roller remains perpendicular to the tangent of the reference tow path. During this process, the tow material is deposited by relying on its in-plane bending deformation. However, this can lead to fiber buckling on the inner side and straightening on the outer side of the curved tow path. This kind of problematic is further discussed by Kim, Potter and Weaver [2], who proposed an alternative method called continuous tow shearing (CTS). Three-dimensional printing techniques appear to be very promising in the context of composite manufacturing because of the possibility to represent complex geometries, and offer cost-effective solutions for small-scale production. For example, fusion deposition modeling (FDM) is a process where the material is extruded from a nozzle, enabling the formation of layers. Usually, these kinds of processes use thermoplastic polymers and the fiber reinforcements are constrained by the size of the nozzle. The printing process uses materials that are commonly used for conventional composites, including carbon fibers, glass fibers, Kevlar and natural fibers. More details about the kinds of materials that can be considered within this type of process can be found in the work by Zhuo et al. [3]. Various engineering applications can benefit from the use of VAT composites. For example, Brooks et al. [4] used these materials for aircraft wings to optimize the structural weight and fuel consumption. VATs can also be applied in the aerospace context for air launch vehicle applications as discussed by Grenoble et al. [5]. By placing fibers along curvilinear patterns, it is possible to enhance the stiffness of classical composites. In a more generic view, this leads to a wider structural design space and a broader tailoring choice for various kinds of problems. For example, it can be useful to orient the fibers along specific patterns to improve the search for an optimal solution within an optimization context. This can also be convenient when the analyzed structure presents geometrical discontinuities such as cut-outs. One of the main disadvantages of VATs is the increased complexity of analysis. This is due to the greater number of unknown variables that need to be considered, which can lead to fiber patterns during the optimization process that cannot be actually manufactured. The combination of complexities related to shell geometry and VAT composites highlights the necessity of accurate yet computationally efficient numerical tools for the mechanical investigation of VAT shells.

Various methods have been employed for such analyses. Some of these approaches are briefly discussed in the following text to establish the foundational basis upon which the proposed work is developed. The theoretical concept of shell structures is quite old and has been developed under various theories over the last century. Love [6] provided a general overview of plate vibration and flexure theories, followed by their extension to the case of shells. Additionally, Donnell [7] proposed an exact solution to study the stability of cylindrical walls, demonstrating that experimentally obtained buckling loads are smaller than theoretical ones. A fundamental contribution in shell theory was given by Flügge [8], who studied the distribution of stresses in shells considering several geometries, loads and boundary conditions. Among the various techniques that were used to study shells structures, the finite element method (FEM) has been instrumental in advancing this field. A general overview of the application of the FEM to shells was given by Lindberg, Olson and Cowper [9]. As stated in their work, two approaches can be followed to predict
shells' mechanical behavior through the FEM. The former is based on the idea of using flat elements to discretize a curved surface. This method was followed by Clough and Johnson [10], who performed a static analysis of shell structures through flat triangular elements. The second possible approach consists of developing curved shell elements that allow one to better approximate the geometry of the analyzed structure. As an example of this latter approach, Bogner, Fox and Schmit [11] proposed a cylindrical element to perform a static analysis of a pinched cylinder. The multiple attempts to develop a quadrilateral shell element showed that, in most cases, the elements resulted in excessive stiffness. This happens because of the combination of shear and membrane locking that usually characterizes these elements [12]. Bathe and Dvorkin [13] developed reliable shell elements with four and eight nodes through the use of a mixed interpolation of tensorial components. Subsequently, the same approach was extended by Bucalem and Bathe [14] to elements with nine and sixteen nodes. Among the several theories developed in the context of structural mechanics, Carrera's unified formulation (CUF) has shown good results while maintaining a reasonable computational cost. According to CUF, the preliminary choice of expansion functions along the thickness allows one to obtain a wide family of hierarchic theories within a unique formulation (see Carrera [15]). CUF was initially applied to shell structures through a Navier-type closed-form solution. This approach was used by Carrera et al. [16] to perform a free vibration analysis of cylindrical and spherical shells. Giunta et al. [17] used the same method in a static analysis context. Subsequently, CUF was extended to thermo-mechanical analyses as shown by Giunta et al. [18]. In this work, isotropic and composite beams subjected to thermal loads were studied, utilizing a meshless collocation method with Wendland's radial basis functions. Brischetto, Polit and Carrera [19] proposed an eight-node element with nine degrees of freedom per node to contrast the shear and membrane locking phenomena. Cinefra and Carrera [20] extended the use of CUF to a static analysis of cylindrical structures in the FEM context. Cinefra, Chinosi and Della Croce [21] investigated the performance of the proposed nine-node shell element through a static analysis of classical problems like the pinched shell and the Scordelis–Lo roof. Cinefra, Kumar and Carrera [22] introduced Reissner's mixed variational theorem within CUF shell elements to perform static analyses. Kumar et al. [23] performed a free vibration analysis of doubly curved shells subjected to delamination. Shells of revolution were studied through CUF by Carrera and Zozulya [24–26]. An important contribution in shell theory was given by Tornabene, Liverani and Caligiana [27], who applied the generalized differential quadrature (GDQ) method with the first-order shear deformation theory for a free vibration analysis of anisotropic shells of revolution. Subsequently, this approach was extended to CUF by Tornabene et al. [28] for a static analysis of doubly curved shells.

The potential for placing fibers along curvilinear paths within a shell structural domain is promising for enhancing the mechanical properties of a composite. As shown by Wu et al. [29], new manufacturing technologies could efficiently integrate the VAT concept within a shell geometry. To the best of the authors' knowledge, Hyer and Charette [30] were among the pioneers in introducing the VAT concept. In their work, they utilized the finite element method (FEM) to model the variation in the path of continuous fibers, assuming a constant fiber angle within each finite element. In some more recent developments, Tornabene, Fantuzzi and Bacciocchi [31] used CUF for the static analysis of curvilinear fibers shells, and the numerical problems were solved through the GDQ method. Sánchez-Majano et al. [32] studied the static behavior of VAT cylindrical shells through CUF. Sciascia et al. [33] studied VAT shells considering buckling, free vibrations and prestressed vibrations. A seven-parameter shell formulation was used within the principle of virtual work to obtain an equivalent single-layer model for the VAT cylindrical shells stability study by He et al. [34]. Numerous studies have demonstrated that utilizing curvilinear
fibers can increase the fundamental frequency. For example, Blom et al. [35] maximized the first natural frequency of VSC conical shells, considering manufacturing constraints. Carvalho, Sohouli et al. [36] used a genetic algorithm and shell elements based on first-order shear deformation theory for the maximization of the fundamental frequency. Among the various optimization methods that could be applied to VAT structures, the multi-scale two-level (MS2L) approach stands out for its ability to efficiently handle complex designs by simultaneously considering global and local optimization criteria. This method divides the optimization process into two distinct phases. Initially, the composite is represented as a homogeneous anisotropic plate to determine the optimal distribution of polar parameters, which are the mechanical design variables. In the subsequent phase, the focus shifts to determining the optimal stacking sequence based on the mechanical properties distribution derived from the first phase. More details about the MS2L approach can be found in the studies by Catapano et al. [37], Montemurro and Catapano [38,39] and Fiordilino et al. [40], where both stiffness and buckling optimization problems were solved.

Despite the large amount of articles that have been published on shells, the literature regarding the application of CUF to the free vibration analysis of VAT shells of revolution offers some room for contributions. Toward this end, this article presents a new development of CUF for the dynamic study of curvilinear fiber shells with single and double curvature. Specifically, this work extends the proposed approach to surfaces with non-constant coefficients of the first quadratic form, enabling the consideration of more complex shell geometries. The structure of this paper is as follows. Section 2 introduces the CUF approach and demonstrates its application to VAT shells using a pure displacement formulation. Section 3 considers two cylindrical shells and a spherical segment in order to perform an eigenvalue free vibration analysis considering various boundary conditions. The results are compared with reference solutions obtained in Abaqus using three-dimensional (3D) elements. Finally, Section 4 shows the conclusive remarks.

## 2. Carrera's Unified Formulation

Shells are 3D objects defined by two closely spaced surfaces with the distance between them significantly smaller than the other dimensions of the shell. The set of midpoints between these surfaces forms the middle surface of the shell. The thickness of the shell at any point is the distance between the surfaces measured along a normal to the middle surface at that point. In order to formulate the shell kinematics, a curvilinear reference system is considered: $\alpha$ and $\beta$ are curvilinear axes that lay on the shell mid-surface, while the through-the-thickness axis (or z-axis) is oriented along the thickness of the shell. The thickness dimension is indicated as $h$, while $l_\alpha$ and $l_\beta$ refer to the characteristics lengths along $\alpha$ and $\beta$ axes, respectively. The mathematical description of the shell starts from the definition of the radius vector, which points at the undeformed middle surface of the shell, starting from a global Cartesian reference system defined by X, Y and Z axes (as shown in Figure 1):

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Shell geometry and reference system.</figcaption></figure>
$$ \mathbf{r} = \mathbf{r}(\alpha, \beta). \qquad (1) $$

The infinitesimal variation in the vector **r** can be defined as

$$ d\mathbf{r} = \mathbf{r}_{,\alpha} d\alpha + \mathbf{r}_{,\beta} d\beta, \qquad (2) $$

where subscripts preceded by a comma refer to the derivative along the curvilinear reference system axes. The norms of the derivative vectors can be expressed as

$$ \begin{aligned} |\mathbf{r}_{,\alpha}| &= A, \\ |\mathbf{r}_{,\beta}| &= B, \end{aligned} \qquad (3) $$

where $|\square|$ is a vector norm operator. The following unit vectors are introduced:

$$ \begin{aligned} \mathbf{i}_{\alpha} &= \mathbf{r}_{,\alpha} / A, \\ \mathbf{i}_{\beta} &= \mathbf{r}_{,\beta} / B, \\ \mathbf{i}_{n} &= \mathbf{i}_{\alpha} \times \mathbf{i}_{\beta} / \sin(\chi), \end{aligned} \qquad (4) $$

where $\mathbf{i}_{\alpha}$ and $\mathbf{i}_{\beta}$ are the unit tangent vectors along $\alpha$ and $\beta$ axes, respectively. By its definition, $\mathbf{i}_{n}$ is the unit vector along the thickness and is orthogonal to both $\mathbf{i}_{\alpha}$ and $\mathbf{i}_{\beta}$, being $\times$ the vector cross product operator. The angle between the mid-surface axes $\alpha$ and $\beta$ is indicated as $\chi$. The first quadratic form of the surface corresponds to the right-hand side of the following equation:

$$ d\mathbf{r} \cdot d\mathbf{r} = A^2 d\alpha^2 + 2AB \cos(\chi) d\alpha d\beta + B^2 d\beta^2, \qquad (5) $$

where $\cdot$ is vector scalar product operator. The second quadratic form is related to the problem of finding the curvature of a generic curve on a surface. For the sake of brevity, only the final formula of the normal curvature of the surface is presented:

$$ \frac{1}{R} = - \frac{L d\alpha^2 + 2 M d\alpha d\beta + N d\beta^2}{A^2 d\alpha^2 + 2 A B \cos(\chi) d\alpha d\beta + B^2 d\beta^2}, \qquad (6) $$

where the numerator of the right-hand side of the equation is the second quadratic form of a surface. The following terms are also introduced:

$$ \begin{aligned} L &= \mathbf{r}_{,\alpha\alpha} \mathbf{i}_n, \\ M &= \mathbf{r}_{,\alpha\beta} \mathbf{i}_n, \\ N &= \mathbf{r}_{,\beta\beta} \mathbf{i}_n. \end{aligned} \qquad (7) $$

The curvatures of $\alpha$ and $\beta$ lines can be obtained by considering $\beta$ and $\alpha$ as constants, respectively:

$$ \begin{aligned} \frac{1}{R_{\alpha}} &= -\frac{L}{A^2}, \\ \frac{1}{R_{\beta}} &= -\frac{N}{B^2}. \end{aligned} \qquad (8) $$

The assumption that $\alpha$ and $\beta$ are lines of principal curvature of the surface allows one to simplify the shell equations. It can be shown that the necessary and sufficient conditions for $\alpha$ and $\beta$ to be lines of principal curvature are as follows:

$$ \begin{aligned} \cos(\chi) &= 0, \\ M &= 0. \end{aligned} \qquad (9) $$
The location of a generic point within the shell domain can be written as

$$ \mathbf{R}(\alpha, \beta, z) = \mathbf{r}(\alpha, \beta) + z\mathbf{i}_n(\alpha, \beta). \qquad (10) $$

The square of the modulus of an infinitesimal change in the **R** vector can be written as follows:

$$ d\mathbf{R} \cdot d\mathbf{R} = (d\mathbf{r} + z d\mathbf{i}_n + \mathbf{i}_n dz)(d\mathbf{r} + z d\mathbf{i}_n + \mathbf{i}_n dz) = H_{\alpha}^2 d\alpha^2 + H_{\beta}^2 d\beta^2 + H_{z}^2 dz^2, \quad (11) $$

where

$$ \begin{aligned} H_{\alpha} &= A(1 + z/R_{\alpha}), \\ H_{\beta} &= B(1 + z/R_{\beta}), \\ H_{z} &= 1. \end{aligned} \qquad (12) $$

Subsequently, the infinitesimal surface and volume of the fundamental shell element can be written as

$$ \begin{aligned} d\Omega &= H_{\alpha}H_{\beta}d\alpha d\beta, \\ dV &= H_{\alpha}H_{\beta}H_{z}d\alpha d\beta dz. \end{aligned} \qquad (13) $$

In this work, only structures with constant curvature along $\alpha$ and $\beta$ are considered. Despite this assumption, the governing equations are written in a generic form and can be used when $R_{\alpha}$ and $R_{\beta}$ are varying within the shell mid-surface. More details about the shell geometrical formulation can be found in the studies by Reddy [41], Leissa [42], Kraus [43] and Washizu [44]. The displacement field is expressed as

$$ \mathbf{u} = \begin{cases} u_{\alpha} \\ u_{\beta} \\ u_{z} \end{cases}. \qquad (14) $$

The strain vector in Voigt's notation can be divided into the mid-surface and transverse components:

$$ \epsilon_p = \begin{cases} \epsilon_{\alpha\alpha} \\ \epsilon_{\beta\beta} \\ \epsilon_{\alpha\beta} \end{cases}, \quad \epsilon_n = \begin{cases} \epsilon_{\alpha z} \\ \epsilon_{\beta z} \\ \epsilon_{\alpha z} \end{cases}. \qquad (15) $$

The hypothesis of small displacements allows for using a linear strains-displacements relation:

$$ \begin{aligned} \epsilon_p &= (\mathbf{D}_p + \mathbf{A}_p)\mathbf{u}, \\ \epsilon_n &= (\mathbf{D}_{n\Omega} + \mathbf{A}_n + \mathbf{D}_{nz})\mathbf{u}, \end{aligned} \qquad (16) $$

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

$$ \mathbf{D}_p = \begin{bmatrix} \frac{\partial \alpha}{H_{\alpha}} & 0 & 0 \\ 0 & \frac{\partial \beta}{H_{\beta}} & 0 \\ \frac{\partial \beta}{H_{\beta}} & \frac{\partial \alpha}{H_{\alpha}} & 0 \end{bmatrix}, \quad \mathbf{D}_{n\Omega} = \begin{bmatrix} 0 & 0 & \frac{\partial \alpha}{H_{\alpha}} \\ 0 & 0 & \frac{\partial \beta}{H_{\beta}} \\ 0 & 0 & 0 \end{bmatrix}, \quad \mathbf{D}_{nz} = \begin{bmatrix} \partial z & 0 & 0 \\ 0 & \partial z & 0 \\ 0 & 0 & \partial z \end{bmatrix}. \qquad (17) $$
The matrices $A_p$ and $A_n$ are written as follows:

$$A_p = \begin{bmatrix} 0 & \frac{A_{\alpha\beta}}{BH_{\alpha}} & \frac{A}{H_{\alpha}R_{\alpha}} \\ \frac{B_{\alpha\alpha}}{AH_{\beta}} & 0 & \frac{B}{H_{\beta}R_{\beta}} \\ -\frac{H_{\alpha\beta}}{H_{\alpha}H_{\beta}} & -\frac{H_{\beta\alpha}}{H_{\alpha}H_{\beta}} & 0 \end{bmatrix}, \quad A_{n\Omega} = \begin{bmatrix} -\frac{A}{H_{\alpha}R_{\alpha}} & 0 & 0 \\ 0 & -\frac{B}{H_{\beta}R_{\beta}} & 0 \\ 0 & 0 & 0 \end{bmatrix}. \quad (18)$$

The mid-surface and transverse stress components can be written as

$$\sigma_p = \begin{Bmatrix} \sigma_{\alpha\alpha} \\ \sigma_{\beta\beta} \\ \sigma_{\alpha\beta} \end{Bmatrix}, \quad \sigma_n = \begin{Bmatrix} \sigma_{\alpha z} \\ \sigma_{\beta z} \\ \sigma_{\alpha z} \end{Bmatrix}. \quad (19)$$

Hooke's law reads as follows:

$$\begin{align} \sigma_p &= \tilde{C}_{pp}\epsilon_p + \tilde{C}_{pn}\epsilon_n, \tag{20} \\ \sigma_n &= \tilde{C}_{np}\epsilon_p + \tilde{C}_{nn}\epsilon_n, \nonumber \end{align}$$

where the terms $\tilde{C}_{pp}$, $\tilde{C}_{pn}$, $\tilde{C}_{np}$ and $\tilde{C}_{nn}$ are the components of the material stiffness matrix.

## 2.1. Variable Angle-Tow Composite Shells

When considering VAT shells, the components of the material stiffness matrix vary along the mid-surface directions. It is possible to write an equation that allows one to rotate the material stiffness matrix C according to a specific angle $\theta$ around the z-axis as

$$\tilde{C} = TCT^T, \quad (21)$$

where C is the material stiffness matrix in the material reference system while $\tilde{C}$ is the same matrix after a rotation. Matrix T represents the rotation matrix, which depends on the angle $\theta$. For the sake of brevity, the components of C and T are not reported here; see the study by Reddy [41]. A linear variation law can be expressed as

$$\theta(\alpha) = \Phi + T_0 + \frac{T_1 - T_0}{d} |\alpha'|. \quad (22)$$

The angle $\Phi$ describes the original direction along which $\theta$ varies, and $\alpha'$ is a generic spatial variable obtained as $\alpha' = x' \cos(\Phi) + y' \sin(\Phi)$. $x'$ and $y'$ denote the axes of the fiber reference system. $T_0$ is the starting angle of the fiber when $\alpha' = 0$, while $T_1$ is the angle of the fiber when $\alpha' = d$. Figure 2 represents an example of the fibers' path.

<figure><img src="image_4.png" /><figcaption>Figure 2. Example of in-plane fibers' path.</figcaption></figure>
As shown in Figure 2, the fibers' angle is always measured with respect to the x' axis for all of the analysis cases. The variation direction of θ can be x', y' or a combination of the two; for this reason, it is specified case by case. In this work, the following notation (based upon the above introduced parameters) is used in order to describe the in-plane fibers' path: Φ < T₀, T₁ >. Further details about the fiber variation law can be found in the work by Gürdal et al. [45].

## 2.2. Variational Formulation

The principle of virtual displacements (PVD) is considered in order to obtain the problem governing equations resulting in only displacements being considered as primary variables. In the context of free vibration analyses, the following variational statement applies:

$$ \delta L_{int} + \delta L_{ine} = 0. \qquad (23) $$

The virtual internal work can be expressed as

$$ \delta L_{int} = \int_{\Omega} \int_{h} \left( \delta \epsilon_{pG}^{T} \sigma_{pH} + \delta \epsilon_{nG}^{T} \sigma_{nH} \right) dz d\Omega, \qquad (24) $$

where the subscript 'G' refers to the components obtained from geometrical relations in Equation (16), subscript 'H' refers to the components obtained from Hooke's law in Equation (20), Ω is the middle surface of the shell and subscript 'T' refers to the transpose operator. As far as the virtual work of the inertial forces is concerned, the following equation can be written:

$$ \delta L_{ine} = \int_{\Omega} \int_{h} \delta \mathbf{u}^{T} \rho \ddot{\mathbf{u}} dz d\Omega, \qquad (25) $$

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

## 2.3. Kinematic Assumption and Finite Element Approximation

In order to express the primary unknowns, CUF uses an axiomatic approach along the through-the-thickness direction (see the study by Carrera [15]). The generic unknown component $f = f(\alpha, \beta, z)$ can be approximated as

$$ f(\alpha, \beta, z) = F_{\tau}(z)g_{\tau}(\alpha, \beta), \quad \tau = 0, 1, \dots, N, \qquad (26) $$

where f is a displacement component in a PVD context. $F_{\tau}$ is an ansatz function along the thickness and $g_{\tau}$ is an unknown two-dimensional function accounting for the variation in the mid-surface. According to Einstein's notation, a twice-repeated index implies a sum over that index range. Finally, N is the approximation order. Both N and $F_{\tau}$ can be imposed a priori. This feature of CUF allows one to obtain multiple theories in the same formulation.

According to the choice of $F_{\tau}$, it is possible to obtain equivalent single-layer (ESL) or layer-wise (LW) models. In the ESL case, the number of unknowns does not depend on the number of layers of the shell. For this reason, the total stiffness matrix terms are obtained as a weighted average of the layers' stiffness along the thickness. Taylor's polynomials are the most common choice as ansatz function for developing an ESL model:

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

ESL models are able to predict the general response of relatively thin laminates and are characterized by a reduced computational cost. On the other hand, they are not able to accurately predict the behavior of thick shells, especially in the case of a high degree of anisotropy. Furthermore, since these models are based on $C^{\infty}$ approximation functions, ESL approaches cannot correctly describe the zigzag displacements effect (nonetheless, it
could be possible to include this feature by adding Murakami's function, as explained in the study by Carrera [46], to the approximation set of functions).

In order to improve the accuracy of the model, it is possible to consider an LW approach, where the kinematics of each layer are formulated independently. In this case, the number of unknowns depends on the number of layers and it is necessary to impose the continuity of the problem's main unknowns at an interface between consecutive layers. LW models usually use Lagrange or Legendre polynomials to approximate the unknown fields. In an LW approach, the approximation along the thickness direction reads as follows:

$$f^k(\alpha, \beta, z) = F_b(z)g_b^k(\alpha, \beta) + F_r(z)g_r^k(\alpha, \beta) + F_t(z)g_t^k(\alpha, \beta), \quad r = 2, \dots, N, \quad (28)$$

where the superscript 'k' refers to a specific layer of the shell, with k ranging between one and $N_l$, the latter being the total number of layers. Subscripts 't' and 'b' refer to the top and bottom faces of the generic layer, respectively. In the case of Legendre polynomials, the through-the-thickness approximating functions are

$$F_t(z(\zeta_k)) = \frac{P_0 + P_1}{2}, \quad F_b(z(\zeta_k)) = \frac{P_0 - P_1}{2}, \quad F_r(z(\zeta_k)) = P_r - P_{r-2}, \quad r = 2, \dots, N, \quad (29)$$

where $P_i = P_i(\zeta_k)$ is the ith-order Legendre polynomial defined in the domain of the k-th layer and $-1 \le \zeta_k \le 1$. LW models are able to predict the zigzag through-the-thickness behavior of the displacement field. Nevertheless, LW models demand a higher computational cost since they provide an independent approximation for every layer of the shell.

An FEM solution calls for shape functions being introduced into the approximation. In the case of a bi-dimensional model, Equation (26) becomes

$$f(\alpha, \beta, z) = F_{\tau}(z)N_{i}(\alpha, \beta)g_{\tau i}, \quad \tau = 0, 1, \dots, N, \quad i = 1, \dots, N_{n}, \quad (30)$$

where $N_i$ represents the shape functions that are used for the approximation of the unknowns into the mid-surface of the shell, and $N_n$ is equal to the number of nodes used for the domain discretization. Classical Lagrange shape functions are used. They are not presented here for the sake of brevity but they can be found in the study by Bathe [47].

## 2.4. Acronym System

An acronym system is introduced in order to identify all the derived theories as shown in Figure 3.

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

The first letter refers to the approximation level that is used: 'E' corresponds to ESL models while 'L' corresponds to LW models. The second letter refers to the variational statement: 'D' stands for PVD. The last number refers to the expansion order used along
the structure thickness. The first number, when present, refers to the number of virtual layers that have been used for an LW model to represent every single physical layer; if the number at the beginning of the acronym is not present, it is assumed that only one virtual layer has been used for each physical layer.

As an example of the used acronym system, in EDN models, the displacement field is expressed in the following form:

$$ 
\begin{aligned}
u_{\alpha} &= u_{\alpha 0} + u_{\alpha 1}z + u_{\alpha 2}z^2 + \cdots + u_{\alpha N}z^N, \\
u_{\beta} &= u_{\beta 0} + u_{\beta 1}z + u_{\beta 2}z^2 + \cdots + u_{\beta N}z^N, \\
u_{z} &= u_{z 0} + u_{z 1}z + u_{z 2}z^2 + \cdots + u_{z N}z^N.
\end{aligned} 
\quad (31) $$

In a vectorial 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 (32) $$

where $F_{\tau} = z^{\tau}$ and $\mathbf{u}_{\tau} = \mathbf{u}_{\tau}(\alpha, \beta)$. Also, the first-order shear deformation theory is indicated as FSDT and is obtained as a particular case of an ED1 solution. Specifically, FSDT is obtained through the penalization of the $u_{z1}$ term, and the material stiffness matrix, as classically performed, is reduced to account for a plane stress state and resolve thickness locking.

In LW theories, $N$ refers to the approximation order used in every layer:

$$ \mathbf{u}^k = F_0 \mathbf{u}_0^k + F_1 \mathbf{u}_1^k + \cdots + 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 (33) $$

It can be observed that ESL theories can be considered as a particular case of LW ones. While, in the first case, the integration along the thickness is performed in order to represent composite properties through an equivalent single layer, for the second case, the integration is computed layer by layer. This allows one to represent the kinematic of each layer separately for LW models. Unless otherwise stated, LDN solutions are obtained with Lagrange polynomials with equally spaced nodes.

## 2.5. Locking Correction

The locking phenomenon is common for shells' finite elements and it is well known in the literature. Its correction is performed through the mixed interpolation of tensorial components (MITC) method. This method is based on the idea of considering the interpolation of the strain components on specific points of the element domain, which are called "tying points", and it is implemented by considering the strain tensor in the local coordinate system of each element, defined by $\xi$, $\eta$ and $\zeta$ axes. The location of the tying points depends on the strain component being interpolated and on the number of nodes of the element. In this work, nine-node elements using the MITC approach are used. They are named MITC9 shell elements. Three different patterns are used for the tying points in this case, as shown in Figure 4.

<figure><img src="image_6.png" /><figcaption>Figure 4. Tying points for the CUF MITC9 shell element.</figcaption></figure>
Classical Lagrange polynomials are used as interpolation functions and they are defined referring to the points of each interpolation scheme. These function are grouped into the following vectors:

$$
\begin{align}
\mathbf{N}_{m1} &= [N_{A1}, N_{B1}, N_{C1}, N_{D1}, N_{E1}, N_{F1}], \nonumber \\
\mathbf{N}_{m2} &= [N_{A2}, N_{B2}, N_{C2}, N_{D2}, N_{E2}, N_{F2}], \nonumber \\
\mathbf{N}_{m3} &= [N_{P}, N_{Q}, N_{R}, N_{S}]. \tag{34}
\end{align}
$$

The strain components are interpolated in the following way:

$$
\epsilon_p = \begin{bmatrix} \epsilon_{\xi\xi} \\ \epsilon_{\eta\eta} \\ \epsilon_{\xi\eta} \end{bmatrix} = \begin{bmatrix} \mathbf{N}_{m1} & 0 & 0 \\ 0 & \mathbf{N}_{m2} & 0 \\ 0 & 0 & \mathbf{N}_{m3} \end{bmatrix} \begin{bmatrix} \epsilon_{\xi\xi m1} \\ \epsilon_{\eta\eta m2} \\ \epsilon_{\xi\eta m3} \end{bmatrix} = \mathbf{N1} \begin{bmatrix} \epsilon_{\xi\xi m1} \\ \epsilon_{\eta\eta m2} \\ \epsilon_{\xi\eta m3} \end{bmatrix}, \quad (35)
$$

$$
\epsilon_n = \begin{bmatrix} \epsilon_{\xi z} \\ \epsilon_{\eta z} \\ \epsilon_{zz} \end{bmatrix} = \begin{bmatrix} \mathbf{N}_{m1} & 0 & 0 \\ 0 & \mathbf{N}_{m2} & 0 \\ 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} \epsilon_{\xi z m1} \\ \epsilon_{\eta z m2} \\ \epsilon_{zz} \end{bmatrix} = \mathbf{N2} \begin{bmatrix} \epsilon_{\xi z m1} \\ \epsilon_{\eta z m2} \\ \epsilon_{zz} \end{bmatrix}, \quad (36)
$$

where the subscripts *m1*, *m2* and *m3* indicate that the quantities are calculated in the points (A1, B1, C1, D1, E1, F1), (A2, B2, C2, D2, E2, F2) and (P, Q, R, S), respectively.

2.6. FEM Matrices Expression

In the PVD case, the primary unknown field is the displacements field. Considering Equation (30), the displacements field can be written as follows:

$$
\mathbf{u} = F_{\tau} N_{i} \begin{cases} q_{\alpha\tau i} \\ q_{\beta\tau i} \\ q_{\zeta\tau i} \end{cases} = F_{\tau} N_{i} \mathbf{q}_{\tau i}. \qquad (37)
$$

Through the substitution of Equations (25), (24), (16), (20) and (37), into Equation (23),
the weak form of the PVD method can be obtained:

$$
\begin{equation}
\begin{split}
& \int_{\Omega} \delta \boldsymbol{q}_{\tau i}^{T} [(\mathbf{N}_1 \mathbf{m}_p)^T (N_i \mathbf{I}) \tilde{\mathbf{Z}}_{pp}^{\tau s} (\mathbf{N}_{1o} \mathbf{D}_p) (N_j \mathbf{I}) + (\mathbf{N}_1 \mathbf{m}_p)^T (N_i \mathbf{I}) \tilde{\mathbf{Z}}_{pp}^{\tau s} (\mathbf{N}_{1o} \mathbf{A}_p) (N_j \mathbf{I}) + \\
& \qquad + (\mathbf{N}_1 \mathbf{m}_p)^T (N_i \mathbf{I}) \tilde{\mathbf{Z}}_{pn}^{\tau s} (\mathbf{N}_{2o} \mathbf{D}_{n\Omega}) (N_j \mathbf{I}) + (\mathbf{N}_1 \mathbf{m}_p)^T (N_i \mathbf{I}) \tilde{\mathbf{Z}}_{pn}^{\tau s} (\mathbf{N}_{2o} \mathbf{A}_n) (N_j \mathbf{I}) + \\
& \qquad + (\mathbf{N}_1 \mathbf{m}_p)^T (N_i \mathbf{I}) \tilde{\mathbf{Z}}_{pn}^{\tau s,z} (\mathbf{N}_{2o} (N_j \mathbf{I}) + (\mathbf{N}_1 \mathbf{m}_p)^T (N_i \mathbf{I}) \tilde{\mathbf{Z}}_{pp}^{\tau s} (\mathbf{N}_{1o} \mathbf{D}_p) (N_j \mathbf{I}) + \\
& \qquad + (\mathbf{N}_1 \mathbf{m}_p)^T (N_i \mathbf{I}) \tilde{\mathbf{Z}}_{pp}^{\tau s} (\mathbf{N}_{1o} \mathbf{A}_p) (N_j \mathbf{I}) + (\mathbf{N}_1 \mathbf{m}_p)^T (N_i \mathbf{I}) \tilde{\mathbf{Z}}_{pn}^{\tau s} (\mathbf{N}_{2o} \mathbf{D}_{n\Omega}) (N_j \mathbf{I}) + \\
& \qquad + (\mathbf{N}_1 \mathbf{m}_p)^T (N_i \mathbf{I}) \tilde{\mathbf{Z}}_{pn}^{\tau s} (\mathbf{N}_{2o} \mathbf{A}_n) (N_j \mathbf{I}) + (\mathbf{N}_1 \mathbf{m}_p)^T (N_i \mathbf{I}) \tilde{\mathbf{Z}}_{pn}^{\tau s,z} (\mathbf{N}_{2o} (N_j \mathbf{I}) + \\
& \qquad + (\mathbf{N}_2 \mathbf{m}_p)^T (N_i \mathbf{I}) \tilde{\mathbf{Z}}_{np}^{\tau s} (\mathbf{N}_{1o} \mathbf{D}_p) (N_j \mathbf{I}) + (\mathbf{N}_2 \mathbf{m}_p)^T (N_i \mathbf{I}) \tilde{\mathbf{Z}}_{np}^{\tau s} (\mathbf{N}_{1o} \mathbf{A}_p) (N_j \mathbf{I}) + \\
& \qquad + (\mathbf{N}_2 \mathbf{m}_p)^T (N_i \mathbf{I}) \tilde{\mathbf{Z}}_{np}^{\tau s} (\mathbf{N}_{2o} \mathbf{D}_{n\Omega}) (N_j \mathbf{I}) + (\mathbf{N}_2 \mathbf{m}_p)^T (N_i \mathbf{I}) \tilde{\mathbf{Z}}_{np}^{\tau s} (\mathbf{N}_{2o} \mathbf{A}_n) (N_j \mathbf{I}) + \\
& \qquad + (\mathbf{N}_2 \mathbf{m}_p)^T (N_i \mathbf{I}) \tilde{\mathbf{Z}}_{nn}^{\tau s} (\mathbf{N}_{2o} (N_j \mathbf{I}) + (\mathbf{N}_2 \mathbf{m}_p)^T (N_i \mathbf{I}) \tilde{\mathbf{Z}}_{np}^{\tau s} (\mathbf{N}_{1o} \mathbf{D}_p) (N_j \mathbf{I}) + \\
& \qquad + (\mathbf{N}_2 \mathbf{m}_p)^T (N_i \mathbf{I}) \tilde{\mathbf{Z}}_{np}^{\tau s} (\mathbf{N}_{1o} \mathbf{A}_p) (N_j \mathbf{I}) + (\mathbf{N}_2 \mathbf{m}_p)^T (N_i \mathbf{I}) \tilde{\mathbf{Z}}_{nn}^{\tau s} (\mathbf{N}_{2o} \mathbf{D}_{n\Omega}) (N_j \mathbf{I}) + \\
& \qquad + (\mathbf{N}_2 \mathbf{m}_p)^T (N_i \mathbf{I}) \tilde{\mathbf{Z}}_{nn}^{\tau s} (\mathbf{N}_{2o} \mathbf{A}_n) (N_j \mathbf{I}) + (\mathbf{N}_2 \mathbf{m}_p)^T (N_i \mathbf{I}) \tilde{\mathbf{Z}}_{nn}^{\tau s} (\mathbf{N}_{2o} (N_j \mathbf{I}) + \\
& \qquad + \mathbf{N}_2 (N_i \mathbf{I}) \tilde{\mathbf{Z}}_{np}^{\tau,zs} (\mathbf{N}_{1o} \mathbf{D}_p) (N_j \mathbf{I}) + \mathbf{N}_2 (N_i \mathbf{I}) \tilde{\mathbf{Z}}_{np}^{\tau,zs} (\mathbf{N}_{1o} \mathbf{A}_p) (N_j \mathbf{I}) + \\
& \qquad + \mathbf{N}_2 (N_i \mathbf{I}) \tilde{\mathbf{Z}}_{np}^{\tau,zs} (\mathbf{N}_{2o} \mathbf{D}_{n\Omega}) (N_j \mathbf{I}) + \mathbf{N}_2 (N_i \mathbf{I}) \tilde{\mathbf{Z}}_{np}^{\tau,zs} (\mathbf{N}_{2o} \mathbf{A}_n) (N_j \mathbf{I}) + \\
& \qquad + \mathbf{N}_2 (N_i \mathbf{I}) \tilde{\mathbf{Z}}_{nn}^{\tau,zs} (\mathbf{N}_{2o} (N_j \mathbf{I}) + \mathbf{N}_2 (N_i \mathbf{I}) \tilde{\mathbf{Z}}_{nn}^{\tau,zs} (\mathbf{N}_{2o} \mathbf{A}_n) (N_j \mathbf{I}) + \\
& \qquad + \mathbf{N}_2 (N_i \mathbf{I}) \tilde{\mathbf{Z}}_{nn}^{\tau,zs} (\mathbf{N}_{2o} (N_j \mathbf{I})] \boldsymbol{q}_{sj} d\Omega = - \int_{\Omega} \delta \boldsymbol{q}_{\tau i}^{T} (N_i \mathbf{I}) \rho E_{\tau s} (N_j \mathbf{I}) \ddot{\boldsymbol{q}}_{sj} d\Omega,
\end{split}
\tag{38}
\end{equation}
$$
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 (39) $$

$$ (E_{\tau s}, E_{\tau_{z}s}, E_{\tau_{s,z}}, E_{\tau_{z}s_{z}}) = \int_{h} (F_{\tau}F_{s}, F_{\tau_{z}}F_{s}, F_{\tau}F_{s_{z}}, F_{\tau_{z}}F_{s_{z}}) dz . \quad (40) $$

Indices $\tau, i$ and $m$ refer to virtual variations, while $s, j$ and $o$ refer to real quantities. In a compact vectorial form, Equation (38) reads as follows:

$$ \delta q_{\tau i}^{T} \mathbf{K}^{\tau s i j} q_{s j} = -\delta q_{\tau i}^{T} \mathbf{M}^{\tau s i j} \ddot{\mathbf{q}}_{s j}, \quad (41) $$

where $\mathbf{K}^{\tau s i j}$ and $\mathbf{M}^{\tau s i j}$ are the $3 \times 3$ fundamental nuclei (FN) of the stiffness and mass matrices, respectively. The components of the FN are reported in Appendix A. Through the loops on the indices $\tau, s, i$ and $j$, it is possible to build the stiffness and mass matrices of the whole shell. The mid-surface integrals are computed through Gauss quadrature with a $3 \times 3$ grid of points to integrate all the terms of the FN. The algorithm is implemented in Python and the eigenvalue problem described by Equation (38) is solved using the scipy.linalg.eigh function from the SciPy library. This function utilizes LAPACK drivers to solve standard or generalized eigenvalue problems for complex Hermitian or real symmetric matrices.

### 3. Numerical Results and Discussion

Three benchmark problems are considered in this section: a monolayer cylindrical shell, a monolayer cylinder and a bi-layer spherical segment. Parametric studies are performed considering different values of thickness as the other geometrical parameter constant. Material properties are represented in Table 1. The subscript “L” refers to the longitudinal direction of the fibers, while subscript “T” refers to the transverse fibers’ direction.

<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> [GPa]</th><th>G<sub>TT</sub> [GPa]</th><th>ν<sub>LT</sub>, ν<sub>TT</sub></th></tr></thead><tbody><tr><td>1, 2, 3</td><td>50</td><td>10</td><td>5</td><td>5</td><td>0.25</td></tr></tbody></table>

Reference solutions are represented by Abaqus 3D models, where quadratic solid elements with reduced integration and three degrees of freedom per node (C3D20R) are used. In Abaqus 3D solutions, the orientation is assumed to remain constant within each element and corresponds to the fiber angle at the centroid of the respective element. Consequently, a refined mesh is required to achieve accurate results. The material orientation is assigned on an element-by-element basis using Python scripting (version 3.9.10) in Abaqus 2018. A frequency step is defined to compute the first fifteen eigenvalues, with the default Lanczos eigensolver employed. The analysis utilizes three CPUs for parallel processing to enhance computational efficiency. For CUF solutions, nine-node square elements are used (QUAD9). The software Gmsh is used to generate the mesh and represent the modal shapes during post-processing. For each case study, an initial convergence analysis is performed to determine the suitable mesh for both CUF and Abaqus solutions. Percentage errors are evaluated as follows:

$$ \delta_{\text{err}\%} = \frac{|\phi_{i, \text{CUF}} - \phi_{i, \text{Abaqus3D}}|}{|\phi_{i, \text{Abaqus3D}}|} \cdot 100, \quad (42) $$

where $\phi_i$ is the modal frequency corresponding to the $i$-th mode. The effective mass percentage is used to obtain a deeper physical understanding of the modal behavior of the analyzed structures. Indeed, this value represents the amount of the total mass that
participates in a given mode of vibration. This parameter is determined through the Abaqus 3D model and can be calculated as follows:

$$m_{\text{eff},i} = \frac{m_{\text{eff},i}}{m_{\text{tot}}} \cdot 100, \quad i = X, Y, Z. \qquad (43)$$

Here, $m_{\text{eff},i}$ represents the effective mass in one of the directions defined by the axes of the global Cartesian reference system, while $m_{\text{tot}}$ denotes the total mass of the shell. To ensure that the number of considered modes is sufficient for describing the modal behavior of the structure, the sum of the effective masses along each direction should closely approximate the total mass of the system. However, since this is not the primary objective of this work, only the first six natural frequencies are considered for all the studied cases. Further details regarding the computation of effective masses are omitted for the sake of brevity but can be found in the Abaqus Theory Manual [48].

### 3.1. Monolayer Cylindrical Shell

The first case corresponds to a monolayer portion of the cylindrical surface with the following dimensions: $R_α = 1$ m, $l_α = R_απ/4$ m, $l_β = 2l_α$. The reference systems that have been used to describe the problem are represented in Figure 5.

<figure><img src="image_7.png" /><figcaption>Figure 5. Geometry and reference system, case 1.</figcaption></figure>

The radius vector is defined as

$$\mathbf{r} = \begin{cases} R_\alpha \sin(\gamma) \\ Y \\ R_\alpha \cos(\gamma) \end{cases}. \qquad (44)$$

Considering $α = R_αγ$ and $β = Y$, the coefficients of the first quadratic form are defined as follows:

$$A = 1, \quad B = 1, \qquad (45)$$

and, as a consequence,

$$H_α = \left(1 + \frac{z}{R_α}\right), \quad H_β = 1. \qquad (46)$$

It is assumed that the fibers' angle is a function of $y'$. For this problem, axes $x'$ and $y'$ of the local reference system of the fibers' path are coincident with axes $α$ and $β$ of the curvilinear reference system of the shell. Therefore, the characteristic length of Equation (22) is set to $d = l_β$. The angle variational law is expressed as $90 < 0, 90 >$, and it is taken from Viglietti et al. [49], where it is applied on a rectangular plate for vibration analyses.
The fiber orientation within the shell is illustrated in Figure 6, where the axis labeled "1" represents the tangent direction of a fiber at a specific point.

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Fibers' orientation, case 1.</figcaption></figure>

The shell is fully clamped with density $ρ = 1540$ kg/m³. The Abaqus reference solution contains 12 elements along the thickness, while 40 and 80 elements are used along $α$ and $β$, respectively. For the CUF results, a 10 × 20 mesh is considered. Table 2 shows the number of degrees of freedom (DoF) for the considered solutions. FSDT shows the smallest number of DoF since it is obtained through a manipulation of the ED1 model.

<table><caption>**Table 2.** Number of degrees of freedom, case 1.</caption><thead><tr><td></td><td>DOF</td></tr></thead><tbody><tr><td>Abaqus 3D</td><td>503,355</td></tr><tr><td>2LD4</td><td>23,247</td></tr><tr><td>2LD2</td><td>12,915</td></tr><tr><td>ED6</td><td>18,081</td></tr><tr><td>ED4</td><td>12,915</td></tr><tr><td>ED2</td><td>7749</td></tr><tr><td>FSDT</td><td>5166</td></tr></tbody></table>

It is possible to observe that the Abaqus reference solution is characterized by a number of DoF which is at least one order of magnitude bigger than CUF solutions. This shows that the developed models guarantee a reduction in the computational cost in terms of DoF. Table 3 shows the first six modal frequencies for different theories in the case $l_α/h = 100$.

For this case, all the theories show a comparable level of accuracy. This happens because a thin shell is considered and a linear law is sufficient for approximating the displacements along the thickness, as shown by the FSDT theory, which is characterized by a maximum error of 0.29% for the sixth frequency. When the aspect ratio is reduced, the thickness of the shell is increased and the differences among the various theories become more evident. Table 4 shows the first six modal frequencies considering $l_α/h = 10$.
<table><caption>**Table 3.** Natural frequencies (10²·rad/s), $l_\alpha/h = 100$, case 1.</caption><thead><tr><th></th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th></tr></thead><tbody><tr><td>Abaqus 3D</td><td>10.843</td><td>13.862</td><td>16.123</td><td>18.316</td><td>20.670</td><td>21.881</td></tr><tr><td>2LD4</td><td>10.858</td><td>13.911</td><td>16.145</td><td>18.374</td><td>20.701</td><td>21.946</td></tr><tr><td>2LD2</td><td>10.859</td><td>13.912</td><td>16.145</td><td>18.375</td><td>20.702</td><td>21.947</td></tr><tr><td>ED6</td><td>10.858</td><td>13.911</td><td>16.145</td><td>18.374</td><td>20.701</td><td>21.946</td></tr><tr><td>ED4</td><td>10.858</td><td>13.911</td><td>16.145</td><td>18.374</td><td>20.701</td><td>21.946</td></tr><tr><td>ED2</td><td>10.861</td><td>13.918</td><td>16.149</td><td>18.386</td><td>20.708</td><td>21.963</td></tr><tr><td>FSDT</td><td>10.852</td><td>13.893</td><td>16.142</td><td>18.366</td><td>20.700</td><td>21.944</td></tr></tbody></table>

<table><caption>**Table 4.** Natural frequencies (10²·rad/s), $l_\alpha/h = 10$, case 1.</caption><thead><tr><th></th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th></tr></thead><tbody><tr><td>Abaqus 3D</td><td>40.262</td><td>52.831</td><td>63.309</td><td>65.599</td><td>77.539</td><td>81.103</td></tr><tr><td>2LD4</td><td>40.283</td><td>52.854</td><td>63.353</td><td>65.669</td><td>77.597</td><td>81.175</td></tr><tr><td>2LD2</td><td>40.309</td><td>52.900</td><td>63.429</td><td>65.787</td><td>77.723</td><td>81.371</td></tr><tr><td>ED6</td><td>40.283</td><td>52.854</td><td>63.353</td><td>65.669</td><td>77.597</td><td>81.176</td></tr><tr><td>ED4</td><td>40.286</td><td>52.858</td><td>63.362</td><td>65.685</td><td>77.613</td><td>81.202</td></tr><tr><td>ED2</td><td>40.596</td><td>53.404</td><td>64.197</td><td>66.943</td><td>78.937</td><td>83.210</td></tr><tr><td>FSDT</td><td>40.551</td><td>53.374</td><td>64.144</td><td>66.830</td><td>78.881</td><td>83.132</td></tr></tbody></table>

It is possible to observe that, in this case, the FSDT model shows the worst approximation and the percentage error on the sixth mode grows to 2.50%. This error can be reduced to 0.09% with 2LD4 and ED6 models. Moreover, the error of FSDT grows when higher frequencies are considered. This is less evident for higher-order theories, which are able to better predict the reference solution even for higher frequencies. The first and the second modes show two and three semi-waves along the $\alpha$-axis, respectively, and one semi-wave along $\beta$. The third mode exhibits two semi-waves along both $\alpha$ and $\beta$ directions. The fourth mode shows three semi-waves along $\alpha$ and two semi-waves along $\beta$. The fifth mode presents two semi-waves along $\alpha$ and three semi-waves along $\beta$. Finally, the sixth mode shows three semi-waves along $\alpha$ and $\beta$. Table 5 shows the percentage of effective masses for the first six modal frequencies considering $l_\alpha/h = 10$.

<table><caption>**Table 5.** Effective mass percentage $m_{\text{eff\%},i}$, $l_\alpha/h = 10$, case 1.</caption><thead><tr><th></th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th></tr></thead><tbody><tr><td>$m_{\text{eff\%},X}$</td><td>4.7</td><td>1.0</td><td>-*</td><td>7.8</td><td>0.3</td><td>1.7</td></tr><tr><td>$m_{\text{eff\%},Y}$</td><td>-*</td><td>-*</td><td>-*</td><td>-*</td><td>-*</td><td>0.2</td></tr><tr><td>$m_{\text{eff\%},Z}$</td><td>34.8</td><td>3.9</td><td>11.7</td><td>0.3</td><td>0.2</td><td>-*</td></tr></tbody></table>

As evident from the table, the first mode is the most significant among the presented modes, as it involves the highest effective mass percentage along the Z-axis. This behavior can be attributed to the boundary conditions of the analyzed shell: since the structure is clamped and has a limited curvature, the first mode typically corresponds to a predominant displacement along the axis perpendicular to the boundaries, as observed in this case. The displacement of the structure mass along the Z-axis is most prominent in the first three modes. In contrast, the fourth mode exhibits a predominant displacement along the X-axis, while the fifth mode does not demonstrate a principal displacement direction compared to the previous ones. The sixth mode shows a slightly higher effective mass percentage along the X-axis. In this case, the first, second and third modes show one, two and three semi-waves along $\beta$, respectively, and one semi-wave along $\alpha$. The number of semi-waves

<footnote>* Values below $10^{-1}$ are considered as negligible and marked with “-”.</footnote>
increases in both directions for higher modes. Table 6 shows the first six modal frequencies for $l_\alpha/h = 5$. The FSDT maximum error increases to 3.13% for the sixth mode and can be reduced to 0.4% via a 2LD4 model. In this case, the first three modes and the fifth one are similar to the $l_\alpha/h = 10$ shell. The fourth mode shows two semi-waves along $\alpha$ and one semi-wave along $\beta$. The sixth mode displays four semi-waves along $\beta$ and one along $\alpha$.

<table><caption>**Table 6.** Natural frequencies ($10^2 \cdot \text{rad/s}$), $l_\alpha/h = 5$, case 1.</caption><thead><tr><th></th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th></tr></thead><tbody><tr><td>Abaqus 3D</td><td>55.388</td><td>70.438</td><td>86.799</td><td>93.100</td><td>103.030</td><td>108.270</td></tr><tr><td>2LD4</td><td>55.402</td><td>70.459</td><td>86.828</td><td>93.155</td><td>103.049</td><td>108.317</td></tr><tr><td>2LD2</td><td>55.545</td><td>70.674</td><td>87.128</td><td>93.550</td><td>103.195</td><td>108.755</td></tr><tr><td>ED6</td><td>55.403</td><td>70.460</td><td>86.829</td><td>93.157</td><td>103.050</td><td>108.318</td></tr><tr><td>ED4</td><td>55.425</td><td>70.494</td><td>86.881</td><td>93.239</td><td>103.080</td><td>108.401</td></tr><tr><td>ED2</td><td>56.714</td><td>72.344</td><td>89.328</td><td>95.873</td><td>104.149</td><td>111.735</td></tr><tr><td>FSDT</td><td>56.610</td><td>72.277</td><td>89.261</td><td>95.692</td><td>103.999</td><td>111.657</td></tr></tbody></table>

## 3.2. Monolayer Cylinder

The second case corresponds to a cylindrical surface with a single layer. The following dimensions are considered: $R_\alpha = 1$ m, $l_\alpha = 2R_\alpha\pi$, $l_\beta = 2l_\alpha$. The same reference systems, radius vector and coefficients of the first quadratic form of the middle surface as for the previous case are used. Also, the same $90 < 0, 90 >$ variation law of the previous case is considered for the fibers' angle. The fiber orientation within the cylinder is illustrated in Figure 7, where the axis labeled "1" represents the tangent direction of a fiber at a specific point.

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Fibers' orientation, case 2.</figcaption></figure>
For this problem, axes $x'$ and $y'$ of the local reference system of the fibers' path are coincident with axes $\alpha$ and $\beta$ of the curvilinear reference system of the shell. In this case, $l_{\beta}$ is considered as the characteristic length in Equation (22) ($d = l_{\beta}$). As for the previous case, the Abaqus reference solution contains 40 elements along $\alpha$, 80 elements along $\beta$ and 12 elements along the thickness. For CUF results, 8 elements are used in the $\alpha$ direction, while 16 elements are used along $\beta$. The cylinder is clamped at its ends ($\beta = 0$ and $\beta = l_{\beta}$). Since the structure is axisymmetric, frequencies with a multiplicity higher than one are expected. The corresponding modal shapes are the same, apart from a rotation around the main axis of the cylinder. For the sake of brevity, repeated modes are presented only once. The first six modal frequencies are shown in Table 7 in the case of $l_{\alpha}/h = 100$.

<table><caption>**Table 7.** Natural frequencies (10· rad/s), $l_{\alpha}/h = 100$, case 2.</caption><thead><tr><th></th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th></tr></thead><tbody><tr><td>Abaqus 3D</td><td>21.656</td><td>24.166</td><td>32.512</td><td>43.473</td><td>46.874</td><td>48.037</td></tr><tr><td>2LD4</td><td>21.737</td><td>24.172</td><td>32.531</td><td>45.266</td><td>46.876</td><td>47.853</td></tr><tr><td>2LD2</td><td>21.738</td><td>24.172</td><td>32.532</td><td>45.269</td><td>46.876</td><td>47.855</td></tr><tr><td>ED6</td><td>21.737</td><td>24.172</td><td>32.531</td><td>45.266</td><td>46.876</td><td>47.853</td></tr><tr><td>ED4</td><td>21.737</td><td>24.172</td><td>32.531</td><td>45.266</td><td>46.876</td><td>47.853</td></tr><tr><td>ED2</td><td>21.750</td><td>24.173</td><td>32.550</td><td>45.301</td><td>46.876</td><td>47.873</td></tr><tr><td>FSDT</td><td>21.760</td><td>24.155</td><td>32.590</td><td>45.277</td><td>46.874</td><td>47.928</td></tr></tbody></table>

In this case, the best approximation is given by 2LD4 and ED6 models, which both show a maximum error of 4.12% for the fourth modal frequency. Since a thin shell is considered, transverse stresses do not play an important role. Therefore, classical and low-order theories provide good results. The first frequency corresponds to a radial mode with two semi-waves in the circumferential direction and one semi-wave in the longitudinal direction. The second mode exhibits bending with one semi-wave along the cylinder axis. The third mode is a radial mode characterized by two semi-waves in the circumferential direction and two longitudinal semi-waves. The fourth mode involves radial displacements with three circumferential semi-waves and one semi-wave along the cylinder axis. The fifth mode represents a radial mode with no semi-waves in the circumferential direction and one semi-wave in the longitudinal direction. Finally, the sixth mode features two semi-waves in the circumferential direction and three semi-waves in the longitudinal direction. Table 8 shows the first six modal frequencies for $l_{\alpha}/h = 10$.

<table><caption>**Table 8.** Natural frequencies (10· rad/s), $l_{\alpha}/h = 10$, case 2.</caption><thead><tr><th></th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th></tr></thead><tbody><tr><td>Abaqus 3D</td><td>26.336</td><td>47.354</td><td>57.253</td><td>88.047</td><td>94.851</td><td>103.540</td></tr><tr><td>2LD4</td><td>26.338</td><td>47.364</td><td>57.258</td><td>88.065</td><td>94.879</td><td>103.581</td></tr><tr><td>2LD2</td><td>26.346</td><td>47.365</td><td>57.284</td><td>88.067</td><td>94.934</td><td>103.585</td></tr><tr><td>ED6</td><td>26.338</td><td>47.364</td><td>57.258</td><td>88.065</td><td>94.879</td><td>103.581</td></tr><tr><td>ED4</td><td>26.340</td><td>47.364</td><td>57.263</td><td>88.066</td><td>94.889</td><td>103.581</td></tr><tr><td>ED2</td><td>26.397</td><td>47.374</td><td>57.473</td><td>88.075</td><td>95.344</td><td>103.621</td></tr><tr><td>FSDT</td><td>26.405</td><td>47.256</td><td>57.537</td><td>87.897</td><td>95.515</td><td>103.286</td></tr></tbody></table>

The best approximation is given by 2LD4 and ED6 models, which show a maximum error of 0.04% for the sixth frequency. The modal shapes obtained through the ED4 model for $l_{\alpha}/h = 10$ are shown in Figure 8. These shapes match the ones obtained in Abaqus, which are not reported here for the sake of brevity. Mode one, three and five are bending modes. Modes two and four exhibit both translation along the Y-axis and rotation about the same axis. Finally, in the sixth mode, radial and longitudinal displacements can be
observed. The deformation associated with the sixth mode is concentrated near the bottom of the cylinder, close to $y = 0$. This behavior can be explained by considering the fiber orientation: at the bottom, the fibers are aligned axially, whereas at the top, they are oriented circumferentially. As a result, the structure exhibits greater stiffness in the radial direction at $y = l_β$.

<figure><img src="image_10.png" /><figcaption>**Figure 8.** First six eigenmodes, $l_α/h = 10$, case 2.</figcaption></figure>

Table 9 shows the percentage of effective masses for the first six modal frequencies considering $l_α/h = 10$.

<table><caption>**Table 9.** Effective mass percentage $m_{eff\%,i}$, $l_α/h = 10$, case 2.</caption><thead><tr><th></th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th></tr></thead><tbody><tr><td>$m_{eff\%,X}$</td><td>0.6</td><td>-*</td><td>2.7</td><td>-*</td><td>1.8</td><td>-*</td></tr><tr><td>$m_{eff\%,Y}$</td><td>-*</td><td>1.7</td><td>-*</td><td>59.4</td><td>-*</td><td>9.7</td></tr><tr><td>$m_{eff\%,Z}$</td><td>72.1</td><td>-*</td><td>-*</td><td>-*</td><td>9.2</td><td>-*</td></tr></tbody></table>

In this case, the predominant modes are the first and the fourth. The first mode involves significant displacement of the mass along the Z-axis. The fourth mode, as previously discussed, couples longitudinal displacement along the Y-axis with rotation about the same axis. Nevertheless, the table shows that the majority of the mass moves in the Y direction, indicating that the fourth mode is primarily a longitudinal mode. The same considerations can be applied to the sixth mode. Table 10 shows the first six modal fre-

<footnote>* Values below $10^{-1}$ are considered as negligible and marked with “-”.</footnote>
quencies for $l_\alpha/h = 5$. In this case, classical and lower-order theories are not sufficient for matching the modal frequencies of the reference solution because a thick shell is considered. For example, the FSDT shows the inversion of the fifth and sixth vibration modes due to its limitations in capturing the complex couplings inherent to these modes. Specifically, the fifth mode exhibits torsional-radial coupling, while the sixth mode involves a complex interaction between torsional and bending effects. Accurately modeling these phenomena requires the inclusion of higher-order deformation effects, which are accounted for in more advanced theories, enabling more precise and reliable predictions. It is also possible to notice that the difference between 2LD4 and ED6 theories is not evident because a monolayer structure is considered. Indeed, both these theories show the best approximation, which is characterized by a maximum error of 0.06% for the sixth natural frequency. The first five modes of this case are similar to the ones of the $l_\alpha/h = 10$ cylinder. The sixth mode shows radial displacements with zero semi-waves along $\alpha$ and two semi-waves along $\beta$.

<table><caption>**Table 10.** Natural frequencies (10· rad/s), $l_\alpha/h = 5$, case 2.</caption><thead><tr><th></th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th></tr></thead><tbody><tr><td>Abaqus 3D</td><td>29.694</td><td>47.943</td><td>64.012</td><td>88.485</td><td>105.250</td><td>106.100</td></tr><tr><td>2LD4</td><td>29.697</td><td>47.958</td><td>64.017</td><td>88.495</td><td>105.278</td><td>106.165</td></tr><tr><td>2LD2</td><td>29.728</td><td>47.964</td><td>64.109</td><td>88.504</td><td>105.447</td><td>106.200</td></tr><tr><td>ED6</td><td>29.697</td><td>47.958</td><td>64.017</td><td>88.495</td><td>105.278</td><td>106.165</td></tr><tr><td>ED4</td><td>29.705</td><td>47.959</td><td>64.041</td><td>88.497</td><td>105.323</td><td>106.171</td></tr><tr><td>ED2</td><td>29.874</td><td>47.994</td><td>64.540</td><td>88.528</td><td>106.222</td><td>106.339</td></tr><tr><td>FSDT</td><td>29.972</td><td>47.984</td><td>64.819</td><td>88.273</td><td>106.695</td><td>106.618</td></tr></tbody></table>

### 3.3. Bi-Layer Spherical Segment

The last case corresponds to a spherical segment with two layers of equal thickness. This kind of surface is obtained by cutting a sphere with two parallel planes, with the first one passing by the sphere center. The geometry and reference systems of this case are shown in Figure 9. The following dimensions are considered: $R = R_\alpha = R_\beta = 1$ m, $l_\alpha = R\pi/6$, $l_\beta = 2R\pi$.

The radius vector is defined as

$$ \mathbf{r} = \begin{cases} R \sin(\gamma) \\ R \cos(\gamma) \sin(\theta) \\ R \cos(\gamma) \cos(\theta) \end{cases}. \qquad (47) $$

Considering $\alpha = R\gamma$ and $\beta = R\theta$, the coefficients of the first quadratic form are defined as follows:

$$ A = 1, \quad B = \cos(\gamma). \qquad (48) $$

As a consequence,

$$ H_{\alpha} = \left(1 + \frac{z}{R}\right), \quad H_{\beta} = \cos(\gamma)\left(1 + \frac{z}{R}\right). \qquad (49) $$

The angle variation of the fibers is such that continuity along the circumferential direction is ensured. It is assumed that the fibers' angle is a function of $y'$. Axes $x'$ and $y'$ of the local fibers' reference system are parallel to $\alpha$ and $\beta$ and their origin $o'$ is placed in correspondence of $\alpha = l_\alpha/2$ and $\beta = l_\beta/2$. Therefore, the characteristic length in Equation (22) is $d = l_\beta/2$. The stacking sequence can be expressed as $[90 < 0, -45 > / 90 < 0, 45 >]$. The fiber orientation within the sphere is illustrated in Figure 10, where the axis labeled "1" represents the tangent direction of a fiber at a specific point.
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Geometry and reference system, case 3.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Stacking sequence, case 3.</figcaption></figure>

The Abaqus reference solution uses 18 elements along $\alpha$, 216 elements along $\beta$ and 12 elements along thickness. CUF solutions use three elements in the $\alpha$ direction and 36 elements in the $\beta$ direction. The spherical segment is clamped at the lower edge ($\alpha = 0$), while the top edge is free ($\alpha = l_{\alpha}$). As already explained for the previous case, only the first frequency of each couple is presented in the following tables. The first six modal frequencies for $l_{\alpha}/h = 100$ are presented in Table 11.

<table><caption>**Table 11.** Natural frequencies (10·rad/s), $l_{\alpha}/h = 100$, case 3.</caption><thead><tr><th></th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th></tr></thead><tbody><tr><td>Abaqus 3D</td><td>51.095</td><td>57.514</td><td>59.244</td><td>73.896</td><td>85.980</td><td>91.735</td></tr><tr><td>LD4</td><td>51.666</td><td>57.842</td><td>60.297</td><td>75.704</td><td>86.271</td><td>94.557</td></tr><tr><td>LD2</td><td>51.667</td><td>57.843</td><td>60.298</td><td>75.707</td><td>86.272</td><td>94.561</td></tr><tr><td>ED6</td><td>51.666</td><td>57.843</td><td>60.297</td><td>75.705</td><td>86.272</td><td>94.557</td></tr><tr><td>ED4</td><td>51.666</td><td>57.843</td><td>60.297</td><td>75.705</td><td>86.272</td><td>94.558</td></tr><tr><td>ED2</td><td>51.674</td><td>57.849</td><td>60.309</td><td>75.722</td><td>86.278</td><td>94.582</td></tr><tr><td>FSDT</td><td>51.606</td><td>57.714</td><td>60.280</td><td>75.711</td><td>86.084</td><td>94.573</td></tr></tbody></table>

It is possible to observe that the theories show a good approximation of the reference results. In this case, all the modes show radial displacements along the $\beta$-axis. The first three modes show six, five and seven semi-waves, respectively. The fourth, fifth and sixth modes show eight, four and nine semi-waves, respectively. The percentage error of FSDT is 3.09% for the sixth mode. For $l_{\alpha}/h = 10$, the first six modal frequencies are displayed in Table 12.
<table><caption>**Table 12.** Natural frequencies (10²· rad/s), $l_\alpha/h = 10$, case 3.</caption><thead><tr><th></th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th></tr></thead><tbody><tr><td>Abaqus 3D</td><td>18.724</td><td>19.391</td><td>24.738</td><td>25.494</td><td>32.710</td><td>35.354</td></tr><tr><td>LD4</td><td>18.789</td><td>19.489</td><td>24.917</td><td>25.555</td><td>32.978</td><td>35.473</td></tr><tr><td>LD2</td><td>18.795</td><td>19.501</td><td>24.937</td><td>25.558</td><td>33.009</td><td>35.475</td></tr><tr><td>ED6</td><td>18.789</td><td>19.489</td><td>24.918</td><td>25.555</td><td>32.979</td><td>35.475</td></tr><tr><td>ED4</td><td>18.790</td><td>19.491</td><td>24.921</td><td>25.557</td><td>32.984</td><td>35.480</td></tr><tr><td>ED2</td><td>18.836</td><td>19.562</td><td>25.046</td><td>25.595</td><td>33.182</td><td>35.541</td></tr><tr><td>FSDT</td><td>18.808</td><td>19.554</td><td>25.056</td><td>25.573</td><td>33.213</td><td>35.533</td></tr></tbody></table>

As previously noted, classical theories, and low-order theories in general, provide less accurate approximations of natural frequencies, particularly in cases with a low side-to-thickness ratio. Conversely, layer-wise theories offer the best approximation and align more closely with the Abaqus solution for higher frequencies. The modal shapes obtained through the ED4 model for $l_\alpha/h = 10$ are shown in Figures 11. These shapes correspond to those obtained in Abaqus, which are not included here for the sake of brevity.

<figure><img src="image_13.png" /><figcaption>**Figure 11.** First six eigenmodes, $l_\alpha/h = 10$, case 3.</figcaption></figure>

Table 13 shows the percentage of effective masses for the first six modal frequencies, considering $l_\alpha/h = 10$.

In this case, the predominant mode is the sixth, which exhibits significant displacement of the structural mass along the Z-axis. Table 14 provides the first six modal frequencies for $l_\alpha/h = 5$.
<table><caption>**Table 13.** Effective mass percentage $m_{\text{eff}\%i, j}$, $l_{\alpha}/h = 10$, case 3.</caption><thead><tr><th></th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th></tr></thead><tbody><tr><th>$m_{\text{eff}\%X}$</th><td>-*</td><td>-*</td><td>-*</td><td>-*</td><td>-*</td><td>0.3</td></tr><tr><th>$m_{\text{eff}\%Y}$</th><td>0.1</td><td>-*</td><td>-*</td><td>-*</td><td>-*</td><td>0.2</td></tr><tr><th>$m_{\text{eff}\%Z}$</th><td>-*</td><td>-*</td><td>-*</td><td>-*</td><td>-*</td><td>43.4</td></tr></tbody></table>

<table><caption>**Table 14.** Natural frequencies ($10^2 \cdot \text{rad/s}$), $l_{\alpha}/h = 5$, case 3.</caption><thead><tr><th></th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th></tr></thead><tbody><tr><th>Abaqus 3D</th><td>25.989</td><td>28.627</td><td>32.082</td><td>36.316</td><td>42.619</td><td>51.291</td></tr><tr><th>LD4</th><td>26.068</td><td>28.699</td><td>32.207</td><td>36.407</td><td>42.808</td><td>51.430</td></tr><tr><th>LD2</th><td>26.097</td><td>28.721</td><td>32.267</td><td>36.417</td><td>42.905</td><td>51.440</td></tr><tr><th>ED6</th><td>26.068</td><td>28.701</td><td>32.208</td><td>36.414</td><td>42.811</td><td>51.440</td></tr><tr><th>ED4</th><td>26.074</td><td>28.709</td><td>32.217</td><td>36.430</td><td>42.824</td><td>51.459</td></tr><tr><th>ED2</th><td>26.271</td><td>28.853</td><td>32.597</td><td>36.575</td><td>43.442</td><td>51.645</td></tr><tr><th>FSDT</th><td>26.261</td><td>28.862</td><td>32.631</td><td>36.572</td><td>43.524</td><td>51.620</td></tr></tbody></table>

Since a thick shell is considered, the effect of transverse stresses is not negligible, which causes the classical and lower-order theories to be less inaccurate. This can be observed for FSDT, which has an error as high as 2.12% for the fifth mode. This error can be reduced to 0.45% and 0.44% considering ED6 and LD4, respectively. In this case, the difference between ESL and LW models is more evident, since a multilayer shell is considered. Indeed, the possibility of obtaining a layer-wise representation of the laminate allows one to better predict the displacement field of each layer and, as a consequence, have more accurate values of the natural frequencies. The first mode is similar in both the $l_{\alpha}/h = 10$ and the current case. The second mode shows two semi-waves in the circumferential direction. The third mode exhibits four semi-waves along $\beta$. The fourth mode shows simple bending along the X-axis. The fifth mode displays five semi-waves along $\beta$. Finally, the sixth mode shows an extension of the structure on the Y-Z plane, which corresponds to a radial mode with zero circumferential semi-waves.

## 4. Conclusions

This paper presents a numerical methodology for the free vibration analysis of VAT shells. The proposed approach is based on the use of the CUF within an FE context, where the governing equations are derived through the PVD, considering displacements as primary variables. The main innovation is the adaptation of this framework to shell structures characterized by a mid-surface with non-constant first quadratic form coefficients. As can be observed in Equation (18), non-constant A and B coefficients generate new terms in the $A_p$ matrix, which has a direct influence on geometric relations. This also implies more complex terms in the FN of the stiffness matrix, since the derivatives $A_{r\beta}$, $B_{r\alpha}$, $H_{\alpha\beta}$ and $H_{\beta\alpha}$ need to be integrated. Considering variable coefficients of the first quadratic form enables the study of more complex geometries, laying the foundation for further research on variable curvature shells. Moreover, the formulation of the FN components is also applicable to shells with varying curvature radii. Three different cases of analysis are presented, and the results obtained through various CUF theories are validated with Abaqus 3D reference solutions. The variation in the fibers' angle is described through linear laws, but the presented formulation is general and can be adapted to more complex fiber distributions. ESL and LW approaches are both considered and their performance is compared to FSDT. The latter offers a good balance between accuracy and computational cost; however, it provides less accurate predictions of the behavior

<footnote>* Values below $10^{-1}$ are considered as negligible and marked with “-”.</footnote>
of thick shells. This is evident in the second analysis case, where FSDT demonstrates a significant loss of accuracy in predicting high-frequency modes when the side-to-thickness ratio is equal to five. On the other hand, ESL and LW theories yield the best match of the reference 3D solution, independently from the shell geometry or fiber variational law. Moreover, the differences between ESL and LW theories are more evident in cases where shells with multiple layers are considered. Indeed, for those cases, a layer-wise description of the shell kinematics allows one to obtain more accurate displacement fields and, as a consequence, modal frequencies. Considering the same expansion order N along the thickness, when multiple layers are taken into account, LW theories exhibit higher computational costs compared to ESL theories. Consequently, the choice between using LW and ESL models should be guided by a balance between accuracy and computational cost. In conclusion, the application of PVD within CUF has shown the possibility of modeling VAT shells for free vibration analyses within an FE context. One of the main drawbacks of the presented approach is the need to account for the mid-surface equation to accurately represent the shell geometry. Consequently, the analysis process depends on the analytical representation of the surface, increasing the complexity of its application to generic geometries. Additionally, the significant increase in DOF for LW models is not always justified, particularly for thin plates. In such cases, the FSDT model remains advantageous, offering a better balance between computational cost and result accuracy. However, the FSDT is incorporated into the presented formulation and can be utilized as a first-order ESL model. The potential of this approach is not limited to the shells presented in this work. Potential prospects foresee the use of this framework for studying more complex VAT shells, characterized by a non-constant curvature radius. Moreover, more complex variational statements, like the Reissner's mixed variational theorem, can be developed into this framework. This would enable the inclusion of transverse stresses as primary variables, further enhancing the accuracy of the solution. Additionally, the proposed approach can be extended to multi-field analyses, as the coupling between thermal, mechanical and electrical fields is of significant interest in various modern engineering applications, including aircraft and spacecraft design, aeroelasticity, wind turbine blades, actuators and sensors. In these contexts, additional advantages can be achieved by the ability to study surfaces with variable curvature radii and generic mid-surface geometries.

**Author Contributions:** Methodology, D.A.I., G.G., L.K. and M.M.; software, D.A.I.; validation, D.A.I.; writing—original draft, D.A.I.; writing—review and editing, G.G., L.K. 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:** Data will be made available on request.

**Conflicts of Interest:** The authors declare no conflicts of interest.
Appendix A

This appendix reports the FN of the structure stiffness and mass matrices. In the case of an orthotropic material, the components of the stiffness matrix FN can be written as follows:

$$
\begin{equation}
\begin{aligned}
K_{xx}^{\tau sij} = \int_{\Omega} & \left[ \left( \mathbf{N}_{m1}\mathbf{N}_{i_{\alpha}m1}\tilde{Z}_{pp11}^{\tau s} \frac{H_{\beta}}{H_{\alpha}} + \mathbf{N}_{m3}\mathbf{N}_{i_{\beta}m3}\tilde{Z}_{pp16}^{\tau s} \right) \mathbf{N}_{o1}\mathbf{N}_{j_{\alpha}o1} \right. \\
& \left. + \left( \mathbf{N}_{m1}\mathbf{N}_{i_{\alpha}m1}\tilde{Z}_{pp16}^{\tau s} + \mathbf{N}_{m3}\mathbf{N}_{i_{\beta}m3}\tilde{Z}_{pp66}^{\tau s} \frac{H_{\alpha}}{H_{\beta}} \right) \mathbf{N}_{o3}\mathbf{N}_{j_{\beta}o3} \right. \\
& \left. + \left( \mathbf{N}_{m1}\mathbf{N}_{i_{\alpha}m1}\tilde{Z}_{pp12}^{\tau s} \frac{B_{r\alpha}}{A} + \mathbf{N}_{m3}\mathbf{N}_{i_{\beta}m3}\tilde{Z}_{pp26}^{\tau s} \frac{B_{r\alpha} H_{\alpha}}{AH_{\beta}} \right) \mathbf{N}_{o2}\mathbf{N}_{j_{o2}} \right. \\
& \left. - \left( \mathbf{N}_{m1}\mathbf{N}_{i_{\alpha}m1}\tilde{Z}_{pp16}^{\tau s} \frac{H_{\alpha,\beta}}{H_{\alpha}} + \mathbf{N}_{m3}\mathbf{N}_{i_{\beta}m3}\tilde{Z}_{pp66}^{\tau s} \frac{H_{\alpha,\beta}}{H_{\beta}} \right) \mathbf{N}_{o3}\mathbf{N}_{j_{o3}} \right. \\
& \left. + \left( \mathbf{N}_{m2}\mathbf{N}_{im2}\tilde{Z}_{pp12}^{\tau s} \frac{B_{r\alpha}}{A} - \mathbf{N}_{m3}\mathbf{N}_{im3}\tilde{Z}_{pp16}^{\tau s} \frac{H_{\alpha,\beta}}{H_{\alpha}} \right) \mathbf{N}_{o1}\mathbf{N}_{j_{\alpha}o1} \right. \\
& \left. + \left( \mathbf{N}_{m2}\mathbf{N}_{im2}\tilde{Z}_{pp26}^{\tau s} \frac{B_{r\alpha} H_{\alpha}}{AH_{\beta}} - \mathbf{N}_{m3}\mathbf{N}_{im3}\tilde{Z}_{pp66}^{\tau s} \frac{H_{\alpha,\beta}}{H_{\beta}} \right) \mathbf{N}_{o3}\mathbf{N}_{j_{\beta}o3} \right. \\
& \left. + \left( \mathbf{N}_{m2}\mathbf{N}_{im2}\tilde{Z}_{pp22}^{\tau s} \frac{B_{r\alpha}^2 H_{\alpha}}{A^2 H_{\beta}} - \mathbf{N}_{m3}\mathbf{N}_{im3}\tilde{Z}_{pp26}^{\tau s} \frac{B_{r\alpha} H_{\alpha,\beta}}{AH_{\beta}} \right) \mathbf{N}_{o2}\mathbf{N}_{j_{o2}} \right. \\
& \left. - \left( \mathbf{N}_{m2}\mathbf{N}_{im2}\tilde{Z}_{pp26}^{\tau s} \frac{B_{r\alpha} H_{\alpha,\beta}}{AH_{\beta}} - \mathbf{N}_{m3}\mathbf{N}_{im3}\tilde{Z}_{pp66}^{\tau s} \frac{H_{\alpha,\beta}^2}{H_{\alpha}H_{\beta}} \right) \mathbf{N}_{o3}\mathbf{N}_{j_{o3}} \right. \\
& \left. + \mathbf{N}_{m1}\mathbf{N}_{im1}\tilde{Z}_{nn44}^{\tau s} \frac{A^2 H_{\beta}}{R_{\alpha}^2 H_{\alpha}} \mathbf{N}_{o1}\mathbf{N}_{j_{o1}} - \mathbf{N}_{m1}\mathbf{N}_{im1}\tilde{Z}_{nn44}^{\tau s_{z}} \frac{AH_{\beta}}{R_{\alpha}} \mathbf{N}_{o1}\mathbf{N}_{j_{o1}} - \mathbf{N}_{m1}\mathbf{N}_{im1}\tilde{Z}_{nn44}^{\tau z_{s}} \frac{AH_{\beta}}{R_{\alpha}} \mathbf{N}_{o1}\mathbf{N}_{j_{o1}} \right. \\
& \left. + \mathbf{N}_{m1}\mathbf{N}_{im1}\tilde{Z}_{nn44}^{\tau z_{s z}} H_{\alpha} H_{\beta} \mathbf{N}_{o1}\mathbf{N}_{j_{o1}} \right] d\alpha d\beta,
\end{aligned}
\tag{A1}
\end{equation}
$$

$$
\begin{equation}
\begin{aligned}
K_{xy}^{\tau sij} = \int_{\Omega} & \left[ \left( \mathbf{N}_{m3}\mathbf{N}_{im3}\tilde{Z}_{pp66}^{\tau s} \frac{H_{\alpha,\beta} H_{\beta,\alpha}}{H_{\alpha} H_{\beta}} - \mathbf{N}_{m2}\mathbf{N}_{im2}\tilde{Z}_{pp26}^{\tau s} \frac{B_{r\alpha} H_{\beta,\alpha}}{A H_{\beta}} \right) \mathbf{N}_{o3}\mathbf{N}_{j_{o3}} \right. \\
& \left. + \left( -\mathbf{N}_{m3}\mathbf{N}_{im3}\tilde{Z}_{pp16}^{\tau s} \frac{H_{\alpha,\beta} A_{r\beta}}{B H_{\alpha}} + \mathbf{N}_{m2}\mathbf{N}_{im2}\tilde{Z}_{pp12}^{\tau s} \frac{B_{r\alpha} A_{r\beta}}{A B} \right) \mathbf{N}_{o1}\mathbf{N}_{j_{o1}} \right. \\
& \left. - \left( \mathbf{N}_{m1}\mathbf{N}_{i_{\alpha}m1}\tilde{Z}_{pp16}^{\tau s} \frac{H_{\beta,\alpha}}{H_{\alpha}} + \mathbf{N}_{m3}\mathbf{N}_{i_{\beta}m3}\tilde{Z}_{pp66}^{\tau s} \frac{H_{\beta,\alpha}}{H_{\beta}} \right) \mathbf{N}_{o3}\mathbf{N}_{j_{o3}} \right. \\
& \left. + \left( \mathbf{N}_{m1}\mathbf{N}_{i_{\alpha}m1}\tilde{Z}_{pp11}^{\tau s} \frac{H_{\beta} A_{r\beta}}{H_{\alpha} B} + \mathbf{N}_{m3}\mathbf{N}_{i_{\beta}m3}\tilde{Z}_{pp16}^{\tau s} \frac{A_{r\beta}}{B} \right) \mathbf{N}_{o1}\mathbf{N}_{j_{o1}} \right. \\
& \left. + \left( -\mathbf{N}_{m3}\mathbf{N}_{im3}\tilde{Z}_{pp66}^{\tau s} \frac{H_{\alpha,\beta}}{H_{\alpha}} + \mathbf{N}_{m2}\mathbf{N}_{im2}\tilde{Z}_{pp26}^{\tau s} \frac{B_{r\alpha}}{A} \right) \mathbf{N}_{o3}\mathbf{N}_{j_{\alpha}o3} \right. \\
& \left. + \left( -\mathbf{N}_{m3}\mathbf{N}_{im3}\tilde{Z}_{pp26}^{\tau s} \frac{H_{\alpha,\beta}}{H_{\beta}} + \mathbf{N}_{m2}\mathbf{N}_{im2}\tilde{Z}_{pp22}^{\tau s} \frac{B_{r\alpha} H_{\alpha}}{A H_{\beta}} \right) \mathbf{N}_{o2}\mathbf{N}_{j_{\beta}o2} \right. \\
& \left. + \left( \mathbf{N}_{m1}\mathbf{N}_{i_{\alpha}m1}\tilde{Z}_{pp12}^{\tau s} + \mathbf{N}_{m3}\mathbf{N}_{i_{\beta}m3}\tilde{Z}_{pp26}^{\tau s} \frac{H_{\alpha}}{H_{\beta}} \right) \mathbf{N}_{o2}\mathbf{N}_{j_{\beta}o2} \right. \\
& \left. + \left( \mathbf{N}_{m1}\mathbf{N}_{i_{\alpha}m1}\tilde{Z}_{pp16}^{\tau s} \frac{H_{\beta}}{H_{\alpha}} + \mathbf{N}_{m3}\mathbf{N}_{i_{\beta}m3}\tilde{Z}_{pp66}^{\tau s} \right) \mathbf{N}_{o3}\mathbf{N}_{j_{\alpha}o3} \right. \\
& \left. + \mathbf{N}_{m1}\mathbf{N}_{im1}\tilde{Z}_{nn45}^{\tau s} \frac{A B}{R_{\alpha} R_{\beta}} \mathbf{N}_{o2}\mathbf{N}_{j_{o2}} - \mathbf{N}_{m1}\mathbf{N}_{im1}\tilde{Z}_{nn45}^{\tau s_{z}} \frac{A H_{\beta}}{R_{\alpha}} \mathbf{N}_{o2}\mathbf{N}_{j_{o2}} - \mathbf{N}_{m1}\mathbf{N}_{im1}\tilde{Z}_{nn45}^{\tau z_{s}} \frac{B H_{\alpha}}{R_{\beta}} \mathbf{N}_{o2}\mathbf{N}_{j_{o2}} \right. \\
& \left. + \mathbf{N}_{m1}\mathbf{N}_{im1}\tilde{Z}_{nn45}^{\tau z_{s z}} H_{\alpha} H_{\beta} \mathbf{N}_{o2}\mathbf{N}_{j_{o2}} \right] d\alpha d\beta,
\end{aligned}
\tag{A2}
\end{equation}
$$
$$
\begin{aligned}
K_{xz}^{\tau sij} = & \int_{\Omega} \left[ \left( -\mathbf{N}_{m3}\mathbf{N}_{im3}\tilde{Z}_{pp16}^{\tau s} \frac{AH_{\alpha,\beta}}{H_{\alpha}R_{\alpha}} + \mathbf{N}_{m2}\mathbf{N}_{im2}\tilde{Z}_{pp12}^{\tau s} \frac{AB_{\alpha,\alpha}}{AR_{\alpha}} \right) \mathbf{N}_{o1}\mathbf{N}_{jo1} \right. \\
& \left. + \left( -\mathbf{N}_{m3}\mathbf{N}_{im3}\tilde{Z}_{pp26}^{\tau s} \frac{H_{\alpha,\beta}B}{R_{\beta}H_{\beta}} + \mathbf{N}_{m2}\mathbf{N}_{im2}\tilde{Z}_{pp22}^{\tau s} \frac{BB_{\alpha,\alpha}H_{\alpha}}{AR_{\beta}H_{\beta}} \right) \mathbf{N}_{o2}\mathbf{N}_{jo2} \right. \\
& \left. + \left( \mathbf{N}_{m1}\mathbf{N}_{i,\alpha}m1\tilde{Z}_{pp11}^{\tau s} \frac{AH_{\beta}}{R_{\alpha}H_{\alpha}} + \mathbf{N}_{m3}\mathbf{N}_{i,\beta}m3\tilde{Z}_{pp16}^{\tau s} \frac{A}{R_{\alpha}} \right) \mathbf{N}_{o1}\mathbf{N}_{jo1} \right. \\
& \left. + \left( \mathbf{N}_{m1}\mathbf{N}_{i,\alpha}m1\tilde{Z}_{pp12}^{\tau s} \frac{B}{R_{\beta}} + \mathbf{N}_{m3}\mathbf{N}_{i,\beta}m3\tilde{Z}_{pp26}^{\tau s} \frac{BH_{\alpha}}{R_{\beta}H_{\beta}} \right) \mathbf{N}_{o2}\mathbf{N}_{jo2} \right. \\
& \left. + \left( -\mathbf{N}_{m3}\mathbf{N}_{im3}\tilde{Z}_{pn36}^{\tau s,z}H_{\alpha,\beta} + \mathbf{N}_{m2}\mathbf{N}_{im2}\tilde{Z}_{pn23}^{\tau s,z} \frac{B_{\alpha,\alpha}H_{\alpha}}{A} \right) \mathbf{N}_{j} \right. \\
& \left. + \left( \mathbf{N}_{m1}\mathbf{N}_{i,\alpha}m1\tilde{Z}_{pn13}^{\tau s,z}H_{\beta} + \mathbf{N}_{m3}\mathbf{N}_{i,\beta}m3\tilde{Z}_{pn36}^{\tau s,z}H_{\alpha} \right) \mathbf{N}_{j} \right. \\
& \left. - \mathbf{N}_{m1}\mathbf{N}_{im1}\tilde{Z}_{nn44}^{\tau s} \frac{AH_{\beta}}{R_{\alpha}H_{\alpha}} \mathbf{N}_{o1}\mathbf{N}_{j,\alpha}o1 - \mathbf{N}_{m1}\mathbf{N}_{im1}\tilde{Z}_{nn45}^{\tau s} \frac{A}{R_{\alpha}} \mathbf{N}_{o2}\mathbf{N}_{j,\beta}o2 + \mathbf{N}_{m1}\mathbf{N}_{im1}\tilde{Z}_{nn44}^{\tau s,z}H_{\beta}\mathbf{N}_{o1}\mathbf{N}_{j,\alpha}o1 \right. \\
& \left. + \mathbf{N}_{m1}\mathbf{N}_{im1}\tilde{Z}_{nn45}^{\tau s,z}H_{\alpha}\mathbf{N}_{o2}\mathbf{N}_{j,\beta}o2 \right] d\alpha d\beta,
\end{aligned}
\quad (A3)
$$

$$
\begin{aligned}
K_{yx}^{\tau sij} = & \int_{\Omega} \left[ \left( \mathbf{N}_{m3}\mathbf{N}_{im3}\tilde{Z}_{pp66}^{\tau s} \frac{H_{\alpha,\beta}H_{\beta,\alpha}}{H_{\alpha}H_{\beta}} - \mathbf{N}_{m1}\mathbf{N}_{im1}\tilde{Z}_{pp16}^{\tau s} \frac{A_{\alpha,\beta}H_{\alpha,\beta}}{BH_{\alpha}} \right) \mathbf{N}_{o3}\mathbf{N}_{jo3} \right. \\
& \left. + \left( -\mathbf{N}_{m3}\mathbf{N}_{im3}\tilde{Z}_{pp26}^{\tau s} \frac{H_{\beta,\alpha}B_{\alpha}}{AH_{\beta}} + \mathbf{N}_{m1}\mathbf{N}_{im1}\tilde{Z}_{pp12}^{\tau s} \frac{B_{\alpha,\alpha}A_{\beta}}{AB} \right) \mathbf{N}_{o2}\mathbf{N}_{jo2} \right. \\
& \left. - \left( \mathbf{N}_{m3}\mathbf{N}_{i,\alpha}m3\tilde{Z}_{pp66}^{\tau s} \frac{H_{\alpha,\beta}}{H_{\alpha}} + \mathbf{N}_{m2}\mathbf{N}_{i,\beta}m2\tilde{Z}_{pp26}^{\tau s} \frac{H_{\alpha,\beta}}{H_{\beta}} \right) \mathbf{N}_{o3}\mathbf{N}_{jo3} \right. \\
& \left. + \left( \mathbf{N}_{m3}\mathbf{N}_{i,\alpha}m3\tilde{Z}_{pp26}^{\tau s} \frac{B_{\alpha,\alpha}}{A} + \mathbf{N}_{m2}\mathbf{N}_{i,\beta}m2\tilde{Z}_{pp22}^{\tau s} \frac{B_{\alpha,\alpha}H_{\alpha}}{AH_{\beta}} \right) \mathbf{N}_{o2}\mathbf{N}_{jo2} \right. \\
& \left. + \left( -\mathbf{N}_{m3}\mathbf{N}_{im3}\tilde{Z}_{pp16}^{\tau s} \frac{H_{\beta,\alpha}}{H_{\alpha}} + \mathbf{N}_{m1}\mathbf{N}_{im1}\tilde{Z}_{pp11}^{\tau s} \frac{A_{\alpha,\beta}H_{\beta}}{BH_{\alpha}} \right) \mathbf{N}_{o1}\mathbf{N}_{j,\alpha}o1 \right. \\
& \left. + \left( -\mathbf{N}_{m3}\mathbf{N}_{im3}\tilde{Z}_{pp66}^{\tau s} \frac{H_{\beta,\alpha}}{H_{\beta}} + \mathbf{N}_{m1}\mathbf{N}_{im1}\tilde{Z}_{pp16}^{\tau s} \frac{A_{\alpha,\beta}}{B} \right) \mathbf{N}_{o3}\mathbf{N}_{j,\beta}o3 \right. \\
& \left. + \left( \mathbf{N}_{m3}\mathbf{N}_{i,\alpha}m3\tilde{Z}_{pp16}^{\tau s} \frac{H_{\beta}}{H_{\alpha}} + \mathbf{N}_{m2}\mathbf{N}_{i,\beta}m2\tilde{Z}_{pp12}^{\tau s} \right) \mathbf{N}_{o1}\mathbf{N}_{j,\alpha}o1 \right. \\
& \left. + \left( \mathbf{N}_{m3}\mathbf{N}_{i,\alpha}m1\tilde{Z}_{pp66}^{\tau s} + \mathbf{N}_{m2}\mathbf{N}_{i,\beta}m2\tilde{Z}_{pp26}^{\tau s} \frac{H_{\alpha}}{H_{\beta}} \right) \mathbf{N}_{o3}\mathbf{N}_{j,\beta}o3 \right. \\
& \left. + \mathbf{N}_{m2}\mathbf{N}_{im2}\tilde{Z}_{nn45}^{\tau s} \frac{AB}{R_{\alpha}R_{\beta}} \mathbf{N}_{o1}\mathbf{N}_{jo1} - \mathbf{N}_{m2}\mathbf{N}_{im2}\tilde{Z}_{nn45}^{\tau s,z} \frac{AH_{\beta}}{R_{\alpha}} \mathbf{N}_{o1}\mathbf{N}_{jo1} - \mathbf{N}_{m2}\mathbf{N}_{im2}\tilde{Z}_{nn45}^{\tau s,z} \frac{BH_{\alpha}}{R_{\beta}} \mathbf{N}_{o1}\mathbf{N}_{jo1} \right. \\
& \left. + \mathbf{N}_{m2}\mathbf{N}_{im2}\tilde{Z}_{nn45}^{\tau s,z}H_{\alpha}H_{\beta}\mathbf{N}_{o1}\mathbf{N}_{jo1} \right] d\alpha d\beta,
\end{aligned}
\quad (A4)
$$
$$
\begin{aligned}
K_{yy}^{\tau sij} = & \int_{\Omega} \left[ \left( \mathbf{N}_{m3}\mathbf{N}_{im3}\tilde{Z}_{pp66}^{\tau s} \frac{H_{\beta,\alpha}^2}{H_{\alpha}H_{\beta}} - \mathbf{N}_{m1}\mathbf{N}_{im1}\tilde{Z}_{pp16}^{\tau s} \frac{A_{\beta}H_{\beta,\alpha}}{BH_{\alpha}} \right) \mathbf{N}_{o3}\mathbf{N}_{jo3} \right. \\
& + \left( -\mathbf{N}_{m3}\mathbf{N}_{im3}\tilde{Z}_{pp16}^{\tau s} \frac{A_{\beta}H_{\beta,\alpha}}{BH_{\alpha}} + \mathbf{N}_{m1}\mathbf{N}_{im1}\tilde{Z}_{pp11}^{\tau s} \frac{A_{\beta}^2H_{\beta}}{B^2H_{\alpha}} \right) \mathbf{N}_{o1}\mathbf{N}_{jo1} \\
& - \left( \mathbf{N}_{m3}\mathbf{N}_{i,\alpha}m3\tilde{Z}_{pp66}^{\tau s} \frac{H_{\beta,\alpha}}{H_{\alpha}} + \mathbf{N}_{m2}\mathbf{N}_{i,\beta}m2\tilde{Z}_{pp26}^{\tau s} \frac{H_{\beta,\alpha}}{H_{\beta}} \right) \mathbf{N}_{o3}\mathbf{N}_{jo3} \\
& + \left( \mathbf{N}_{m3}\mathbf{N}_{i,\alpha}m3\tilde{Z}_{pp16}^{\tau s} \frac{A_{\beta}H_{\beta}}{BH_{\alpha}} + \mathbf{N}_{m2}\mathbf{N}_{i,\beta}m2\tilde{Z}_{pp12}^{\tau s} \frac{A_{\beta}}{B} \right) \mathbf{N}_{o1}\mathbf{N}_{jo1} \\
& + \left( -\mathbf{N}_{m3}\mathbf{N}_{im3}\tilde{Z}_{pp66}^{\tau s} \frac{H_{\beta,\alpha}}{H_{\alpha}} + \mathbf{N}_{m1}\mathbf{N}_{im1}\tilde{Z}_{pp16}^{\tau s} \frac{A_{\beta}H_{\beta}}{BH_{\alpha}} \right) \mathbf{N}_{o3}\mathbf{N}_{jo3} \\
& + \left( -\mathbf{N}_{m3}\mathbf{N}_{im3}\tilde{Z}_{pp26}^{\tau s} \frac{H_{\beta,\alpha}}{H_{\beta}} + \mathbf{N}_{m1}\mathbf{N}_{im1}\tilde{Z}_{pp12}^{\tau s} \frac{A_{\beta}}{B} \right) \mathbf{N}_{o2}\mathbf{N}_{j,\beta o2} \\
& + \left( \mathbf{N}_{m3}\mathbf{N}_{i,\alpha}m3\tilde{Z}_{pp66}^{\tau s} \frac{H_{\beta}}{H_{\alpha}} + \mathbf{N}_{m2}\mathbf{N}_{i,\beta}m2\tilde{Z}_{pp26}^{\tau s} \right) \mathbf{N}_{o3}\mathbf{N}_{j,\alpha o3} \\
& + \left( \mathbf{N}_{m3}\mathbf{N}_{i,\alpha}m3\tilde{Z}_{pp26}^{\tau s} + \mathbf{N}_{m2}\mathbf{N}_{i,\beta}m2\tilde{Z}_{pp22}^{\tau s} \frac{H_{\alpha}}{H_{\beta}} \right) \mathbf{N}_{o2}\mathbf{N}_{j,\beta o2} \\
& + \mathbf{N}_{m2}\mathbf{N}_{im2}\tilde{Z}_{nn55}^{\tau s} \frac{B^2H_{\alpha}}{R_{\beta}^2H_{\beta}} \mathbf{N}_{o2}\mathbf{N}_{jo2} - \mathbf{N}_{m2}\mathbf{N}_{im2}\tilde{Z}_{nn55}^{\tau s} \frac{BH_{\alpha}}{R_{\beta}} \mathbf{N}_{o2}\mathbf{N}_{jo2} - \mathbf{N}_{m2}\mathbf{N}_{im2}\tilde{Z}_{nn55}^{\tau s} \frac{BH_{\alpha}}{R_{\beta}} \mathbf{N}_{o2}\mathbf{N}_{jo2} \\
& + \mathbf{N}_{m2}\mathbf{N}_{im2}\tilde{Z}_{nn55}^{\tau s} H_{\alpha}H_{\beta} \mathbf{N}_{o2}\mathbf{N}_{jo2} \left. \right] d\alpha d\beta,
\end{aligned}
\quad (A5)
$$

$$
\begin{aligned}
K_{yz}^{\tau sij} = & \int_{\Omega} \left[ \left( -\mathbf{N}_{m3}\mathbf{N}_{im3}\tilde{Z}_{pp16}^{\tau s} \frac{AH_{\beta,\alpha}}{R_{\alpha}H_{\alpha}} + \mathbf{N}_{m1}\mathbf{N}_{im1}\tilde{Z}_{pp11}^{\tau s} \frac{AA_{\beta}H_{\beta}}{BR_{\alpha}H_{\alpha}} \right) \mathbf{N}_{o1}\mathbf{N}_{jo1} \right. \\
& + \left( -\mathbf{N}_{m3}\mathbf{N}_{im3}\tilde{Z}_{pp26}^{\tau s} \frac{BH_{\beta,\alpha}}{R_{\beta}H_{\beta}} + \mathbf{N}_{m1}\mathbf{N}_{im1}\tilde{Z}_{pp12}^{\tau s} \frac{A_{\beta}}{R_{\beta}} \right) \mathbf{N}_{o2}\mathbf{N}_{jo2} \\
& + \left( \mathbf{N}_{m3}\mathbf{N}_{i,\alpha}m3\tilde{Z}_{pp16}^{\tau s} \frac{AH_{\beta}}{R_{\alpha}H_{\alpha}} + \mathbf{N}_{m2}\mathbf{N}_{i,\beta}m2\tilde{Z}_{pp12}^{\tau s} \frac{A_{\beta}}{R_{\alpha}} \right) \mathbf{N}_{o1}\mathbf{N}_{jo1} \\
& + \left( \mathbf{N}_{m3}\mathbf{N}_{i,\alpha}m3\tilde{Z}_{pp26}^{\tau s} \frac{B}{R_{\beta}} + \mathbf{N}_{m2}\mathbf{N}_{i,\beta}m2\tilde{Z}_{pp22}^{\tau s} \frac{BH_{\alpha}}{R_{\beta}H_{\beta}} \right) \mathbf{N}_{o2}\mathbf{N}_{jo2} \\
& - \left( \mathbf{N}_{m3}\mathbf{N}_{im3}\tilde{Z}_{pn36}^{\tau s} H_{\beta,\alpha} + \mathbf{N}_{m1}\mathbf{N}_{im1}\tilde{Z}_{pn13}^{\tau s} \frac{A_{\beta}H_{\beta}}{B} \right) N_j \\
& + \left( \mathbf{N}_{m3}\mathbf{N}_{i,\alpha}m3\tilde{Z}_{pn36}^{\tau s} H_{\beta} + \mathbf{N}_{m2}\mathbf{N}_{i,\beta}m2\tilde{Z}_{pn23}^{\tau s} H_{\alpha} \right) N_j \\
& - \mathbf{N}_{m2}\mathbf{N}_{im2}\tilde{Z}_{nn45}^{\tau s} \frac{B}{R_{\beta}} \mathbf{N}_{o1}\mathbf{N}_{j,\alpha o1} - \mathbf{N}_{m2}\mathbf{N}_{im2}\tilde{Z}_{nn55}^{\tau s} \frac{BH_{\alpha}}{R_{\beta}H_{\beta}} \mathbf{N}_{o2}\mathbf{N}_{j,\beta o2} + \mathbf{N}_{m2}\mathbf{N}_{im2}\tilde{Z}_{nn45}^{\tau s} H_{\beta} \mathbf{N}_{o1}\mathbf{N}_{j,\alpha o1} \\
& \left. + \mathbf{N}_{m2}\mathbf{N}_{im2}\tilde{Z}_{nn55}^{\tau s} H_{\alpha} \mathbf{N}_{o2}\mathbf{N}_{j,\beta o2} \right] d\alpha d\beta,
\end{aligned}
\quad (A6)
$$
$$
\begin{aligned}
K_{zx}^{\tau sij} = & \int_{\Omega} \left[ - \left( \mathbf{N}_{m1}\mathbf{N}_{im1}\tilde{Z}_{pp16}^{\tau s} \frac{AH_{\alpha\beta}}{R_{\alpha}H_{\alpha}} + \mathbf{N}_{m2}\mathbf{N}_{im2}\tilde{Z}_{pp26}^{\tau s} \frac{BH_{\alpha\beta}}{R_{\beta}H_{\beta}} \right) \mathbf{N}_{o3}\mathbf{N}_{jo3} \right. \\
& \left. + \left( \mathbf{N}_{m1}\mathbf{N}_{im1}\tilde{Z}_{pp12}^{\tau s} \frac{B_{\alpha\alpha}}{R_{\alpha}} + \mathbf{N}_{m2}\mathbf{N}_{im2}\tilde{Z}_{pp22}^{\tau s} \frac{BB_{\alpha\alpha}H_{\alpha}}{AR_{\alpha}H_{\beta}} \right) \mathbf{N}_{o2}\mathbf{N}_{jo2} \right. \\
& \left. + \left( \mathbf{N}_{m1}\mathbf{N}_{im1}\tilde{Z}_{pp11}^{\tau s} \frac{AH_{\beta}}{R_{\alpha}H_{\alpha}} + \mathbf{N}_{m2}\mathbf{N}_{im2}\tilde{Z}_{pp12}^{\tau s} \frac{B}{R_{\beta}} \right) \mathbf{N}_{o1}\mathbf{N}_{j_{\alpha}o1} \right. \\
& \left. + \left( \mathbf{N}_{m1}\mathbf{N}_{im1}\tilde{Z}_{pp16}^{\tau s} \frac{A}{R_{\alpha}} + \mathbf{N}_{m2}\mathbf{N}_{im2}\tilde{Z}_{pp26}^{\tau s} \frac{BH_{\alpha}}{R_{\beta}H_{\beta}} \right) \mathbf{N}_{o3}\mathbf{N}_{j_{\beta}o3} \right. \\
& \left. - \left( \mathbf{N}_{m1}\mathbf{N}_{i_{\alpha}m1}\tilde{Z}_{nn44}^{\tau s} \frac{AH_{\beta}}{R_{\alpha}H_{\alpha}} + \mathbf{N}_{m2}\mathbf{N}_{i_{\beta}m2}\tilde{Z}_{nn45}^{\tau s} \frac{A}{R_{\alpha}} \right) \mathbf{N}_{o1}\mathbf{N}_{jo1} \right. \\
& \left. + \left( \mathbf{N}_{m1}\mathbf{N}_{i_{\alpha}m1}\tilde{Z}_{nn44}^{\tau s,z} H_{\beta} + \mathbf{N}_{m2}\mathbf{N}_{i_{\beta}m2}\tilde{Z}_{nn45}^{\tau s,z} H_{\alpha} \right) \mathbf{N}_{o1}\mathbf{N}_{jo1} \right. \\
& \left. - N_{i}\tilde{Z}_{pn36}^{\tau z s} H_{\alpha,\beta} \mathbf{N}_{o3}\mathbf{N}_{jo3} + N_{i}\tilde{Z}_{pn23}^{\tau z s} \frac{B_{\alpha\alpha}H_{\alpha}}{A} \mathbf{N}_{o2}\mathbf{N}_{jo2} + N_{i}\tilde{Z}_{pn13}^{\tau z s} H_{\beta} \mathbf{N}_{o1}\mathbf{N}_{j_{\alpha}o1} \right. \\
& \left. + N_{i}\tilde{Z}_{pn36}^{\tau z s} H_{\alpha} \mathbf{N}_{o3}\mathbf{N}_{j_{\beta}o3} \right] d\alpha d\beta,
\end{aligned}
\quad (A7)
$$

$$
\begin{aligned}
K_{zy}^{\tau sij} = & \int_{\Omega} \left[ - \left( \mathbf{N}_{m1}\mathbf{N}_{i_{\alpha}m1}\tilde{Z}_{nn45}^{\tau s} \frac{B}{R_{\beta}} + \mathbf{N}_{m2}\mathbf{N}_{i_{\beta}m2}\tilde{Z}_{nn55}^{\tau s} \frac{BH_{\alpha}}{R_{\beta}H_{\beta}} \right) \mathbf{N}_{o2}\mathbf{N}_{jo2} \right. \\
& \left. - \left( \mathbf{N}_{m1}\mathbf{N}_{im1}\tilde{Z}_{pp16}^{\tau s} \frac{AH_{\beta\alpha}}{R_{\alpha}H_{\alpha}} + \mathbf{N}_{m2}\mathbf{N}_{im2}\tilde{Z}_{pp26}^{\tau s} \frac{BH_{\beta\alpha}}{R_{\beta}H_{\beta}} \right) \mathbf{N}_{o3}\mathbf{N}_{jo3} \right. \\
& \left. + \left( \mathbf{N}_{m1}\mathbf{N}_{im1}\tilde{Z}_{pp11}^{\tau s} \frac{AA_{\beta}H_{\beta}}{BR_{\alpha}H_{\alpha}} + \mathbf{N}_{m2}\mathbf{N}_{im2}\tilde{Z}_{pp12}^{\tau s} \frac{A_{\beta}}{R_{\beta}} \right) \mathbf{N}_{o1}\mathbf{N}_{jo1} \right. \\
& \left. + \left( \mathbf{N}_{m1}\mathbf{N}_{im1}\tilde{Z}_{pp16}^{\tau s} \frac{AH_{\beta}}{R_{\alpha}H_{\alpha}} + \mathbf{N}_{m2}\mathbf{N}_{im2}\tilde{Z}_{pp26}^{\tau s} \frac{B}{R_{\beta}} \right) \mathbf{N}_{o3}\mathbf{N}_{j_{\alpha}o3} \right. \\
& \left. + \left( \mathbf{N}_{m1}\mathbf{N}_{im1}\tilde{Z}_{pp12}^{\tau s} \frac{A}{R_{\alpha}} + \mathbf{N}_{m2}\mathbf{N}_{im2}\tilde{Z}_{pp22}^{\tau s} \frac{BH_{\alpha}}{R_{\beta}H_{\beta}} \right) \mathbf{N}_{o2}\mathbf{N}_{j_{\beta}o2} \right. \\
& \left. + \left( \mathbf{N}_{m1}\mathbf{N}_{i_{\alpha}m1}\tilde{Z}_{nn45}^{\tau s,z} H_{\beta} + \mathbf{N}_{m2}\mathbf{N}_{i_{\beta}m2}\tilde{Z}_{nn55}^{\tau s,z} H_{\alpha} \right) \mathbf{N}_{o2}\mathbf{N}_{jo2} \right. \\
& \left. - N_{i}\tilde{Z}_{pn36}^{\tau z s} H_{\beta_{\alpha}} \mathbf{N}_{o3}\mathbf{N}_{jo3} + N_{i}\tilde{Z}_{pn13}^{\tau z s} \frac{A_{\beta}H_{\beta}}{B} \mathbf{N}_{o1}\mathbf{N}_{jo1} + N_{i}\tilde{Z}_{pn36}^{\tau z s} H_{\beta} \mathbf{N}_{o3}\mathbf{N}_{j_{\alpha}o3} \right. \\
& \left. + N_{i}\tilde{Z}_{pn23}^{\tau z s} H_{\alpha} \mathbf{N}_{o2}\mathbf{N}_{j_{\beta}o2} \right] d\alpha d\beta,
\end{aligned}
\quad (A8)
$$

$$
\begin{aligned}
K_{zz}^{\tau sij} = & \int_{\Omega} \left[ \left( \mathbf{N}_{m1}\mathbf{N}_{im1}\tilde{Z}_{pp11}^{\tau s} \frac{A^2 H_{\beta}}{R_{\alpha}^2 H_{\alpha}} + \mathbf{N}_{m2}\mathbf{N}_{im2}\tilde{Z}_{pp12}^{\tau s} \frac{AB}{R_{\alpha}R_{\beta}} \right) \mathbf{N}_{o1}\mathbf{N}_{jo1} \right. \\
& \left. + \left( \mathbf{N}_{m1}\mathbf{N}_{im1}\tilde{Z}_{pp12}^{\tau s} \frac{AB}{R_{\alpha}R_{\beta}} + \mathbf{N}_{m2}\mathbf{N}_{im2}\tilde{Z}_{pp22}^{\tau s} \frac{B^2 H_{\alpha}}{R_{\beta}^2 H_{\beta}} \right) \mathbf{N}_{o2}\mathbf{N}_{jo2} \right. \\
& \left. + \left( \mathbf{N}_{m1}\mathbf{N}_{i_{\alpha}m1}\tilde{Z}_{nn44}^{\tau s} \frac{H_{\beta}}{H_{\alpha}} + \mathbf{N}_{m2}\mathbf{N}_{i_{\beta}m2}\tilde{Z}_{nn45}^{\tau s} \right) \mathbf{N}_{o1}\mathbf{N}_{j_{\alpha}o1} \right. \\
& \left. + \left( \mathbf{N}_{m1}\mathbf{N}_{i_{\alpha}m1}\tilde{Z}_{nn45}^{\tau s} + \mathbf{N}_{m2}\mathbf{N}_{i_{\beta}m2}\tilde{Z}_{nn55}^{\tau s} \frac{H_{\alpha}}{H_{\beta}} \right) \mathbf{N}_{o2}\mathbf{N}_{j_{\beta}o2} \right. \\
& \left. + \left( \mathbf{N}_{m1}\mathbf{N}_{im1}\tilde{Z}_{pn13}^{\tau s,z} \frac{AH_{\beta}}{R_{\alpha}} + \mathbf{N}_{m2}\mathbf{N}_{im2}\tilde{Z}_{pn23}^{\tau s,z} \frac{BH_{\alpha}}{R_{\beta}} \right) N_j \right. \\
& \left. + N_i \tilde{Z}_{pn13}^{\tau z s} \frac{AH_{\beta}}{R_{\alpha}} \mathbf{N}_{o1}\mathbf{N}_{jo1} + N_i \tilde{Z}_{pn23}^{\tau z s} \frac{BH_{\alpha}}{R_{\beta}} \mathbf{N}_{o2}\mathbf{N}_{jo2} + N_i \tilde{Z}_{nn33}^{\tau z s,z} H_{\alpha} H_{\beta} N_j \right] d\alpha d\beta.
\end{aligned}
\quad (A9)
$$
The components of the mass matrix FN can be written as follows:

$$M_{xx}^{\tau sij} = M_{yy}^{\tau sij} = M_{zz}^{\tau sij} = \int_{\Omega} (N_i \rho E_{\tau s} H_{\alpha} H_{\beta} N_j) d\alpha d\beta, \quad (A10)$$

$$M_{xy}^{\tau sij} = M_{xz}^{\tau sij} = M_{yx}^{\tau sij} = M_{yz}^{\tau sij} = M_{zx}^{\tau sij} = M_{zy}^{\tau sij} = 0. \quad (A11)$$

## References

1. Brasington, A.; Sacco, C.; Halbritter, J.; Wehbe, R.; Harik, R. Automated fiber placement: A review of history, current technologies, and future paths forward. *Compos. Part C Open Access* **2021**, 6, 100182. [CrossRef]

2. Kim, C.; Potter, K.; Weaver, M. Continuous tow shearing for manufacturing variable angle tow composites. *Compos. Part A Appl. Sci. Manuf.* **2012**, 43, 1347–1356. [CrossRef]

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

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

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

6. Love, A.E.H. XVI. The small free vibrations and deformation of a thin elastic shell. *Philos. Trans. R. Soc. A* **1888**, 179, 491–546.

7. Donnell, L.H. *Stability of Thin-Walled Tubes Under Torsion*; NASA: Washington, DC, USA, 1933.

8. Flügge, W. *Stresses in Shells*; Springer Science & Business Media: Berlin/Heidelberg, Germany, 2013.

9. Lindberg, G.M.; Olson, M.D.; Cowper, G.R. New developments in the finite element analysis of shells. *Q. Bull. Div. Mech. Eng. Natl. Aeronaut. Establ.* **1969**, 4, 1–38.

10. Clough, R.W.; Johnson, C.P. A finite element approximation for the analysis of thin shells. *Int. J. Solids* **1969**, 4, 43–60. [CrossRef]

11. Bogner, F.K.; Fox, R.L.; Schmit, L.A. A cylindrical shell discrete element. *AIAA J.* **1967**, 5, 745–750. [CrossRef]

12. Maksymyuk, V.A. Locking phenomenon in computational methods of the shell theory. *Int. J. Appl. Mech.* **2020**, 56, 347–350. [CrossRef]

13. Bathe, K.-I.; Dvorkin, E.N. A formulation of general shell elements—The use of mixed interpolation of tensorial components. *Int. J. Numer. Methods Eng.* **1986**, 22, 697–722. [CrossRef]

14. Bathe, K.-J.; Bucalem, M.L. Higher-order MITC general shell elements. *Int. J. Numer. Methods Eng.* **1993**, 36, 3729–3754.

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

16. Carrera, E.; Campisi, C.; Cinefra, M.; Soave, M. Evaluation of various trough the thickness and curvature approximations for free vibrational analysis of cylindrical and spherical shells. *Int. J. Signal Imaging Syst. Eng.* **2008**, 1, 197.

17. Giunta, G.; Biscani, F.; Belouettar, S.; Carrera, E. Hierarchical modelling of doubly curved laminated composite shells under distributed and localised loadings. *Compos. B Eng.* **2011**, 42, 682–691. [CrossRef]

18. Giunta, G.; Melta, N.; Belouettar, S.; Ferreira, A.J.M.; Carrera, E. A thermo-mechanical analysis of isotropic and composite beams via collocation with radial basis functions. *J. Therm. Stresses* **2013**, 36, 1169–1199. [CrossRef]

19. Brischetto, S.; Polit, O.; Carrera, E. Refined shell model for the linear analysis of isotropic and composite elastic structures. *Eur. J. Mech. A-Solids* **2012**, 34, 102–119. [CrossRef]

20. Cinefra, M.; Carrera, E. Shell finite elements with different through-the-thickness kinematics for the linear analysis of cylindrical multilayered structures. *Int. J. Numer. Meth. Eng.* **2013**, 93, 160–182. [CrossRef]

21. Cinefra, M.; Chinosi, C.; Della Croce, L. MITC9 shell elements based on refined theories for the analysis of isotropic cylindrical structures. *Int. J. Numer. Meth. Eng.* **2013**, 20, 91–100. [CrossRef]

22. Cinefra, M.; Kumar, S.K.; Carrera, E. MITC9 Shell elements based on RMVT and CUF for the analysis of laminated composite plates and shells. *Compos. Struct.* **2019**, 209, 383–390. [CrossRef]

23. Kumar, S.K.; Harursampath, D.; Carrera, E.; Cinefra, M.; Valvano, S. Modal analysis of delaminated plates and shells using Carrera Unified Formulation–MITC9 shell element. *Mech. Adv. Mater. Struct.* **2018**, 25, 681–697. [CrossRef]

24. Carrera, E.; Zozulya, V.V. Carrera unified formulation (CUF) for shells of revolution. I. Higher-order theory. *Acta Mech.* **2023**, 234, 109–136. [CrossRef]

25. Carrera, E.; Zozulya, V.V. Carrera unified formulation (CUF) for the shells of revolution. II. Navier close form solutions. *Acta Mech.* **2023**, 234, 137–161. [CrossRef]

26. Carrera, E.; Zozulya, V.V. Carrera unified formulation (CUF) for the shells of revolution. Numerical evaluation. *Mech. Adv. Mater. Struct.* **2024**, 31, 1597–1619. [CrossRef]
27. Tornabene, F.; Liverani, A.; Caligiana, G. General anisotropic doubly-curved shell theory: A differential quadrature solution for free vibrations of shells and panels of revolution with a free-form meridian. *J. Sound Vib.* **2012**, *331*, 4848–4869. [CrossRef]

28. Tornabene, F.; Fantuzzi, N.; Viola, E.; Carrera, E. Static analysis of doubly-curved anisotropic shells and panels using CUF approach, differential geometry and differential quadrature method. *Compos. Struct.* **2014**, *107*, 675–697. [CrossRef]

29. Wu, K.; Tatting, B.; Smith, B.; Stevens, R.; Occhipinti, G.; Swift, J.; Achary, D.; Thornburgh, R. Design and manufacturing of tow-steered composite shells using fiber placement. In Proceedings of the 50th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, Palm Springs, CA, USA, 4–7 May 2009; p. 2700.

30. Hyer, M.W.; Charette, R.F. *Innovative Design of Composite Structures: The Use of Curvilinear Fiber Format in Composite Structure Design*; NASA: Washington, DC, USA, 1990.

31. Tornabene, F.; Fantuzzi, N.; Baciocchi, M. Higher-order structural theories for the static analysis of doubly-curved laminated composite panels reinforced by curvilinear fibers. *Thin-Walled Struct.* **2016**, *102*, 222–245. [CrossRef]

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

33. Sciascia, G.; Oliveri, V.; Weaver, P.M. Eigenfrequencies of prestressed variable stiffness composite shells. *Compos. Struct.* **2021**, *270*, 114019. [CrossRef]

34. He, T.; Huang, W.; Kuang, Z.; Huang, Q.; Hou, R.; Giunta, G.; Belouettar, S.; Hu, H. A seven-parameter high-order finite element model for multi-stable analysis of variable stiffness laminated shells. *Thin-Walled Struct.* **2023**, *190*, 110983. [CrossRef]

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

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

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

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

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

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

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

42. Leissa, A.W. *Vibration of Shells*; Scientific and Technical Information Office, National Aeronautics and Space Administration: Washington, DC, USA, 1973.

43. Kraus, H. *Thin Elastic Shells: An Introduction to the Theoretical Foundations and the Analysis of Their Static and Dynamic Behavior*; John Wiley & Sons: Hoboken, NJ, USA, 1967.

44. Washizu, K. *Variational Methods in Elasticity and Plasticity*; Pergamon Press: Oxford, UK, 1974.

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

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

47. Bathe, K.-J. *Finite Element Procedures*; Prentice Hall: Hoboken, NJ, USA, 2006.

48. ABAQUS. *Abaqus Theory Manual*; Dassault Systèmes: Vélizy-Villacoublay, France, 2003.

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

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