# Vibration and Instability of Rotating Composite Thin-Walled Shafts with Internal Damping

Authors: Ren Yongsheng; Zhang Xingqi; Liu Yanghang; Chen Xiulong
Journal: Shock and Vibration
 (2014)

---

## Abstract

The dynamical analysis of a rotating thin-walled composite shaft with internal damping is carried out analytically. The equations of motion are derived using the thin-walled composite beam theory and the principle of virtual work. The internal damping of shafts is introduced by adopting the multiscale damping analysis method. Galerkin’s method is used to discretize and solve the governing equations. Numerical study shows the effect of design parameters on the natural frequencies, critical rotating speeds, and instability thresholds of shafts.


---

## Body

## 1. Introduction

The rotating composite material shafts are being used as structural elements in many application areas involving the rotating machinery systems. This is likely to contribute to the high strength to weight ratio, lower vibration level, and longer service life of composite materials. A significant weight saving can be achieved by the use of composite materials. Also by appropriate design of the composite layup configuration, orientation, and number of plies the improved performance of the shaft system can be obtained. Furthermore, the use of composite would permit the use of longer shafts in the supercritical range than what is possible with conventional metallic shafts. In the last few years, there exist numerous researches related to predicting critical speeds and natural frequencies of composite shaft. Zinberg and Symonds [1] investigated the critical speeds of rotating anisotropic cylindrical shafts based on an equivalent modulus beam theory (EMBT), and dos Reis et al. [2] evaluated the shaft of Zinberg and Symonds [1] by the finite element method. Kim and Bert [3] adopted the thin- and thick-shell theories of first-order approximation to derive the motion equations of the rotating composite thin-walled shafts. They used this model to obtain a closed form solution for a simply supported drive shaft and to analyze the critical speeds of composite shafts. Singh and Gupta [4] developed two composite spinning shaft models employing EMBT and layerwise beam theory (LBT), respectively. It was shown that a discrepancy exists between the critical speeds obtained from both models for the unsymmetric laminated composite shaft. Chang et al. [5] presented a simple spinning composite shaft model based on a first-order shear deformable beam theory. The finite element method is used here to find the approximate solution of the system. The model was used to analyze the critical speeds, frequencies, mode shapes, and transient response of a particular composite shaft system. Gubran and Gupta [6] presented a modified EMBT model to account for the effects of a stacking sequence and different coupling mechanisms. Song et al. [7] used Rehfield’s thin-walled beam theory [8] that presented a composite thin-walled shaft model. The effects of rotatory inertias, axial edge load, and boundary conditions on the natural frequencies and stability of the system were investigated. Ren et al. [9] proposed another composite thin-walled shaft model by means of the composite thin-walled beam theory, an asymptotically correct theory referred to as variational asymptotically method (VAM). The flexible composite shaft is assumed to support on bearings which are modeled as springs and dampers and containing of the rigid disks mounted on it. The natural frequencies and critical rotating speeds of the rotating composite shaft with the variation of the lamination angle, ratios of length over radius, ratios of radius over thickness, and shear deformation are then analyzed.On the other hand, composite material shafts have higher internal damping than conventional metallic shafts. However, as it has been shown in previous research [10] that internal damping in rotating assemblies may lead to whirl instabilities in high speed rotors. Therefore, accurate prediction of effects of internal damping in composite material rotors is essential. So far, there have been less work related to the stability analysis, particularly effects of internal damping on the stability of a rotating composite shaft [4, 11–14]. Singh and Gupta [4, 11] introduced discrete viscous damping coefficients to account for effect of internal damping. In a similar approach, dynamic instability analysis of composite shaft has been performed by Montagnier and Hochard [13], Mazzei and Scott [12], Kim et al. [14]. In above cases, internal damping terms were included simply in equations of motion of rotating shaft, no internal damping modeling has been described. Sino et al. [15] investigated the stability of an internally damped rotating composite shaft. Internal damping was introduced by the complex constitutive relation of a viscoelastic composite. The shaft was modeled by finite element method. However, only the mechanical coupling effects induced by symmetrical stacking are taken into account in their mode.In the present work, an analytical model applicable to the dynamical analysis of rotating composite thin-walled shafts with internal damping is proposed. This model is based on the composite thin-walled beam theory, referred to as variational asymptotically method (VAM) by Berdichevsky et al. [16]. The internal damping of shaft is introduced via the multiscale damping mechanics [17]. The equations of motion of the composite shafts are derived by the principle of virtual work. Galerkin’s method is used then to discretize and solve the governing equations. The natural frequencies, critical rotating speeds, and instability thresholds are obtained through numerical simulations. The effect of the ply angle and aspect ratio of cross-section are then assessed. The validity of the model is proved by comparing the results with those in literature and convergence examination.


## 2. Model of Composite Shaft

Consider the slender thin-walled composite shaft given in Figure1. The length of the shaft is denoted by L, its thickness is denoted by h, the radius of curvature of the middle is denoted by r, and the maximum cross-sectional dimension is denoted by d. It is assumed that r
≪
L, h
≪
r, and the coordinate ζ is measured along the normal to the middle surface within the limits -
h
/
2
⩽
ζ
⩽
h
/
2. The shaft rotates about its longitudinal x-axis at a constant rate Ω.Figure 1
Geometry and coordinate systems of composite thin-walled shaft.


### 2.1. The Displacement and Strain

The components of the displacement along the Cartesian coordinate (x
,
y
,
z) are expressed as follows [16]:
(1)
u
1
x
,
s
=
u
x
-
y
s
v
′
x
-
z
s
w
′
x
+
g
s
,
x
,
u
2
x
,
s
=
v
x
-
z
s
φ
x
,
u
3
x
,
s
=
w
x
+
y
s
φ
x
,
where the primes in ((1)) denote differentiation with respect to x.The tangential and normal displacementsv
2 and v
ζ can be expressed as follows:
(2)
v
2
=
u
2
d
y
d
s
+
u
3
d
z
d
s
,
v
ζ
=
u
2
d
z
d
s
-
u
3
d
y
d
s
.Based on the expressions shown in ((1)) and ((2)), the in-plane strain components can be written in terms of the displacement variables as follows:
(3)
γ
11
=
∂
u
1
∂
x
,
2
γ
12
=
∂
u
1
∂
s
+
∂
v
2
∂
x
,
γ
22
=
∂
v
2
∂
s
+
v
ζ
r
.


### 2.2. Equations of Motion

The equations of motion of the shaft free vibration can be described by a variational form(4)
∫
0
l
δ
U
s
-
δ
T
s
+
δ
W
s
d
x
=
0
,
where δ
U
s, δ
T
s, and δ
W
s are the variation of the strain energy, the kinetic energy, and the dissipated energy of the cross-section, respectively.The variation of the strain energy of the cross-section can be expressed as follows:(5)
δ
U
s
=
∫
A
δ
ε
x
,
y
T
Q
-
i
j
ε
x
,
y
d
s
d
ζ
,
where ε
x
,
y
T
=
γ
11
γ
22
2
γ
12, A is the cross-sectional area of the shaft, and [
Q
-
i
j
] is equivalent off-axis stiffness matrix.The variation of the dissipated energy of the cross-section can be expressed as follows [17]:
(6)
δ
W
s
=
∫
A
δ
ε
x
,
y
T
η
-
i
j
Q
-
i
j
ε
x
,
y
d
s
d
ζ
,
where [
η
-
i
j
] is equivalent off-axis damping matrix.Further, the variation of the kinetic energy of the cross-section with rotating motion is(7)
δ
T
s
=
∫
A
δ
R
T
d
i
a
g
ρ
R
¨
d
s
d
ζ
,
where [
d
i
a
g
(
ρ
)
] is a diagonal matrix with components equal to the mass density ρ of a ply.The position, velocity, and acceleration vectors for the deformed shaft are described as follows:(8)
R
=
y
+
u
2
i
+
z
+
u
3
j
+
x
+
u
1
k
,
R
˙
=
u
˙
2
-
Ω
z
+
u
3
i
+
u
˙
3
+
Ω
y
+
u
2
j
+
u
˙
1
k
,
R
¨
=
u
¨
2
-
2
Ω
u
˙
3
-
Ω
2
y
+
u
2
i
+
u
¨
3
+
2
Ω
u
˙
2
-
Ω
2
z
+
u
3
j
+
u
¨
1
k
.


### 2.3. Equivalent Cross-Section Stiffness Matrix

For the case of no internal pressure acting on the shaft, ((5)) can be simplified by using free hoop stress resultant (N
22
=
0) assumption as follows:
(9)
δ
U
s
=
1
2
∮
δ
γ
11
δ
γ
12
A
B
B
C
γ
11
γ
12
d
s
,
where ∮
(
·
)
d
s denotes the integral around the loop of the midline cross-section, and the reduced axial, coupling, and shear stiffness A, B, and C can be written as follows:
(10)
A
s
=
A
11
-
A
12
2
A
22
,
B
s
=
2
A
16
-
A
12
A
26
A
22
,
C
s
=
4
A
66
-
A
26
2
A
22
,
A
i
j
=
∑
k
=
1
N
Q
-
i
j
k
z
k
-
z
k
-
1
,
i
,
j
=
1,2
,
6
.In ((9)), δ
U
s can also be expressed with respect to u, v, w, and φ by combining ((1))–((3)) and ((9)). Thus, one has
(11)
δ
U
s
=
δ
Δ
T
K
Δ
,
where Δ is 4 × 1 column matrix of kinematic variables defined as Δ
T
=
(
u
′
φ
′
w
′′
v
′′
) and [
K
] is 4 × 4 symmetric stiffness matrix. Its components k
i
j are given as follows:
(12)
k
11
=
∮
A
-
B
2
C
d
s
+
∮
B
/
C
d
s
2
∮
1
/
C
d
s
,
k
12
=
∮
B
/
C
d
s
∮
1
/
C
d
s
A
e
,
k
13
=
-
∮
A
-
B
2
C
z
d
s
-
∮
B
/
C
d
s
∮
B
/
C
z
d
s
∮
1
/
C
d
s
,
k
14
=
-
∮
A
-
B
2
C
y
d
s
-
∮
B
/
C
d
s
∮
B
/
C
y
d
s
∮
1
/
C
d
s
,
k
22
=
1
∮
1
/
C
d
s
A
e
2
,
k
23
=
-
∮
B
/
C
z
d
s
∮
1
/
C
d
s
A
e
,
k
24
=
-
∮
B
/
C
y
d
s
∮
1
/
C
d
s
A
e
,
k
33
=
∮
A
-
B
2
C
z
2
d
s
+
∮
B
/
C
z
d
s
2
∮
1
/
C
d
s
,
k
34
=
-
∮
A
-
B
2
C
y
z
d
s
-
∮
B
/
C
y
d
s
∮
B
/
C
z
d
s
∮
1
/
C
d
s
,
k
44
=
∮
A
-
B
2
C
y
2
d
s
+
∮
Γ
B
/
C
y
d
s
2
∮
1
/
C
d
s
,
A
e
=
1
2
∮
y
d
z
d
s
-
z
d
y
d
s
d
s
.


### 2.4. Equivalent Cross-Section Damping Matrix

Similar to the derivation of the previous cross-section stiffness formulations, the variation of the dissipated energy of the cross-section in terms of the strainsγ
11 and γ
12 can be modeled as follows:
(13)
δ
W
s
=
∮
γ
11
γ
12
A
d
B
d
C
d
D
d
γ
11
γ
12
d
s
,
where
(14)
A
d
=
A
d
11
+
A
12
2
A
22
2
A
d
22
-
A
12
A
22
A
d
12
+
A
d
21
,
B
d
=
A
d
16
+
A
d
61
+
2
A
12
A
26
A
22
2
A
d
22
-
A
26
A
22
A
d
12
+
A
d
21
-
A
12
A
22
A
d
26
+
A
d
62
,
C
d
=
4
A
d
66
+
A
26
2
A
22
2
A
d
22
-
A
26
A
22
A
d
26
+
A
d
62
,
A
d
i
j
=
∫
-
h
/
2
h
/
2
ψ
-
i
l
Q
-
l
j
d
ζ
=
2
∑
k
=
1
N
/
2
ψ
-
i
l
k
Q
-
l
j
k
h
k
-
h
k
-
1
,
22222222222222.2222222222
i
,
j
,
l
=
1,2
,
6
.The variation of the dissipated energy can be also expressed in terms of the kinematic variables as follows:(15)
δ
W
s
=
δ
Δ
T
C
Δ
,
where [
C
] is 4 × 4 symmetric damping matrix. The formulation of its components c
i
j is analogous to stiffness components k
i
j as shown in ((12)), but the terms A, B, and C in ((10)) should be replaced by the terms A
d, B
d, and C
d, respectively.


### 2.5. Equivalent Cross-Section Mass

Substituting ((1)) into ((8)) and in view of ((7)), the variation of the kinetic energy of the cross-section can be obtained as follows:
(16)
δ
T
s
=
-
I
1
δ
u
+
I
2
δ
v
+
I
3
δ
w
+
I
4
δ
φ
,
where
(17)
I
1
=
b
1
u
¨
,
I
2
=
b
1
v
¨
-
2
Ω
w
˙
-
Ω
2
v
-
b
2
2
Ω
φ
˙
+
Ω
2
-
b
3
φ
¨
-
Ω
2
φ
,
I
3
=
b
1
w
¨
+
2
Ω
v
˙
-
Ω
2
w
+
b
2
φ
¨
-
Ω
2
φ
-
b
3
2
Ω
φ
˙
+
Ω
2
,
I
4
=
b
2
w
¨
+
2
Ω
v
˙
-
Ω
2
w
-
b
3
v
¨
-
2
Ω
w
˙
-
Ω
2
v
+
b
4
+
b
5
φ
¨
-
Ω
2
φ
,
b
1
=
∫
A
ρ
d
s
d
ζ
,
b
2
=
∫
A
ρ
y
d
s
d
ζ
,
b
3
=
∫
A
ρ
z
d
s
d
ζ
,
b
4
=
∫
A
ρ
y
2
d
s
d
ζ
,
b
5
=
∫
A
ρ
z
2
d
s
d
ζ
.


### 2.6. Approximate Solution Method

In order to find the approximate solution of the rotating composite shaft, the quantitiesu
(
x
,
t
), v
(
x
,
t
), w
(
x
,
t
), and φ
(
x
,
t
) are assumed in the form
(18)
u
x
,
t
=
∑
j
=
1
N
A
j
α
j
x
e
i
λ
t
,
φ
x
,
t
=
∑
j
=
1
N
B
j
θ
j
x
e
i
λ
t
,
v
x
,
t
=
∑
j
=
1
N
C
j
ψ
j
x
e
i
λ
t
,
w
x
,
t
=
∑
j
=
1
N
D
j
ψ
j
x
e
i
λ
t
,
where α
j
(
x
), θ
j
(
x
), and ψ
j
(
x
)are mode shape functions which fulfill all the boundary conditions of the composite shaft; λ is complex eigenvalues of the system; A
j, B
j, C
j, and D
j are undetermined constants; and i
=
-
1.Substituting ((18)) into the governing equations of motion equations ((4))–((6)) and applying Galerkin’s procedure, the following governing equations in matrix form can be found:
(19)
δ
U
T
-
λ
2
M
+
i
λ
G
+
i
λ
C
+
K
U
=
0
,
where U
T
=
(
A
1
,
A
2
,
…
,
A
N
,
B
1
,
B
2
,
…
,
B
N
,
C
1
,
C
2
,
…
,
C
N
,
D
1
,
D
2
,
…
,
D
N
) is a constant vector, [
M
] is the mass matrix; [
G
] is the gyroscopic matrix, [
C
] is the damping matrix, and [
K
] is the stiffness matrix which also includes contribution from the centrifugal forces. The detailed expressions of these matrices are as follows:(20)
M
=
-
b
1
H
i
j
0
0
0
0
-
b
4
+
b
5
L
i
j
-
b
2
M
i
j
b
3
M
i
j
0
b
2
Q
i
j
b
1
R
i
j
0
0
-
b
3
Q
i
j
0
b
1
R
i
j
,
G
=
0
0
0
0
0
0
0
0
0
0
0
2
b
1
R
i
j
Ω
0
0
-
2
b
1
R
i
j
Ω
0
,
C
=
c
11
E
i
j
c
12
F
i
j
c
13
G
i
j
c
14
G
i
j
c
12
I
i
j
c
22
J
i
j
c
23
K
i
j
c
24
K
i
j
c
13
N
i
j
c
23
O
i
j
c
33
P
i
j
c
34
P
i
j
c
14
N
i
j
c
24
O
i
j
c
34
P
i
j
c
44
P
i
j
,
K
=
k
11
E
i
j
k
12
F
i
j
k
13
G
i
j
k
14
G
i
j
k
12
I
i
j
k
22
J
i
j
+
b
4
+
b
5
L
i
j
Ω
2
k
23
K
i
j
k
24
K
i
j
k
13
N
i
j
k
23
O
i
j
k
33
P
i
j
-
b
1
R
i
j
Ω
2
k
34
P
i
j
k
14
N
i
j
k
24
O
i
j
k
34
P
i
j
k
44
P
i
j
-
b
1
R
i
j
Ω
2
,where(21)
E
i
j
=
∫
0
L
α
i
α
j
′′
d
x
,
F
i
j
=
∫
0
L
α
i
θ
j
′′
d
x
,
G
i
j
=
∫
0
L
α
i
ψ
j
′′′
d
x
,
H
i
j
=
∫
0
L
α
i
α
j
d
x
,
I
i
j
=
∫
0
L
θ
i
α
j
′′
d
x
,
J
i
j
=
∫
0
L
θ
i
θ
j
′′
d
x
,
K
i
j
=
∫
0
L
θ
i
ψ
j
′′′
d
x
,
L
i
j
=
∫
0
L
θ
i
θ
j
d
x
,
M
i
j
=
∫
0
L
θ
i
ψ
j
d
x
,
N
i
j
=
∫
0
L
ψ
i
α
j
′′′
d
x
,
O
i
j
=
∫
0
L
ψ
i
θ
j
′′′
d
x
,
P
i
j
=
∫
0
L
ψ
i
ψ
j
′′′′
d
x
,
Q
i
j
=
∫
0
L
ψ
i
θ
j
d
x
,
R
i
j
=
∫
0
L
ψ
i
ψ
j
d
x
,
i
,
j
=
1
,
…
,
N
.From ((19)), one can obtain the following complex eigenvalue problem:
(22)
det
⁡
-
λ
2
M
+
i
λ
G
+
C
+
K
=
0
.Complex eigenvalueλ can be expressed in the form
(23)
λ
=
σ
+
i
ω
.The damping natural frequency or whirl frequency of the system is the imaginary partω, whereas its real part σ gives the decay or growth of the amplitude of vibration. A negative value of σ indicates a stable motion, whereas a positive value indicates an unstable motion, growing exponentially in time.


## 2.1. The Displacement and Strain

The components of the displacement along the Cartesian coordinate (x
,
y
,
z) are expressed as follows [16]:
(1)
u
1
x
,
s
=
u
x
-
y
s
v
′
x
-
z
s
w
′
x
+
g
s
,
x
,
u
2
x
,
s
=
v
x
-
z
s
φ
x
,
u
3
x
,
s
=
w
x
+
y
s
φ
x
,
where the primes in ((1)) denote differentiation with respect to x.The tangential and normal displacementsv
2 and v
ζ can be expressed as follows:
(2)
v
2
=
u
2
d
y
d
s
+
u
3
d
z
d
s
,
v
ζ
=
u
2
d
z
d
s
-
u
3
d
y
d
s
.Based on the expressions shown in ((1)) and ((2)), the in-plane strain components can be written in terms of the displacement variables as follows:
(3)
γ
11
=
∂
u
1
∂
x
,
2
γ
12
=
∂
u
1
∂
s
+
∂
v
2
∂
x
,
γ
22
=
∂
v
2
∂
s
+
v
ζ
r
.


## 2.2. Equations of Motion

The equations of motion of the shaft free vibration can be described by a variational form(4)
∫
0
l
δ
U
s
-
δ
T
s
+
δ
W
s
d
x
=
0
,
where δ
U
s, δ
T
s, and δ
W
s are the variation of the strain energy, the kinetic energy, and the dissipated energy of the cross-section, respectively.The variation of the strain energy of the cross-section can be expressed as follows:(5)
δ
U
s
=
∫
A
δ
ε
x
,
y
T
Q
-
i
j
ε
x
,
y
d
s
d
ζ
,
where ε
x
,
y
T
=
γ
11
γ
22
2
γ
12, A is the cross-sectional area of the shaft, and [
Q
-
i
j
] is equivalent off-axis stiffness matrix.The variation of the dissipated energy of the cross-section can be expressed as follows [17]:
(6)
δ
W
s
=
∫
A
δ
ε
x
,
y
T
η
-
i
j
Q
-
i
j
ε
x
,
y
d
s
d
ζ
,
where [
η
-
i
j
] is equivalent off-axis damping matrix.Further, the variation of the kinetic energy of the cross-section with rotating motion is(7)
δ
T
s
=
∫
A
δ
R
T
d
i
a
g
ρ
R
¨
d
s
d
ζ
,
where [
d
i
a
g
(
ρ
)
] is a diagonal matrix with components equal to the mass density ρ of a ply.The position, velocity, and acceleration vectors for the deformed shaft are described as follows:(8)
R
=
y
+
u
2
i
+
z
+
u
3
j
+
x
+
u
1
k
,
R
˙
=
u
˙
2
-
Ω
z
+
u
3
i
+
u
˙
3
+
Ω
y
+
u
2
j
+
u
˙
1
k
,
R
¨
=
u
¨
2
-
2
Ω
u
˙
3
-
Ω
2
y
+
u
2
i
+
u
¨
3
+
2
Ω
u
˙
2
-
Ω
2
z
+
u
3
j
+
u
¨
1
k
.


## 2.3. Equivalent Cross-Section Stiffness Matrix

For the case of no internal pressure acting on the shaft, ((5)) can be simplified by using free hoop stress resultant (N
22
=
0) assumption as follows:
(9)
δ
U
s
=
1
2
∮
δ
γ
11
δ
γ
12
A
B
B
C
γ
11
γ
12
d
s
,
where ∮
(
·
)
d
s denotes the integral around the loop of the midline cross-section, and the reduced axial, coupling, and shear stiffness A, B, and C can be written as follows:
(10)
A
s
=
A
11
-
A
12
2
A
22
,
B
s
=
2
A
16
-
A
12
A
26
A
22
,
C
s
=
4
A
66
-
A
26
2
A
22
,
A
i
j
=
∑
k
=
1
N
Q
-
i
j
k
z
k
-
z
k
-
1
,
i
,
j
=
1,2
,
6
.In ((9)), δ
U
s can also be expressed with respect to u, v, w, and φ by combining ((1))–((3)) and ((9)). Thus, one has
(11)
δ
U
s
=
δ
Δ
T
K
Δ
,
where Δ is 4 × 1 column matrix of kinematic variables defined as Δ
T
=
(
u
′
φ
′
w
′′
v
′′
) and [
K
] is 4 × 4 symmetric stiffness matrix. Its components k
i
j are given as follows:
(12)
k
11
=
∮
A
-
B
2
C
d
s
+
∮
B
/
C
d
s
2
∮
1
/
C
d
s
,
k
12
=
∮
B
/
C
d
s
∮
1
/
C
d
s
A
e
,
k
13
=
-
∮
A
-
B
2
C
z
d
s
-
∮
B
/
C
d
s
∮
B
/
C
z
d
s
∮
1
/
C
d
s
,
k
14
=
-
∮
A
-
B
2
C
y
d
s
-
∮
B
/
C
d
s
∮
B
/
C
y
d
s
∮
1
/
C
d
s
,
k
22
=
1
∮
1
/
C
d
s
A
e
2
,
k
23
=
-
∮
B
/
C
z
d
s
∮
1
/
C
d
s
A
e
,
k
24
=
-
∮
B
/
C
y
d
s
∮
1
/
C
d
s
A
e
,
k
33
=
∮
A
-
B
2
C
z
2
d
s
+
∮
B
/
C
z
d
s
2
∮
1
/
C
d
s
,
k
34
=
-
∮
A
-
B
2
C
y
z
d
s
-
∮
B
/
C
y
d
s
∮
B
/
C
z
d
s
∮
1
/
C
d
s
,
k
44
=
∮
A
-
B
2
C
y
2
d
s
+
∮
Γ
B
/
C
y
d
s
2
∮
1
/
C
d
s
,
A
e
=
1
2
∮
y
d
z
d
s
-
z
d
y
d
s
d
s
.


## 2.4. Equivalent Cross-Section Damping Matrix

Similar to the derivation of the previous cross-section stiffness formulations, the variation of the dissipated energy of the cross-section in terms of the strainsγ
11 and γ
12 can be modeled as follows:
(13)
δ
W
s
=
∮
γ
11
γ
12
A
d
B
d
C
d
D
d
γ
11
γ
12
d
s
,
where
(14)
A
d
=
A
d
11
+
A
12
2
A
22
2
A
d
22
-
A
12
A
22
A
d
12
+
A
d
21
,
B
d
=
A
d
16
+
A
d
61
+
2
A
12
A
26
A
22
2
A
d
22
-
A
26
A
22
A
d
12
+
A
d
21
-
A
12
A
22
A
d
26
+
A
d
62
,
C
d
=
4
A
d
66
+
A
26
2
A
22
2
A
d
22
-
A
26
A
22
A
d
26
+
A
d
62
,
A
d
i
j
=
∫
-
h
/
2
h
/
2
ψ
-
i
l
Q
-
l
j
d
ζ
=
2
∑
k
=
1
N
/
2
ψ
-
i
l
k
Q
-
l
j
k
h
k
-
h
k
-
1
,
22222222222222.2222222222
i
,
j
,
l
=
1,2
,
6
.The variation of the dissipated energy can be also expressed in terms of the kinematic variables as follows:(15)
δ
W
s
=
δ
Δ
T
C
Δ
,
where [
C
] is 4 × 4 symmetric damping matrix. The formulation of its components c
i
j is analogous to stiffness components k
i
j as shown in ((12)), but the terms A, B, and C in ((10)) should be replaced by the terms A
d, B
d, and C
d, respectively.


## 2.5. Equivalent Cross-Section Mass

Substituting ((1)) into ((8)) and in view of ((7)), the variation of the kinetic energy of the cross-section can be obtained as follows:
(16)
δ
T
s
=
-
I
1
δ
u
+
I
2
δ
v
+
I
3
δ
w
+
I
4
δ
φ
,
where
(17)
I
1
=
b
1
u
¨
,
I
2
=
b
1
v
¨
-
2
Ω
w
˙
-
Ω
2
v
-
b
2
2
Ω
φ
˙
+
Ω
2
-
b
3
φ
¨
-
Ω
2
φ
,
I
3
=
b
1
w
¨
+
2
Ω
v
˙
-
Ω
2
w
+
b
2
φ
¨
-
Ω
2
φ
-
b
3
2
Ω
φ
˙
+
Ω
2
,
I
4
=
b
2
w
¨
+
2
Ω
v
˙
-
Ω
2
w
-
b
3
v
¨
-
2
Ω
w
˙
-
Ω
2
v
+
b
4
+
b
5
φ
¨
-
Ω
2
φ
,
b
1
=
∫
A
ρ
d
s
d
ζ
,
b
2
=
∫
A
ρ
y
d
s
d
ζ
,
b
3
=
∫
A
ρ
z
d
s
d
ζ
,
b
4
=
∫
A
ρ
y
2
d
s
d
ζ
,
b
5
=
∫
A
ρ
z
2
d
s
d
ζ
.


## 2.6. Approximate Solution Method

In order to find the approximate solution of the rotating composite shaft, the quantitiesu
(
x
,
t
), v
(
x
,
t
), w
(
x
,
t
), and φ
(
x
,
t
) are assumed in the form
(18)
u
x
,
t
=
∑
j
=
1
N
A
j
α
j
x
e
i
λ
t
,
φ
x
,
t
=
∑
j
=
1
N
B
j
θ
j
x
e
i
λ
t
,
v
x
,
t
=
∑
j
=
1
N
C
j
ψ
j
x
e
i
λ
t
,
w
x
,
t
=
∑
j
=
1
N
D
j
ψ
j
x
e
i
λ
t
,
where α
j
(
x
), θ
j
(
x
), and ψ
j
(
x
)are mode shape functions which fulfill all the boundary conditions of the composite shaft; λ is complex eigenvalues of the system; A
j, B
j, C
j, and D
j are undetermined constants; and i
=
-
1.Substituting ((18)) into the governing equations of motion equations ((4))–((6)) and applying Galerkin’s procedure, the following governing equations in matrix form can be found:
(19)
δ
U
T
-
λ
2
M
+
i
λ
G
+
i
λ
C
+
K
U
=
0
,
where U
T
=
(
A
1
,
A
2
,
…
,
A
N
,
B
1
,
B
2
,
…
,
B
N
,
C
1
,
C
2
,
…
,
C
N
,
D
1
,
D
2
,
…
,
D
N
) is a constant vector, [
M
] is the mass matrix; [
G
] is the gyroscopic matrix, [
C
] is the damping matrix, and [
K
] is the stiffness matrix which also includes contribution from the centrifugal forces. The detailed expressions of these matrices are as follows:(20)
M
=
-
b
1
H
i
j
0
0
0
0
-
b
4
+
b
5
L
i
j
-
b
2
M
i
j
b
3
M
i
j
0
b
2
Q
i
j
b
1
R
i
j
0
0
-
b
3
Q
i
j
0
b
1
R
i
j
,
G
=
0
0
0
0
0
0
0
0
0
0
0
2
b
1
R
i
j
Ω
0
0
-
2
b
1
R
i
j
Ω
0
,
C
=
c
11
E
i
j
c
12
F
i
j
c
13
G
i
j
c
14
G
i
j
c
12
I
i
j
c
22
J
i
j
c
23
K
i
j
c
24
K
i
j
c
13
N
i
j
c
23
O
i
j
c
33
P
i
j
c
34
P
i
j
c
14
N
i
j
c
24
O
i
j
c
34
P
i
j
c
44
P
i
j
,
K
=
k
11
E
i
j
k
12
F
i
j
k
13
G
i
j
k
14
G
i
j
k
12
I
i
j
k
22
J
i
j
+
b
4
+
b
5
L
i
j
Ω
2
k
23
K
i
j
k
24
K
i
j
k
13
N
i
j
k
23
O
i
j
k
33
P
i
j
-
b
1
R
i
j
Ω
2
k
34
P
i
j
k
14
N
i
j
k
24
O
i
j
k
34
P
i
j
k
44
P
i
j
-
b
1
R
i
j
Ω
2
,where(21)
E
i
j
=
∫
0
L
α
i
α
j
′′
d
x
,
F
i
j
=
∫
0
L
α
i
θ
j
′′
d
x
,
G
i
j
=
∫
0
L
α
i
ψ
j
′′′
d
x
,
H
i
j
=
∫
0
L
α
i
α
j
d
x
,
I
i
j
=
∫
0
L
θ
i
α
j
′′
d
x
,
J
i
j
=
∫
0
L
θ
i
θ
j
′′
d
x
,
K
i
j
=
∫
0
L
θ
i
ψ
j
′′′
d
x
,
L
i
j
=
∫
0
L
θ
i
θ
j
d
x
,
M
i
j
=
∫
0
L
θ
i
ψ
j
d
x
,
N
i
j
=
∫
0
L
ψ
i
α
j
′′′
d
x
,
O
i
j
=
∫
0
L
ψ
i
θ
j
′′′
d
x
,
P
i
j
=
∫
0
L
ψ
i
ψ
j
′′′′
d
x
,
Q
i
j
=
∫
0
L
ψ
i
θ
j
d
x
,
R
i
j
=
∫
0
L
ψ
i
ψ
j
d
x
,
i
,
j
=
1
,
…
,
N
.From ((19)), one can obtain the following complex eigenvalue problem:
(22)
det
⁡
-
λ
2
M
+
i
λ
G
+
C
+
K
=
0
.Complex eigenvalueλ can be expressed in the form
(23)
λ
=
σ
+
i
ω
.The damping natural frequency or whirl frequency of the system is the imaginary partω, whereas its real part σ gives the decay or growth of the amplitude of vibration. A negative value of σ indicates a stable motion, whereas a positive value indicates an unstable motion, growing exponentially in time.


## 3. Numerical Results

The numerical calculations are performed by considering the shaft made of graphite-epoxy whose elastic characteristics are listed in Table1. The shaft has rectangular cross-section of fixed geometrical characteristics width a
=
0.32 m, length L
=
4.5952 m, and wall thickness h
=
0.01016 m, whose layup is [
θ
]
16 with clamped-free boundary conditions.Table 1
Mechanical properties of composite material [18].
ρ (kg/m3)
E
11 (GPa)
E
22 (GPa)
G
12 (GPa)
υ
12
η
l
1 (%)
η
l
2 (%)
η
l
6 (%)
1672
25.8
8.7
3.5
0.34
0.65
2.34
2.89In order to examine the influence of the number of mode shape functions used in the solution of the equation on the accuracy of the results, the numerical results of natural frequency are shown in Tables2 and 3 and modal damping in Tables 4 and 5 for an increasing number of mode shape functions, where L, a, and b are the length, width, and height, respectively. From these tables, it can be seen that to obtain accurate results of the first two natural frequencies and dampings, no more than five mode shape functions are required. This indicates clearly that the convergence of the present model is quite good.Table 2
Modal frequencies of cantilever composite box beam:L
/
a
=
14
.
36
,
a
/
b
=
5
,
[
0
]
16.
Mode
Natural frequency (Hz)
Present
Reference [18]
N
=
1
N
=
3
N
=
5
First flapping
3.12
3.13
3.13
3.1
Second flapping
—
19.02
19.02
19.8
First sweeping
10.80
10.81
10.81
11.0
Second sweeping
—
68.38
68.38
65.6
First torsional
37.83
37.84
37.84
37.7
Second torsional
112.91
113.01
113.01
113.3Table 3
Modal frequencies of cantilever composite box beam:L
/
a
=
14
.
36
,
a
/
b
=
5
,
[
90
]
16.
Mode
Natural frequency (Hz)
Present
Reference [18]
N
=
1
N
=
3
N
=
5
First flapping
1.80
1.81
1.81
1.8
Second flapping
—
11.36
11.36
11.5
First sweeping
5.67
6.20
6.20
6.5
Second sweeping
—
39.21
39.21
39.7
First torsional
37.83
37.84
37.84
37.7
Second torsional
112.90
113.01
113.01
113.3Table 4
Modal dampings of cantilever composite box beam:L
/
a
=
14
.
36
,
a
/
b
=
5
,
[
0
]
16.
Mode
Damping
Present
Reference [18]
N
=
1
N
=
3
N
=
5
First flapping
0.64
0.65
0.65
0.65
Second flapping
—
0.69
0.69
0.67
First sweeping
0.68
0.68
0.68
0.68
Second sweeping
—
0.85
0.85
0.9
First torsional
2.89
2.89
2.89
2.89
Second torsional
2.33
2.92
2.92
2.89Table 5
Modal dampings of cantilever composite box beam:L
/
a
=
14
.
36
,
a
/
b
=
5
,
[
90
]
16.
Mode
Damping
Present
Reference [18]
N
=
1
N
=
3
N
=
5
First flapping
2.39
2.38
2.38
2.35
Second flapping
—
2.37
2.37
2.35
First sweeping
3.03
2.53
2.53
2.35
Second sweeping
—
2.36
2.36
2.37
First torsional
2.89
2.89
2.89
2.89
Second torsional
2.82
2.92
2.92
2.89A comparison of predictions using the present model with those obtained in [18] is also shown in Tables 2, 3, 4, and 5. A perfect agreement of numerical results with those in [18] can be seen.Figure2 shows the variation of the first two flexural natural frequencies versus rotating speed for various ply angles. As it can be seen, because of the nonsymmetry of the shaft cross-section (a
/
b
≠
1), the standstill flexural frequencies about the two principal axes (flapping and sweeping denote bending about the y- and z-axis, resp.) are unequal. The behaviors of the flapping and sweeping bending frequencies versus rotating speed are very different. In fact, due to the existence of the Coriolis effect the coupling between flapping and sweeping bending is induced, the first decreases until it becomes zero, while the second continues to increase. It is observed that instead of a rotating speed, there is a whole domain of rotating speed in which the first flapping frequency does not exist. In this domain, the flapping frequency becomes imaginary value, implying that the shaft becomes unstable. When the ply angle is decreased, in addition to shift of instability domain towards larger rotating speeds, it is also observed that the domain of instability is enlarged.Figure 2
The variation of natural frequencies with the rotating speed for various ply angles (a
/
b
=
3.6, first two flexural modes).Figure3 shows the variation of the first two flexural dampings versus rotating speed for various ply angles. It can be seen clearly that as the rotating speed is increased, the damping of flapping bending mode decreases and remains negative for all rotating speed, so the flapping bending mode is stable. From the results of Figure 3, it can be also observed that the dampings corresponding to sweeping bending mode are negative at low rotating speed and increase with increasing rotating speed, and at certain value of rotating speed the dampings vanish and then become positive. Transformation of damping from a negative to a positive value marks the onset of unstable motion. The rotating speed corresponding to zero damping is the threshold of instability of the shaft. The enclosed curves located nearby the threshold of instability represent that the real parts are conjugate. Figure 3 also shows that the enclosed curves are shifted toward larger rotating speed but the extent of the enclosed curve is increased with the decrease of the ply angle.Figure 3
The variation of dampings with the rotating speed for various ply angles (a
/
b
=
3.6, first two flexural modes).Figure4 shows the variation of the first two extension-twist natural frequencies versus rotating speed for various ply angles. As seen in Figure 4, due to the absence of the Coriolis effect, the first extension-twist natural frequency (the twist is dominant) decreases, while the second (the extension is dominant) remains constant at all rotating speeds. From Figure 4 it is seen that the effect of ply angle on the first extension-twist natural frequency is significant and is quite different from the case of flexural mode. The maximum critical speed can be reached when the ply angle θ
=
45°.Figure 4
The variation of natural frequencies with the rotating speed for various ply angles (a
/
b
=
3.6, first two extension-twist modes).Figure5 shows the variation of the first two extension-twist dampings versus rotating speed for various ply angles. It can be observed that the second extension-twist mode is stable whereas the first extension-twist mode is unstable as the rotating speed increases above certain value. The self-excited range is easily identified from the figure by the sign of damping. It may also be noted that the effect of ply angle on the threshold of instability is similar to that previously described for the flexural mode.Figure 5
The variation of dampings with the rotating speed for various ply angles (a
/
b
=
3.6, first two extension-twist modes).Figure6 shows the variation of the first two flexural natural frequencies versus rotating speed for various aspect ratios. The results show that the critical rotating speed increases with the decrease of aspect ratios.Figure 6
The variation of natural frequencies with the rotating speed for various aspect ratios (θ
=
45° first two flexural modes).Figure7 shows the variation of the first two flexural dampings versus rotating speed for various aspect ratios. From the results it can be seen that the threshold of instability increases as aspect ratio decreases.Figure 7
The variation of dampings with the rotating speed for various aspect ratios (θ
=
45° first two flexural modes).Figures8 and 9 present the effect of aspect ratio on the natural frequency-rotating speed curves and damping-rotating speed curves for the extension-twist mode, respectively. The results show that the decrease of aspect ratio yields a significant increase of the critical rotating speed and the threshold of instability.Figure 8
The variation of natural frequencies with the rotating speed for various aspect ratios (θ
=
45°, first two extension-twist modes).Figure 9
The variation of dampings with the rotating speed for various aspect ratios (θ
=
45°, first two extension-twist modes).Figure10 shows the effect of ply angle on the critical rotating speed for the flexural mode. It can be seen that as the ply angle increases, the critical rotating speeds decrease and the maximum critical speed is maximum at θ
=
0°.Figure 10
The variation of critical speeds with ply angle for various aspect ratios (flexural mode).Figure11 shows the effect of ply angle on the threshold of instability for the flexural mode. It is evident that the general effect of the ply angle and aspect ratio on the threshold of instability is similar to that associated with the critical rotating speeds. By comparing Figure 10 with Figure 11, it may be noted that the threshold of instability is larger than the critical rotating speed and the difference between them increases as aspect ratio decreases. This implies that the onset of instability always occurs after the critical rotating speed.Figure 11
The variation of thresholds of instability with ply angle for various aspect ratios (flexural mode).Figures12 and 13 show the variation of the critical rotating speed and threshold of instability for the extension-twist mode, respectively. From these figures it becomes apparent that the maximum ones occur at θ
=
45°.Figure 12
The variation of critical speeds with ply angle for various aspect ratios (extension-twist mode).Figure 13
The variation of thresholds of instability with ply angle for various aspect ratios (extension-twist mode).


## 4. Conclusion

A model was presented for the study of the dynamical behavior of rotating thin-walled composite shaft with internal damping. The presented model was used to predict the natural frequencies, critical rotating speeds, and instability thresholds. Theoretical solutions of the free vibration of the shaft were determined by applying Galerkin’s method. From the present analysis and the numerical results, the following main conclusions were drawn.(1)
The developed model provides means of predicting the natural frequencies, critical rotating speeds, and instability thresholds of rotating composite thin-walled shafts with internal damping.
(2)
The ply angle and aspect ratio affect the vibrational and instability behavior of shaft significantly.
(3)
There is an obvious increase in the critical rotating speeds and instability thresholds as aspect ratio is decreased.
(4)
For the flexural mode, critical rotating speed and threshold of instability have their maximum values atθ
=
0°, while for the extension-twist mode, the maximum ones occur at θ
=
45°.
(5)
The onset of instability always occurs after the critical rotating speed.



---

*Source: 123271-2014-08-17.xml*