## Abstract

This paper deals with the dynamic response of rotors made of anisotropic, laminated composite materials. It is a sequel to the authors’ previous work, which was devoted to the rotordynamics of metallic structures. The used variable kinematic one-dimensional models describe any cross-sectional deformation of the rotor and go beyond the plane strain assumptions of classical Euler–Bernoulli and Timoskenko beam theories. Refined theories are obtained by applying the Carrera unified formulation, which is extended here to the rotordynamics of multilayered composites. The displacement variables over the rotor cross section x-z plane are approximated by x,z polynomials of any order N. Thin-walled cylindrical shafts and boxes are analyzed. These structures are made of unidirectional layers, whose fiber orientation can vary with respect to the rotor–axis as well as in the x-z plane. Several analyses have been carried out to determine the vibrational response as a function of the rotating speed. Classical beam theories are obtained as particular cases and results available in the literature, including shell results, are used to assess the presented theory. The proposed refined models are very effective in investigating the dynamic behavior of laminated composite rotors.

## Introduction

It is well known that composite materials present excellent mechanical properties, such as high specific stiffness and strength, ease of formability, a wide range of operating temperatures, and many others. These properties justify their extensive use in many applications, among which the design of rotors, whose dynamic behavior is worthy of study. For example, interesting experiments were carried out in Refs. [1–3] on either graphite or boron/epoxy shafts. The advantages of orthotropic materials over their isotropic counterparts were pointed out, and useful references were provided for analytical formulations. Despite their limitations, classical beam theories have been extensively used to investigate the critical speeds and instabilities of composite shafts. For instance, with the purpose of proposing an optimization algorithm to define the best lamination scheme, a first shear deformation theory, including the rotatory inertia, was used in Ref. [4]. The optimization was direct in order to maximize the first bending frequency, thus ensuring a sufficient torsional strength by imposing that the lamination angle of a certain number of plies was equal to 45 deg. Later, Bert and Kim developed one-dimensional models, based on Euler–Bernoulli [5] and Bresse–Timoshenko [6] theories, which included bending–twisting coupling. The critical speeds were in good agreement with the results obtained through both shell and experimental approaches. Naturally, when shear effects become important, the Euler–Bernoulli theory does not lead to accurate results. In Ref. [7], Chen and Peng investigated the stability of composite spinning cylinders subjected to compressive loads using the Timoshenko model and the equivalent modulus beam theory (EMBT). The critical speeds and the effects due to a disk location were also studied under these assumptions. A further contribution was made by Chang et al. [8], who developed a simple first-order shear deformation theory for shafts supported by bearings modeled as springs and viscous dampers. The constitutive relations were derived directly from those of three-dimensional continua by means of coordinate transformations.

In order to overcome the limitations of the equivalent modulus approach, Singh and Gupta [9] first proposed a refined beam formulation derived from a layerwise shell theory (LBT) and then, in Ref. [10], presented a modified EMBT that included bending–twisting, shear–normal, and bending–stretching coupling effects. The obtained results were used by Sino et al. [11] to assess their simplified homogenized beam theory (SHBT). In the SHBT, the stiffness parameters are evaluated using an energy formulation that takes into account several contributions, such as Young's modulus, the shear modulus, the distance from the shaft axis, and the thickness of each rotor layer. The results were in good agreement with those presented by Gupta et al. Further attempts have also been made to improve displacement theories. Among these formulations, Librescu et al. provided a higher-order model that incorporates the warping and torsion of thin-walled anisotropic spinning structures. Critical speeds as well as stability were investigated considering thin-walled boxes [12,13] and cylinders [14].

When the assumptions of beam theories are too restrictive, two- and three-dimensional formulations become useful. For instance, in Ref. [15], Kim and Bert compared a number of shell theories (Loo's, Morley's, Love's, Donnell's, and Sanders's theories) to determine the critical speeds of the shaft examined in Ref. [1]. All these models yielded accurate results, except Donnell's theory, which was not effective for long shafts. A further example is Ref. [16], in which Ramezani and Ahmadian combined the layerwise theory and the wave propagation approach for rotating laminated shells under different boundary conditions. Interesting results were presented in Ref. [17], in which the dynamics of rotating cross-ply laminated cylinders was considered. The structure was reinforced with stringers and rings whose stiffness was considered either separately or smeared over the shell surface. As expected, the number of stiffeners had a remarkable effect on the backward and forward frequencies. Moreover, in Ref. [18], a disk–shaft assembly was discretized with shell finite elements, considering a cyclic symmetry of the structure; three analyses were carried out with different numbers of degrees of freedom. This assumption made it possible to limit the computational cost to a great extent. Finally, a comprehensive overview of rotordynamics phenomena with shell approaches can be found in the book by Hua et al. [19].

The present paper has the aim of introducing and comparing a variety of refined one-dimensional beam models in order to study the dynamics of composite rotors. The cross sections of the structures are built with a number of small orthotropic material plates, which are rotated around the rotor–axis of the β angle (see Fig. 1). This approach makes it possible to consider both thin-walled and compact cross sections with symmetric and asymmetric lamination schemes. As far as the thin-walled structures are concerned (cylinders and boxes), higher-order theories are able to detect deformations of the section, which may have an important effect on the overall dynamic behavior. These theories are derived from Carrera's unified formulation (CUF). This procedure makes it possible to consider the order and the type of models as free input parameters. CUF was originally developed for plates and shells [20–22] and was later extended to the beam model [23]. Using one-dimensional finite elements, static [24] as well as free vibration [25] analyses have already been carried out on laminated and sandwich structures, thus demonstrating that refined displacement fields lead to accurate results with a very low computational cost. Recently, Carrera et al. have extended the unified formulation to the study of spinning metallic rotors with thin disks [26] and centrifugally stiffened beams [27].

Fig. 1
Fig. 1
Close modal

In light of the encouraging results, the equations of motion (EOM) of spinning composite structures have been written for the first time. The displacement fields are obtained by adopting Taylor-like expansions. The EOM in the rotating frame are solved using the finite element method. In order to assess the new formulation, several linear analyses have been carried out and the related results have been compared with the solutions presented in the literature.

## The Unified Formulation

The CUF states that the displacement field $u(x,y,z,t)$ is an expansion of generic functions $Fτ(x,z)$ for the vector displacement $uτ(y)$
$u(x,y,z,t)=Fτ(x,z)uτ(y,t) τ=1,2,…,T$
(1)
where T is the number of terms of the expansion and, according to the generalized Einstein's notation, $τ$ indicates summation. In this work, Eq. (1) consists of Taylor's polynomials that are functions of the coordinates of the cross section. For example, the second-order displacement field is
$ux=ux1+x ux2+z ux3+x2 ux4+xz ux5+z2 ux6uy=uy1+x uy2+z uy3+x2 uy4+xz uy5+z2 uy6uz=uz1+x uz2+z uz3+x2 uz4+xz uz5+z2 uz6$
(2)
while the third-order displacement field becomes
$ux=ux1+x ux2+z ux3+x2 ux4+xz ux5+z2 ux6+x3 ux7+x2z ux8+xz2 ux9+z3 ux10uy=uy1+x uy2+z uy3+x2 uy4+xz uy5+z2 uy6+x3 uy7+x2 zuy8+xz2 uy9+z3 uy10uz=uz1+x uz2+z uz3+x2 uz4+xz uz5+z2 uz6+x3 uz7+x2 zuz8+xz2 uz9+z3 uz10$
(3)
A remarkable feature is that classical beam theories are obtainable as particular cases of Taylor expansions. It should be noted that classical theories require reduced material stiffness coefficients to contrast Poisson's locking. Unless otherwise specified, for classical and first-order models Poisson's locking is corrected according to Carrera and Giunta [23]. The stresses and the strains are grouped as follows:
$εp={ɛzz ɛxx ɛxz} T σp={σzz σxx σxz} Tεn={ɛzy ɛxy ɛyy} T σn={σzy σxy σyy} T$
(4)
The subscript “p” stands for terms lying on the cross section, while “n” stands for terms lying on the other planes, which are orthogonal to the cross section. The strain–displacement relations and the Hooke's law are, respectively,
$εp=Dpuεn=(Dny+Dnp)u$
(5)
$σp=C˜ppεp+C˜pnεnσn=C˜npεp+C˜nnεn$
(6)
Both box beams and cylinders can be considered constituted by a certain number of either straight or curved plates of orthotropic material, whose material coordinate systems (1,2,3) generally do not coincide with the physical coordinate system (x,y,z) (see Fig. 1). Using this approach, the matrices of material coefficients of the generic material k are
$C˜ppk=[C˜11kC˜12kC˜14kC˜12kC˜22kC˜24kC˜14kC˜24kC˜44k], C˜pnk=[C˜15kC˜16kC˜13kC˜25kC˜26kC˜23kC˜45kC˜46kC˜43k],C˜nnk=[C˜55kC˜56kC˜35kC˜56kC˜66kC˜36kC˜35kC˜36kC˜33k]$
(7)
For the sake of completeness, the explicit forms of the coefficients of the matrices $C˜$ are reported in Appendix  A. If a classical finite element technique is adopted with the purpose of easily dealing with arbitrary shaped cross sections, the generalized displacement vector becomes
$uτ(y,t)=Ni(y)qτi(t)$
(8)
where $Ni(y)$ are the shape functions and $qτi(t)$ is the nodal displacement vector
$qτi(t)={quxτiquyτiquzτi}T$
(9)

## Rotordynamics Equations in CUF Form

In order to obtain the motion equations of a structure that is rotating about its longitudinal axis with a constant speed Ω, the extended Hamilton's principle is invoked
$δ∫t0t1(T-U+Wb) dt=0$
(10)
where T and U are the kinetic and the potential energies in the rotating reference frame, which were previously derived in Ref. [26]. The term Wb is the contribution due to the Nb bearings that are modeled as springs and viscous dampers, whose expression is
$Wb=∑p=1NbuTKbpu+u·TCbpu$
(11)
Considering the generic pth bearing whose coordinates are (xb, yb, zb), the stiffness and damping coefficients are
$Kbp=[kxxpkxypkxzpkyxpkyypkyzpkzxpkzypkzzp], Cbp=[cxxpcxypcxzpcyxpcyypcyzpczxpczypczzp]$
(12)
Introducing the assumptions (1), (8), and applying the Euler–Lagrange equations, Eq. (10) becomes
$∫t0t1[δqτiTMijτsq··sj+δqτiT(GΩijτs+Gbijτs)q·sj + δqτiT(Kijτs+Kbijτs-KΩijτs)qsj+δqτiTFΩiτ r]dt=0$
(13)
The matrices written in terms of fundamental nuclei are
$Mijτs=Ilij◃(FτρkIFs)▹$
$GΩijτs=Ilij◃(FτρkIFs)▹2Ω$
$Gbijτs=∑p=1NbNi(ybp)Fτ(xbp,zbp)CbpFs(xbp,zbp)Nj(ybp)$
$Kijτs=Ilij◃DnpT(FτI)[C˜npkDp(FsI)+C˜nnkDnp(FsI)] +DpT(FτI)[C˜ppkDp(FsI)+C˜pnkDnp(FsI)]▹ +Ilij,y◃[DnpT(FτI)+DpT(FτI)C˜pnk]Fs▹+ IAy +Ili,yjIAyT◃Fτ[C˜npkDp(FsI)+C˜nnkDnp(FsI)]▹ +Ili,yj,yIAyTIAy◃FτC˜nnkFs▹$
(14)
$Kbijτs=∑p=1NbNi(ybp)Fτ(xbp,zbp)KbpFs(xbp,zbp)Nj(ybp)$
$KΩijτs=Ilij◃(FτρkIFs)▹ΩTΩ$
$FΩiτ=Ili◃Fτρr▹ΩTΩ$
where
$Ω=[00Ω000-Ω00]IAy=[001100010]◃…▹=∫A…dA$
(15)
$(Ili,Ilij,Ilij, y,Ili, yj,Ili, yj, y)=∫l(Ni,NiNj,NiNj, y,Ni, yNj,Ni, yNj, y)dy$
(16)

and $r={xP,0,zP}$ is the distance of a generic point P belonging to the cross section from the neutral axis. For sake of clearness, in Appendix  B, the nine components of the fundamental nucleus of the matrix $Kijτs$ are written in an explicit form.

To obtain the natural frequencies and the normal modes of the rotor, the homogeneous equations are solved assuming a periodic solution $q=qe¯iωt$:
$q¯eiωt[(K+Kp-KΩ)+(GΩ+Gp)iω-(M)ω2]=0$
(17)
The quadratic eigenvalues problem (QEP) of generic order R in Eq. (17) is transformed in a classical linear system of order $2×R$
${Mq··+(GΩ+Gp)q·+(K+Kp-KΩ)q=0-q·+q·=0$
(18)
now, by introducing
$a={qq·} a·={q·q··}$
(19)
the equations of motion assume the following form:
$RT-1iωI=0$
(20)
where
$T-1R=[(K+Kp-KΩ) -1(GΩ+Gp)(K+Kp-KΩ)-1M-I0]$
(21)

The problem in Eq. (20) is in the classical form and it can be solved with the standard eigensolvers.

## Numerical Results

### Cylindrical Composite Shaft.

This section shows the results related to analyses carried out on composite thin-walled structures. The first case concerned the dynamic study of hollow rotating cylinders whose material and geometrical features are summarized in Table 1. The shaft was discretized with seven four-node beam elements and it was supported by one bearing at each end, whose stiffness kxx and kzz was 1740 GN/m. The first critical speeds obtained with the present theories are shown in Table 2 and compared with those found in Ref. [8], adopting the $[90 deg/45 deg/-45 deg/0 deg6/90 deg]$ lamination scheme. The results for both materials are within the interval of reference speeds, and the differences between the three theories (TE2, TE3, and TE4) are essentially negligible. Furthermore, the effects of the aspect ratio for the graphite–epoxy shaft have also been investigated in Table 3. As expected, when the structure becomes thicker and thicker, the critical speeds occur at higher values, and the differences between the three expansions increase; in fact, when the aspect ratio is equal to 2, TE3 and TE4 lead to similar results, whereas TE2 overestimates the critical speed of 10%. In addition, considering the same material, the effects due to the lamination scheme have been evaluated and the results are shown in Table 4. Increasing the lamination angle, with respect to the longitudinal direction, reduces the shaft stiffness and causes a consequent reduction in the critical speeds. In the case in point, some evident discrepancies can be observed between the refined beam models and those found in literature. In fact, when θ is equal to 15 deg, 30 deg, or 45 deg, the present results are closer to speeds obtained with the shell theory SST than those of classical beam models. In order to show the accuracy of the CUF elements, natural frequencies of the boron–epoxy shaft (Table 1), computed with TE elements, are listed in Table 5 and compared with the results found in Ref. [28], where a first-order shear beam theory with a torsion function was conceived, and with those obtained using Di.Qu.M.A.S.P.A.B. [29], a free software based on the differential quadrature method for shells and plates. At least the third-order theory (TE3) is necessary for all the considered cases to obtain accurate results for bending as well as torsional modes. Moreover, it is worth noting that remarkable differences exist between the first natural frequencies of Ref. [28] and those obtained with the other approaches, especially for θ equal to 15 deg, 30 deg, 45 deg, and 60 deg. The use of a constant shear correction factor could be the reason for these important discrepancies, since the shear effects are closely related to the lamination scheme.

Table 1

Material and geometrical features of composite shaft

DimensionsCylindersBox beam
Thickness (h), m0.0013210.01016
Width (c), m0.1269a0.1016
Length (L), m2.4701.016
MaterialsBoron–epoxy [8]Graphite–epoxy [8]
E11, GPa211.0139.0
E22, GPa24.111.0
G23, GPa6.903.78
G31, GPa6.906.05
G12, GPa6.906.05
ρ, kg m81967.01578.0
ν120.3600.313
DimensionsCylindersBox beam
Thickness (h), m0.0013210.01016
Width (c), m0.1269a0.1016
Length (L), m2.4701.016
MaterialsBoron–epoxy [8]Graphite–epoxy [8]
E11, GPa211.0139.0
E22, GPa24.111.0
G23, GPa6.903.78
G31, GPa6.906.05
G12, GPa6.906.05
ρ, kg m81967.01578.0
ν120.3600.313
a

Mean diameter.

Table 2

Critical speeds (rpm) of composite shafts for [90 deg/45 deg/−45 deg/0 deg6/90 deg]

TheoryaBoron–epoxyGraphite–epoxy
SST [8]58725349
EBBT [8]59195302
BTBT [8]57885113
FSDT [8]57625197
TE258085256
TE357665232
TE457245220
TheoryaBoron–epoxyGraphite–epoxy
SST [8]58725349
EBBT [8]59195302
BTBT [8]57885113
FSDT [8]57625197
TE258085256
TE357665232
TE457245220
a

SST: Sanders shell theory. EBBT: Euler–Bernoulli beam theory. BTBT: Bresse–Timoshenko beam theory. FSDT: First-order shear deformation theory.

Table 3

Critical speeds (rpm) of graphite–epoxy shaft for several length-to-diameter ratios

L/D
Theory2515
SST [8]112,40041,6808585
EBBT [8]329,60076,8209072
BTBT [8]176,30054,8308543
FSDT [8]176,30055,7068527
TE2240,90065,5808760
TE3206,40061,2008580
TE4206,40061,2008580

L/D
Theory2515
SST [8]112,40041,6808585
EBBT [8]329,60076,8209072
BTBT [8]176,30054,8308543
FSDT [8]176,30055,7068527
TE2240,90065,5808760
TE3206,40061,2008580
TE4206,40061,2008580
Table 4

Critical speeds (rpm) of composite shaft for various lamination angles

Lamination
Theory0 deg15 deg30 deg45 deg60 deg75 deg90 deg
SST [8]5527436533082386212020201997
EBBT [8]6425539342693171229218851813
FSDT [8]6072533142063124228418901816
TE16270477031202280195018601860
TE26270471031202280195018601860
TE36060417027602160193518601830
TE46060414027302145192018301800

Lamination
Theory0 deg15 deg30 deg45 deg60 deg75 deg90 deg
SST [8]5527436533082386212020201997
EBBT [8]6425539342693171229218851813
FSDT [8]6072533142063124228418901816
TE16270477031202280195018601860
TE26270471031202280195018601860
TE36060417027602160193518601830
TE46060414027302145192018301800
Table 5

Frequencies at standstill of boron–epoxy shaft

Lam.ModeFSDT [28]FSDT [29]TE4TE3TE2TE1aFSDTaEBBTa
0 deg1st deg flex.108.23109.23109.24109.24113.95113.95113.95119.44
2nd deg flex.357.13357.15357.19403.38403.37403.37475.41
1st deg tors.379.13379.14379.14379.14379.14
15 deg1st deg flex.94.8669.5769.8370.0581.7381.7681.76108.39
2nd deg flex.261.23262.15262.81309.33309.47309.47431.46
1st deg tors.412.23412.79412.79413.60412.60
30 deg1st deg flex.76.7145.8246.0546.6052.4152.3352.3374.39
2nd deg flex.180.30181.11183.10205.85205.16205.16296.12
1st deg tors.519.80519.98519.98520.63519.80
45 deg1st deg flex.59.2137.8537.9838.4839.5739.5439.5443.62
2nd deg flex.150.13150.55152.51156.95156.30156.30173.62
1st deg tors.638.00639.43639.44691.04663.24
60 deg1st deg flex.45.8436.4936.5936.7436.7536.7536.7536.89
2nd deg flex.144.11144.30145.10145.35145.14145.14146.84
1st deg tors.534.03540.09540.10675.61586.18
75 deg1st deg flex.40.7438.3338.4638.5838.9939.0039.0039.58
2nd deg flex.149.60149.75150.42153.16153.21153.21157.54
1st deg tors.415.73419.76419.71460.93429.83
90 deg1st deg flex.40.1039.8039.9540.4740.7040.7940.7941.03
2nd deg flex.154.05154.18156.07159.34159.71159.71163.33
1st deg tors.379.14379.14379.14379.14379.14
Lam.ModeFSDT [28]FSDT [29]TE4TE3TE2TE1aFSDTaEBBTa
0 deg1st deg flex.108.23109.23109.24109.24113.95113.95113.95119.44
2nd deg flex.357.13357.15357.19403.38403.37403.37475.41
1st deg tors.379.13379.14379.14379.14379.14
15 deg1st deg flex.94.8669.5769.8370.0581.7381.7681.76108.39
2nd deg flex.261.23262.15262.81309.33309.47309.47431.46
1st deg tors.412.23412.79412.79413.60412.60
30 deg1st deg flex.76.7145.8246.0546.6052.4152.3352.3374.39
2nd deg flex.180.30181.11183.10205.85205.16205.16296.12
1st deg tors.519.80519.98519.98520.63519.80
45 deg1st deg flex.59.2137.8537.9838.4839.5739.5439.5443.62
2nd deg flex.150.13150.55152.51156.95156.30156.30173.62
1st deg tors.638.00639.43639.44691.04663.24
60 deg1st deg flex.45.8436.4936.5936.7436.7536.7536.7536.89
2nd deg flex.144.11144.30145.10145.35145.14145.14146.84
1st deg tors.534.03540.09540.10675.61586.18
75 deg1st deg flex.40.7438.3338.4638.5838.9939.0039.0039.58
2nd deg flex.149.60149.75150.42153.16153.21153.21157.54
1st deg tors.415.73419.76419.71460.93429.83
90 deg1st deg flex.40.1039.8039.9540.4740.7040.7940.7941.03
2nd deg flex.154.05154.18156.07159.34159.71159.71163.33
1st deg tors.379.14379.14379.14379.14379.14
a

Poisson locking correction activated.

### Thin-Walled Composite Shaft.

In the following numerical illustrations, box beams with thin walls, whose dimensions are shown in Table 1, have been studied. The results are displayed for different cross section height-to-width ratios ($R=b/c$, Fig. 2) and lamination angles (θ), in terms of nondimensional frequency ($f¯=f/fn$), as a function of the speed parameter ($Ω¯=Ω/fn$). The value fn is the natural frequency of a cantilever graphite–epoxy box beam with R = 1 and θ = 90 obtained with a 2D finite element solution (fn = 52.649 Hz). Figure 3 shows how the first two frequency ratios vary with the rotational speed for θ equal to 0 deg and 90 deg, using various expansions. As expected, the theory order does not affect the trend of the branches, except that curves start at lower values when the displacement field is enriched, and these differences are more evident for the case in which θ is equal to 0 deg. This fact demonstrates that classical models (and ad-hoc theories!) cannot guarantee the same accuracy for all lamination angles, thus confirming that a simple method to conceive increasingly accurate theories could be a very useful tool. Moreover, in order to evaluate the effect due to the ply angle, the first two frequency ratios are shown in Fig. 4, adopting the TE2 and TE6 expansions. Critical speeds occur at different values, depending on the lamination angle, and maximum and minimum velocities occur for θ = 0 deg and θ = 90 deg, respectively. Similar results were obtained by Librescu et al. in Ref. [12]. In addition, it should be underlined that different lamination schemes determine the occurrence of couplings between different deformation modes. For instance, referring to Fig. 2, two different configurations $[θT/θL/θB/θR]$ were studied: case I$[45 deg/-45 deg/-45 deg/45 deg]$ and case II$[45 deg/-45 deg/45 deg/-45 deg]$. Free vibration analyses were performed to compare the first eight frequencies with those of a 2D finite element solution (Table 6). It should be noted that, for both cases, higher-order theories yield results that are increasingly closer to references. The frequencies (f) and damping (D) ratios were computed as functions of the rotating speed using a number of displacement models (Figs. 5 and 6). Classical theories (EBBT and FSDT, Fig. 5(a)) furnish a qualitatively similar result for the first configuration, but do not predict any instability range. Nevertheless, with the first-order shear deformation theory, a veering of the fourth, fifth, and sixth branches can be detected. Instead, with the TE3 and TE6 expansions, the graphs show dramatic changes, and three instability fields in fact appear within the considered speed interval. This first configuration involves a dominant twist motion and models that are able to describe bending–torsional coupling are therefore needed. In addition, when the expansion order is increased, instability thresholds occur at lower speed values. On the contrary, for the second case, the most important effect is due to shear and, for this reason, the shearable beam model leads to comparable results with those obtained with higher-order theories. For both schemes, the gyroscopic coupling markedly affects the system, causing several veerings of frequency branches.

Fig. 2
Fig. 2
Close modal
Fig. 3
Fig. 3
Close modal
Fig. 4
Fig. 4
Close modal
Fig. 5
Fig. 5
Close modal
Fig. 6
Fig. 6
Close modal
Table 6

Frequencies (Hz) at standstill for box beams R = 1

Theoryf1f2f3f4f5af6bf7f8
[45 deg/−45 deg/−45 deg/45 deg]
2D63.27563.593370.59371.97632.51776.19947.78956.00
TE664.93565.211381.51383.19651.53787.53985.63994.39
TE465.55165.992387.31387.47662.00802.231005.01008.8
TE366.63967.597393.81395.89705.92803.261025.51025.8
TE267.00769.094400.25410.84710.29818.391054.91083.3
FSDT73.36473.490429.38434.77811.121103.61126.6
EBBT74.02274.022455.36455.36897.341239.21239.2
[45 deg/−45 deg/45 deg/−45 deg]
2D60.49566.500358.52391.08628.60775.52927.71995.54
TE661.75668.114367.14402.02646.76785.23958.721033.6
TE462.57268.886372.17406.98652.74795.59974.161049.5
TE362.94469.137376.46411.67674.30797.64992.811072.9
TE268.08975.689408.50451.06723.26803.111079.51175.3
FSDT65.32171.032389.45418.86897.381017.21078.1
EBBT74.01874.021455.34455.35897.381239.11239.2
Theoryf1f2f3f4f5af6bf7f8
[45 deg/−45 deg/−45 deg/45 deg]
2D63.27563.593370.59371.97632.51776.19947.78956.00
TE664.93565.211381.51383.19651.53787.53985.63994.39
TE465.55165.992387.31387.47662.00802.231005.01008.8
TE366.63967.597393.81395.89705.92803.261025.51025.8
TE267.00769.094400.25410.84710.29818.391054.91083.3
FSDT73.36473.490429.38434.77811.121103.61126.6
EBBT74.02274.022455.36455.36897.341239.21239.2
[45 deg/−45 deg/45 deg/−45 deg]
2D60.49566.500358.52391.08628.60775.52927.71995.54
TE661.75668.114367.14402.02646.76785.23958.721033.6
TE462.57268.886372.17406.98652.74795.59974.161049.5
TE362.94469.137376.46411.67674.30797.64992.811072.9
TE268.08975.689408.50451.06723.26803.111079.51175.3
FSDT65.32171.032389.45418.86897.381017.21078.1
EBBT74.01874.021455.34455.35897.381239.11239.2
a

Torsional mode.

b

Axial mode.

In the last illustrative example, a rectangular box beam (R = 0.5) was studied with the purpose of considering another kind of anisotropy. Hereafter, the reference frequency fn is equal to 25.976 Hz. Figure 7(a) has the purpose of showing how the first two frequency ratios vary with the speed parameter when θ is either 0 deg or 15 deg. The analyses were performed adopting four different theories and, for both cases, the curves move toward lower and lower values when the displacement model is enriched. For θ = 0 deg (thin lines), the TE3 branches almost overlap those of the TE6 theory, while discrepancies between these models clearly increase when θ is equal to 15 deg (bold lines). Therefore, in order to ensure a good accuracy for the other lamination schemes, the TE6 expansion was used and the results are shown in Fig. 7(b). It is interesting to note that the curves follow similar paths regardless of the lamination. In fact, after a more or less evident veering of the backward frequency branch, two eigenvalues interact and instability occurs. Furthermore, when θ increases, the instability threshold appears at lower and lower speeds, while the range of the instability field becomes smaller and smaller until it reaches a minimum for θ = 90 deg. The effects of ply angles on the damping are shown in Fig. 7(c).

Fig. 7
Fig. 7
Close modal

## Conclusion

In the present work, Carrera's unified formulation has been used to study the dynamics of composite rotors. By invoking Hamilton's principle, the equations of motion have been derived and solved with the finite element method. Several kinds of rotors have been considered to assess the new theory and the related results have been compared with those found in literature or with 2D finite element solutions. Higher-order elements were tested on shafts with circular and rectangular cross sections. In light of the results, it is possible to draw the following conclusions:

• The use of finite elements based on higher beam models leads to accurate results regardless to the lamination scheme.

• When the lamination configuration involves bending–torsional coupling, higher theories become mandatory to predict both instability thresholds and critical speeds.

• When composite structures are considered, classical theories yield acceptable qualitative results for a limited range of problems.

Future work could be focused on analyzing other aspects of rotating structures, for instance, by introducing the centrifugal stiffening contribution. The structural models proposed in this paper seem able to describe the dynamics of spinning structures with a high deformability, such as thin-walled cylinders and composite thin disks.

### Appendix A

In the following section there are explicitly reported coefficients of the transformed material matrix $C˜$ as a function of stiffness material coefficients $C$ and the two angles θ, β:
$C˜33=C33cos4θ+2(C23+2C66)cos2θsin2θ+C22sin4θC˜32=cos2β(C23(cos4θ+sin4θ)+(C33+C22-4C66)cos2θsin2θ)+sin2β(C13cos2θ+C21sin2θ)C˜31=sin2β(C23(cos4θ+sin4θ)+(C33+C22-4C66)cos2θsin2θ)+cos2β(C13cos2θ+C21sin2θ)C˜36=cosβ((2C66+C23-C33)cos3θsinθ+(-2C66+C22-C23)sin3θcosθ)C˜22=cos4β(C22cos4θ+(4C66+2C23)cos2θsin2θ+C33sin4θ) +cos2βsin2β((4C44+2C21)cos2θ+(4C55+2C13)sin2θ)+C11sin4βC˜12=cos2βsin2β((4C66+2C23)cos2θsin2θ+C33sin4θ-4C55sin2θ+C22cos4θ-4C44cos2θ+C11) +(cos4β+sin4β)(C13sin2θ+C21cos2θ)C˜26=cos3β((2C66+C23-C33)sin3θcosθ+(-2C66+C22-C33)cos3θsinθ) +cosβsin2β((2C44-2C55+C21-C13)cosθsinθ)C˜11=C11cos4β+cos2βsin2β((4C44+2C21)cos2θ+(4C55+2C13)sin2θ) +sin4β(C22cos4θ+C33sin4θ+(4C66+2C23)cos2θsin2θ)C˜16=cos3β(C21-C13)cosθsinθ+cosβsin2β((2C66+C23-C33)cosθsin3θ +(-2C66+C22-C23)cos3θsinθ)C˜44=(cos4β+sin4β)(C44cos2θ+C55sin2θ)+cos2βsin2β(C11+C22cos4θ +(-2C44-2C12)cos2θ+(-2C55-2C13)sin2θ+(4C66+2C23)cos2θsin2θ)C˜45=cos3β(C44-C55)cosθsinθ+cosβsin2β((-2C66+C22-C23)cos3θsinθ +(2C66+C23-C33)cosθsin3θ+(C55-C44-C12+C13))C˜55=cos2β(C55cos2θ+C44sin2θ)+sin2β(C66(cos4θ+sin4θ) +(C33+C22-2C66-2C23)cos2θsin2θ)C˜66=sin2β(C55cos2θ+C44sin2θ)+cos2β(C66(cos4θ+sin4θ) +(C33+C22-2C66-2C23)cos2θsin2θ)C˜34=cosβsinβ((C33+C22-4C66)cos2θsin2θ+C23(cos4β+sin4β)-C13cos2θ-C12sin2θ)C˜35=sinβ((C22-C23-2C66)cosθsin3θ+(-C33+C23+2C66)cos3θsinθ)C˜24=cos3βsinβ(C22cos4θ-(C21+2C44)cos2θ-(4C66+2C23)cos2θsin2θ +(C13+2C55)sin2θ+C33sin4θ)+cosβsin3β((C13+2C55)sin2θ+(C12+2C44)cos2θ-C11)C˜25=sin3β(C21-C13)cosθsinθ+cosβsin3β((2C66+C23-C33)cosθsin3θ +(-2C66-C23+C22)cos3θsinθ+(2C55-2C44)cosθsinθ)C˜14=cosβsin3β(C22cos4θ-(C21+2C44)cos2θ-(4C66+2C23)cos2θsin2θ +(C13+2C55)sin2θ+C33sin4θ)+cos3βsinβ((C13+2C55)sin2θ+(C12+2C44)cos2θ-C11)C˜15=sin3β((2C66+C23-C33)cosθsin3θ+(-2C66-C23+C22)cos3θsinθ)+cos2βsinβ(2C44-2C55+C21-C13)cosθsinθC˜45=sin3β(C44-C55)cosθsinθ+cos2βsinβ((2C66+C23-C33)cosθsin3θ +(-2C66-C23+C22)cos3θsinθ+(C55-C44-C21+C13)cosθsinθ)C˜56=cosβsinβ(C66(cos4θ+sin4θ)+(C22+C33-2C23-2C66)cos2θsin2θ-C44sin2θ-C55cos2θ)$
(A1)

### Appendix B

For a cross section made of nonhomogeneous orthotropic material, the components of the fundamental nucleus $Kijτs$ are written here as
$Kxx=Ili, yj◃FτC˜46Fs,z▹+Ili, yj◃FτC˜26Fs,x▹+Ili, yj, y◃FτC˜66Fs▹ + Ilij◃Fτ,zC˜44Fs,z▹+Ilij◃Fτ,zC˜24Fs,x▹+Ilij, y◃Fτ,zC˜46Fs▹ + Ilij, y◃Fτ,xC˜26Fs▹+Ilij◃Fτ,xC˜24Fs,z▹+Ilij◃Fτ,xC˜22Fs,x▹$
$Kxy=Ili, yj◃FτC˜66Fs,x▹+Ili, yj◃FτC˜56Fs,z▹+Ili, yj, y◃FτC˜36Fs▹ + Ilij◃Fτ,xC˜26Fs,x▹+Ilij◃Fτ,xC˜25Fs,z▹+Ilij◃Fτ,zC˜46Fs,x▹ + Ilij◃Fτ,zC˜45Fs,z▹+Ilij, y◃Fτ,zC˜43Fs▹+Ilij, y◃Fτ,xC˜23Fs▹$
$Kxz=Ili, yj◃FτC˜46Fs,x▹+Ili, yj◃FτC˜16Fs,z▹+Ili, yj, y◃FτC˜56Fs▹ + Ilij◃Fτ,zC˜44Fs,x▹+Ilij◃Fτ,zC˜14Fs,z▹+Ilij◃Fτ,xC˜24Fs,x▹ + Ilij◃Fτ,xC˜21Fs,z▹+Ilij, y◃Fτ,zC˜45Fs▹+Ilij, y◃Fτ,xC˜25Fs▹$
$Kyx=Ilij, y◃Fτ,xC˜66Fs▹+Ilij, y◃Fτ,zC˜56Fs▹+Ili, yj◃FτC˜43Fs,z▹ + Ili, yj◃FτC˜23Fs,x▹+Ili, yj, y◃FτC˜36Fs▹+Ilij◃Fτ,xC˜46Fs,x▹ + Ilij◃Fτ,xC˜26Fs,x▹+Ilij◃Fτ,zC˜45Fs,z▹+Ilij◃Fτ,zC˜25Fs,x▹$
$Kyy=Ilij◃Fτ,xC˜66Fs,x▹+Ilij◃Fτ,xC˜56Fs,z▹+Ilij◃Fτ,zC˜56Fs,x▹ + Ilij◃Fτ,zC˜55Fs,z▹+Ilij, y◃Fτ,xC˜36Fs▹+Ilij, y◃Fτ,zC˜35Fs▹ + Ili, yj◃FτC˜36Fs,x▹+Ili, yj◃FτC˜35Fs,z▹+Ili, yj, y◃FτC˜33Fs▹$
$Kyz=Ilij◃Fτ,xC˜46Fs,x▹+Ilij◃Fτ,xC˜16Fs,z▹+Ilij◃Fτ,zC˜45Fs,x▹ + Ilij◃Fτ,zC˜15Fs,z▹+Ilij, y◃Fτ,xC˜56Fs▹+Ilij, y◃Fτ,zC˜55Fs▹ + Ili, yj◃FτC˜43Fs,x▹+Ili, yj◃FτC˜13Fs,z▹+Ili, yj, y◃FτC˜35Fs▹$
$Kzx=Ili, yj◃FτC˜45Fs,z▹+Ili, yj◃FτC˜25Fs,x▹+Ili, yj, y◃FτC˜56Fs▹ + Ilij◃Fτ,xC˜44Fs,z▹+Ilij◃Fτ,xC˜24Fs,x▹+Ilij◃Fτ,zC˜21Fs,x▹ + Ilij◃Fτ,zC˜14Fs,z▹+Ilij, y◃Fτ,xC˜46Fs▹+Ilij, y◃Fτ,zC˜16Fs▹$
$Kzy=Ili, yj◃FτC˜56Fs,x▹+Ili, yj◃FτC˜55Fs,z▹+Ili, yj, y◃FτC˜35Fs▹ + Ilij◃Fτ,zC˜16Fs,z▹+Ilij◃Fτ,zC˜15Fs,z▹+Ilij◃Fτ,xC˜46Fs,x▹ + Ilij◃Fτ,xC˜45Fs,z▹+Ilij, y◃Fτ,xC˜43Fs▹+Ilij, y◃Fτ,zC˜13Fs▹$
$Kzz=Ili, yj◃FτC˜45Fs,x▹+Ili, yj◃FτC˜15Fs,z▹+Ili, yj, y◃FτC˜55Fs▹ + Ilij◃Fτ,xC˜44Fs,x▹+Ilij◃Fτ,xC˜14Fs,z▹+Ilij◃Fτ,zC˜14Fs,x▹ + Ilij◃Fτ,zC˜11Fs,z▹+Ilij, y◃Fτ,xC˜45Fs▹+Ilij, y◃Fτ,zC˜15Fs▹$

## References

1.
Zinberg
,
H.
, and
Symmonds
,
M. F.
,
1970
, “
The Development of an Advanced Composite Tail Rotor Driveshaft
,”
26th Annual National Forum of the American Helicopter Society
,
Washington, DC, June 16–18
.
2.
Zorzi
,
E.
, and
Giordano
,
J. C.
,
1985
, “
Composite Shaft Rotor Dynamic Evaluation
,”
ASME Design Engineering Conference on Mechanical Vibrations and Noise
, Cincinnati, OH, September 10–13, ASME Paper No. 85-DET-114.
3.
Bauchau
,
O. A.
,
1981
, “
Design, Manufacturing and Testing of High Speed Rotating Graphite/Epoxy Shafts
,” Ph.D. thesis, Department of Aeronautics and Astronautics, MIT, Cambridge, MA.
4.
Bauchau
,
O. A.
,
1983
, “
Optimal Design of High Speed Rotating Graphite/Epoxy Shafts
,”
J. Compos. Mater.
,
17
(
2
), pp.
170
181
.10.1177/002199838301700205
5.
Bert
,
C. W.
, and
Kim
,
C. D.
,
1992
, “
The Effect of Bending–Twisting Coupling on the Critical Speed of a Drive Shaft
,”
6th Japan–U.S. Conference on Composite Materials
,
Orlando, FL, June 22–24
.
6.
Bert
,
C. W.
, and
Kim
,
C. D.
,
1995
, “
Whirling of Composite-Material Driveshafts Including Bending–Twisting Coupling and Transverse Shear Deformation
,”
ASME J. Vib. Acoust.
,
17
, pp.
17
21
.10.1115/1.2873861
7.
Chen
,
L. W.
, and
Peng
,
W. K.
,
1998
, “
The Stability Behavior of Rotating Composite Shafts Under Axial Compressive Loads
,”
Compos. Struct.
,
41
(3–4)
, pp.
253
263
.10.1016/S0263-8223(98)00022-1
8.
Chang
,
M. Y.
,
Chen
,
J. K.
, and
Chang
,
C. Y.
,
2004
, “
A Simple Spinning Laminated Composite Shaft Model
,”
Int. J. Solids Struct.
,
41
(3–4)
, pp.
637
662
.10.1016/j.ijsolstr.2003.09.043
9.
Singh
,
S. P.
, and
Gupta
,
K.
,
1996
, “
Compostite Shaft Rotordynamic Analysis Using a Layerwise Theory
,”
J. Sound Vib.
,
191
(
5
), pp.
739
756
.10.1006/jsvi.1996.0153
10.
Gubran
,
H. B. H.
, and
Gupta
,
K.
,
2005
, “
The Effect of Stacking Sequence and Coupling Mechanisms on the Natural Frequencies of Composite Shafts
,”
J. Sound Vib.
,
282
(1–2)
, pp.
231
248
.10.1016/j.jsv.2004.02.022
11.
Sino
,
R.
,
Baranger
,
T. N.
,
Chatelet
,
E.
, and
Jacquet
,
G.
,
2008
, “
Dynamic Analysis of a Rotating Composite Shaft
,”
Compos. Sci. Technol.
,
68
(
2
), pp.
337
345
.10.1016/j.compscitech.2007.06.019
12.
Song
,
O.
, and
Librescu
,
L.
,
1997
, “
Anisotropy and Structural Coupling on Vibration and Instability of Spinning Thin-Walled Beams
,”
J. Sound Vib.
,
204
(
3
), pp.
477
494
.10.1006/jsvi.1996.0947
13.
Song
,
O.
,
Librescu
,
L.
, and
Jeong
,
N.-H.
,
2000
, “
Vibration and Stability of Pretwisted Spinning Thin-Walled Composite Beams Featuring Bending–Bending Elastic Coupling
,”
J. Sound Vib.
,
237
(
3
), pp.
513
533
.10.1006/jsvi.2000.3100
14.
Na
,
S.
,
Yoon
,
H.
, and
Librescu
,
L.
,
2006
, “
Effect of Taper Ratio on Vibration and Stability of a Composite Thin-Walled Spinning Shaft
,”
Thin-Walled Struct.
,
44
(
3
), pp.
362
371
.10.1016/j.tws.2006.02.007
15.
Kim
,
C. D.
, and
Bert
,
C. W.
,
1993
, “
Critical Speed Analysis of Laminated Composite, Hollow Drive Shafts
,”
Compos. Eng.
,
3
(
7–8
), pp.
633
644
.10.1016/0961-9526(93)90087-Z
16.
Ramezani
,
S.
, and
,
M. T.
,
2009
, “
Free Vibration Analysis of Rotating Laminated Cylindrical Shells Under Different Boundary Conditions Using a Combination of the Layerwise Theory and Wave Propagation Approach
,”
Trans. B: Mech. Eng.
,
16
, pp.
168
176
.
17.
Zhao
,
X.
,
Liew
,
K. M.
, and
Ng
,
T. Y.
,
2002
, “
Vibrations of Rotating Cross-Ply Laminated Circular Cylindrical Shells With Stringer and Ring Stiffeners
,”
Int. J. Solids Struct.
,
39
(
2
), pp.
529
545
.10.1016/S0020-7683(01)00194-9
18.
Chatelet
,
E.
,
Lornage
,
D.
, and
Jacquet-Richardet
,
G.
,
2002
, “
A Three Dimensional Modeling of the Dynamic Behavior of Composite Rotors
,”
Int. J. Rotating Mach.
,
8
(
3
), pp.
185
192
.10.1155/S1023621X02000179
19.
Hua
,
L.
,
Khin-Yong
,
L.
, and
Then-Yong
,
N.
,
2005
,
Rotating Shell Dynamics
,
Elsevier
,
New York
.
20.
Carrera
,
E.
,
2002
, “
Theories and Finite Elements for Multilayered, Anisotropic, Composite Plates and Shells
,”
Arch. Comput. Methods Eng.
,
9
(
2
), pp.
87
140
.10.1007/BF02736649
21.
Carrera
,
E.
,
2003
, “
Theories and Finite Elements for Multilayered Plates and Shells: A Unified Compact Formulation With Numerical Assessment and Benchmarking
,”
Arch. Comput. Methods Eng.
,
10
(
3
), pp.
216
296
.10.1007/BF02736224
22.
Carrera
,
E.
,
2004
, “
Assessment of Theories for Free Vibration Analysis of Homogeneous and Multilayered Plates
,”
Shock Vib.
,
11
(3–4), pp.
261
270
.10.1155/2004/493584
23.
Carrera
,
E.
,
Giunta
,
G.
, and
Petrolo
,
M.
,
2011
,
Beam Structures. Classical and Advanced Theories
,
Wiley
,
New York
.
24.
Carrera
,
E.
,
Filippi
,
M.
, and
Zappino
E.
,
2013
, “
Laminated Beam Analysis by Polynomial, Trigonometric, Exponential and Zig-Zag Theories
,”
Eur. J. Mech. A/Solids
,
41
, pp.
58
69
.10.1016/j.euromechsol.2013.02.006
25.
Carrera
,
E.
,
Filippi
,
M.
, and
Zappino
,
E.
,
2013
, “
Free Vibration Analysis of Laminated Beam by Polynomial, Trigonometric, Exponential and Zig-Zag Theories
,”
J. Compos. Mater
. (in press).10.1177/0021998313497775
26.
Carrera
,
E.
,
Filippi
,
M.
, and
Zappino
,
E.
,
2013
, “
Analysis of Rotor Dynamic by One-Dimensional Variable Kinematic Theories
,”
ASME J. Gas Turbines Power
,
135
(
9
), p.
092501
.10.1115/1.4024381
27.
Carrera
,
E.
Filippi
,
M.
, and
Zappino
,
E.
,
2013
, “
Free Vibration Analysis of Rotating Composite Blades Via Carrera Unified Formulation
Compos. Struct.
,
106
, pp.
317
325
.10.1016/j.compstruct.2013.05.055
28.
Boukhalfa
,
A.
,
2011
, “Dynamic Analysis of a Spinning Laminated Composite-Material Shaft Using the HP-Version of the Finite Element Method,”