Modeling and dynamic analysis of rotating composite shaft

Ren Yongsheng1 , Dai Qiyi2 , Zhang Xingqi3

1, 2, 3College of Mechanical and Electronic Engineering, Shandong University of Science &Technology, Qindao, 266590, China

Journal of Vibroengineering, Vol. 15, Issue 4, 2013, p. 1790-1806.
Received 9 June 2013; accepted 5 November 2013; published 31 December 2013

Copyright © 2013 Vibroengineering This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Creative Commons License
Abstract.

Structural modeling and dynamical analysis of rotating composite shaft are conducted in this paper. A thin-walled composite shaft structure model, which includes the transverse shear deformation of the shaft, rigid disks and the flexible bearings, is presented and then used to predict natural frequencies and dynamical stability. Based on the thin-walled composite beam theory referred to as variational asymptotically method (VAM), the displacement and strain fields of the shaft are described. Hamilton’s principle is employed to derive the equations of motion of the shaft system. Galerkin’s method is used to discretize and solve the governing equations. The validity of the model is proved by comparing the results with those in literatures and convergence examination. The effects of fiber orientation, ratios of length over radius, ratios of radius over thickness and shear deformation on natural frequency and critical speeds are investigated. Finally the unbalance transient responses of the composite shaft system are also given by using the time-integration method.

Keywords: composite shaft, rotating shaft, vibration and stability.

1. Introduction

The rotating shafts made from laminated composite materials 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 a 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 for a specified critical speed than is possible with conventional metallic shafts.

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, the axial edge load and boundary conditions on the natural frequencies and stability of the system were investigated.

In the present study another composite thin-walled shaft model is proposed by means of the composite thin-walled beam theory, an asymptotically correct theory referred to as variational asymptotically method (VAM) by Berdichevsky et al. [9]. Here, however, the original formulation of VAM is expanded and refined to take into account effects of transverse shear deformation. The flexible composite shaft is assumed supporting on bearings which are modeled as springs and dampers and containing of the rigid disks mounted on it. The equations of motion of the composite shaft and the rotor-composite shaft system are derived by the extended Hamilton’s principle. Galerkin’s method is used then to discretize and solve the governing equations. The natural frequencies and critical rotating speeds of the rotating composite shaft with the variation of the lamination angle, the ratios of length over radius, ratios of radius over thickness and shear deformation are then analyzed. The validity of the model is proved by comparing the results with those in literatures and convergence examination. Finally the unbalance transient responses of a rotating composite shaft system are calculated by using the time-integration method.

2. Model and equations

2.1. Strain energy and kinetic energy of composite shaft

In the present study the composite shaft rotating along its longitudinal x-axis at a constant rate Ω is shown in Fig. 1. To describe the motion of the shaft the following coordinate systems are considered: (1) X,Y,Z is inertial reference system whose origin O is located in the geometric center, (2) x,y,z is rotating reference system with the common origin O.I,J,K and i,j,k denote the unit vectors of the reference systems X,Y,Z and x,y,z, respectively. In addition a local coordinate system ξ,s,x is used, where ξ and s are measured along the directions normal and tangent to the middle surface respectively.

The structural modeling is based on the following assumptions: (1) the shaft is characterized as a slender thin-walled elastic cylinder, satisfying dL, hd, hr, where L, h, r and d denote the length, the thickness, the radius of curvature and the maximum cross sectional dimension of the cylinder respectively, (2) transverse shear effects are considered, (3) the loop stress resultant may be ignored, (4) a special laminated composite configuration, referred to as the circumferentially uniform stiffness configuration (CUS) achieved by skewing angle plies with respect to the longitudinal beam axis meeting the conditions θy=θ-y, θz=θ-z, is considered (Fig. 1).

Fig. 1. Composite thin-walled shaft of a circular cross section

Composite thin-walled shaft of a circular cross section

It has been shown [10-12] that VAM is an asymptotically correct theory which can be used effectively for the analysis of tubular composite thin-walled structures. However, the classical VAM does not account for the effects of transverse shear because of which it may generate inaccurate predictions of the rotating composite shaft. Herein the original formulation of VAM is expanded and refined to include effects of transverse shear of the composite shaft.

The displacement field incorporating shear deformation is assumed in the form:

(1)
u 1 x , y , z , t = U 1 x , t - y s θ y x , t - z s θ z x , t + g s , x , t ,
u 2 x , y , z , t = U 2 x , t - z ϕ x , t ,
u 3 x , y , z , t = U 3 x , t + y ϕ x , t ,

in which:

(2)
θ y x , t = U 2 ' x , t - 2 γ z x ,
θ z x , t = U 3 ' x , t - 2 γ y x ,

where U1x,t,U2x,t,U3x,t denote the rigid-body translations along the x-, y- and z-axis, while ϕx,t,θyx,t,θzx,t denote the twist about the z-axis and rotations about the x- and y-axis respectively. γzx and γyx denote the transverse shear strains in the planes xz and xy respectively.

From the classical VAM the warping function gs,x,t is modified as:

(3)
g s , x , t = G s ϕ ' x , t + g 1 s U 1 ' x , t + g 2 s θ y ' x , t + g 3 s θ z ' x , t .

In the above equation the functions g1s,g2s,g3s,Gs are associated with physical behavior for the axial strain, the bending curvatures, and the torsion twist rate respectively. The primes in Eqs. (1) and (2) denote differentiation with respect to x.

Based on the displacement representations (1), (2) and (3), and using the linear strain-displacement relations [9], while referring to the [13], the strain of the composite shaft can be written as:

(4)
γ x x = U 1 ' x , t - y θ y ' x , t - z θ z   ' x , t ,
2 γ x s = d g d s + r n ϕ ' + U 2 ' x , t - θ y x , t d y d s + U 3 ' x , t - θ z x , t d z d s ,
2 γ x ξ = U 2 ' x , t - θ y x , t d z d s - U 3 ' x , t - θ z x , t d y d s   ,

where rn denotes the project of the position vector r of an arbitrary point on the cross-section of the deformed shaft in the normal direction:

(5)
r n = y d z d s - z d y d s .

The position vector r can be expressed as:

(6)
r = y + u 2 i + z + u 3 j + x + u 1 k .

According to the constant rotating rate assumption and using the expressions for the time derivatives of unit vectors, one can obtain the velocity of a generic point as:

(7)
V = r ˙ = u ˙ 2 - z + u 3 i + u ˙ 3 + Ω y + u 2 j + u ˙ 1 k .

The strain energy of the composite shaft Us can be written as:

(8)
U s = 1 2 0 L A   ( σ x x ε x x + σ x s ε x s + σ x ξ ε x ξ ) d A d x ,

where σxx, σxs and σxξ represent the engineering stresses associated with the engineering strains εxx, εxs and εxξ.

Taking variation for the above strain energy expression one obtains:

(9)
δ U s = 0 L A   ( σ x x δ ε x x + σ x s δ ε x s + σ x ξ δ ε x ξ ) d A d x = 0 L F x ' δ U 1 + Q y ' δ U 2 + Q z ' δ U 3 + M x ' δ ϕ + M z ' + Q y δ θ y + ( M y ' + Q z ) δ θ z d x .

Stress resultants Fx, Qy, Qz and stress couples Mx, My, Mz can be defined as:

(10)
F x     = N x x d s ,
M x = N x s r n d s ,     M y = - N x x z d s ,     M z = - N x x y d s ,
Q y = N x s d y d s + N x ξ d z d s d s ,     Q Z = N x s d z d s - N x ξ d y d s d s ,

where Nxx, Nxs and Nxξ are shell stress resultants defined according to the following expressions:

(11)
N x x N x s = A ( s ) B ( s ) 2 B ( s ) 2 C ( s ) 4 γ x x 2 γ x s ,
N x ξ     = D s γ x ξ ,

where:

(12)
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 ,     D s = A 44 - A 45 2 A 55 ,
A i j = k = 1 N Q ¯ i j k Z k - Z k - 1   , i , j = 1,2 , 6 ; i , j = 4,5 .

Parameters As and Bs denote the reduced axial and coupling stiffnesses respectively, while Cs and Ds denote the reduced shear stiffnesses. Aij denote the in-plane stiffnesses components and Q-ijk denote transformed stiffnesses.

Moreover by taking into account of Eqs. (4) and (11), the stress resultants and stress couples given by Eqs. (10) can be expressed in the following form:

(13)
F x = k 11 U 1 ' + k 12 ϕ ' + k 13 θ z ' + k 14 θ y ' + k 15 U 2 ' - θ y + k 16 U 3 ' - θ z ,
M x = k 12 U 1 ' + k 22 ϕ ' + k 23 θ z ' + k 24 θ y ' + k 25 U 2 ' - θ y + k 26 U 3 ' - θ z ,
M y = k 13 U 1 ' + k 23 ϕ ' + k 33 θ z ' + k 34 θ y ' + k 35 U 2 ' - θ y + k 36 U 3 ' - θ z ,
M z = k 14 U 1 ' + k 24 ϕ ' + k 34 θ z ' + k 44 θ y ' + k 45 U 2 ' - θ y + k 46 U 3 ' - θ z ,
Q y = k 51 U 1 ' + k 52 ϕ ' + k 53 θ z ' + k 54 θ y ' + k 55 U 2 ' - θ y + k 56 U 3 ' - θ z ,
Q z = k 61 U 1 ' + k 62 ϕ ' + k 63 θ z ' + k 64 θ y ' + k 65 U 2 ' - θ y + k 66 U 3 ' - θ z ,

where kij (i,j=1,,6) are the stiffness coefficients of the composite shaft, which can be expressed in terms of the cross sectional geometry and material properties as:

(14)
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 ,
k 15 = 1 2 Γ   B d y d s d s ,     k 16 = 1 2 Γ   B d z d s d s ,     k 25 = 1 4 Γ   r n C d y d s d s ,
k 26 = 1 4 Γ   r n C d z d s d s ,     k 35 = - 1 2 Γ   B z d y d s d s ,     k 36 = - 1 2 Γ   B z d z d s d s ,
k 45 = - 1 2 Γ   B y d y d s d s ,     k 46 = - 1 2 Γ   B y d z d s d s ,  
k 55 = Γ   1 4 C d y d s 2 + D d z d s 2 d s ,     k 56 = Γ   1 4 C - D d y d s d z d s d s ,
k 66 = Γ   1 4 C d z d s 2 + D d y d s 2 d s .

The coefficients kij (i,j=1,2,,4) represent the stiffnesses, as described in [9, 12] for the case of a non-shearable thin-walled beam, while the shear contribution to the stiffness coefficients is represented by the new terms kij  i=1,2,,6;j=5,6, kji (i=1,2,,4;j=5,6).

The kinetic energy of the composite shaft Ts can be written as:

(15)
T s = 1 2 0 L A   ρ ( V · V ) d A d x ,

where ρ denotes the mass density of the shaft. In view of Eq. (7) the expression of the kinetic energy can be obtained and taking variation of Ts yields:

(16)
δ T s = - 0 L I 1 δ U 1 + I 2 δ U 2 + I 3 δ U 3 + I 4 δ θ y + I 5 δ θ z + I 6 δ ϕ ,

where:

(17)
I 1 = b 1 U ¨ 1 - b 2 θ ¨ y - b 3 θ ¨ x ,
I 2 = b 1 U ¨ 2 - 2 Ω U ˙ 3 - Ω 2 U 2 - b 2 2 Ω ϕ ˙ + Ω 2 - b 3 ( ϕ ¨ - Ω 2 ϕ ) ,
I 3 = b 1 U ¨ 3 + 2 Ω U ˙ 2 - Ω 2 U 3 + b 2 ϕ ¨ - Ω 2 ϕ - b 3 ( 2 Ω ϕ ˙ + Ω 2 ) ,
I 4 = b 2 U ¨ 1 - b 4 θ ¨ y - b 6 θ ¨ z   ,
I 5 = b 3 U ¨ 3 - b 6 θ ¨ y - b 5 θ ¨ z   ,
I 6 = b 2 U ¨ 3 + 2 Ω U ˙ 2 - Ω 2 U 3 + b 4 + b 5 ϕ ¨ - Ω 2 ϕ - b 3 ( U ¨ 2 - 2 Ω U ˙ 3 - Ω 2 U 2 ) ,
b 1 = A   ρ d A ,     b 2 = A   ρ y d A ,     b 3 = A   ρ z d A ,
b 4 = A   ρ y 2 d A ,     b 5 = A   ρ z 2 d A ,     b 6 = A   ρ y z d A .

2.2. Kinetic energy of rigid disks

According to Eqs. (16) and (17) the expression of variation of the kinetic energy of the rigid disks fixed to the shaft is written as:

(18)
δ T D = - 0 L l = 1 N D I 1 l D δ U 1 + I 2 l D δ U 2 + I 3 l D δ U 3 + I 4 l D δ θ y + I 5 l D δ θ z + I 6 l D δ ϕ Δ x - x D l d x ,

in which:

(19)
I 1 l D = b 1 l D U ¨ 1 - b 2 l D θ ¨ y - b 3 l D θ ¨ Z ,
I 2 l D = b 1 l D ( U ¨ 2 - 2 Ω U ˙ 3 - Ω 2 U 2 ) - b 2 l D 2 Ω ϕ ˙ + Ω 2 - b 3 l D ( ϕ ¨ - Ω 2 ϕ ) ,
I 3 l D = b 1 l D U ¨ 3 + 2 Ω U ˙ 2 - Ω 2 U 3 + b 2 l D ϕ ¨ - Ω 2 ϕ - b 3 l D ( 2 Ω ϕ ˙ + Ω 2 ) ,
I 4 l D = b 2 l D U ¨ 1 - b 4 l D θ ¨ y - b 6 l D θ ¨ Z ,
I 5 l D = b 3 l D U ¨ 3 - b 6 l D θ ¨ y - b 5 l D θ ¨ Z ,
I 6 l D = b 2 l D U ¨ 3 + 2 Ω U ˙ 2 - Ω 2 U 3 + b 4 l D + b 5 l D ϕ ¨ - Ω 2 ϕ - b 3 l D U ¨ 2 - 2 Ω U ˙ 3 - Ω 2 U 2 ,
b 1 l D = A   ρ l D d A ,     b 2 l D = A   ρ l D   y d A ,     b 3 l D = A   ρ l D z d A ,
b 4 l D = A   ρ l D y 2 d A ,     b 5 l D = A   ρ l D z 2 d A ,     b 6 l D = A   ρ l D y z d A ,

where ρlD denotes the mass density of the disk l. The symbol Δ(x-xDl) denotes Dirac delta function, ND the number of the disks mounted on the shaft and xDl the location of the lth disk.

2.3. Work of external forces

The external forces include the centrifugal force resulting from the unbalanced rigid disks and the force from the bearings.

The bearings are considered as the springs and viscous dampers and these elements have linear stiffness and damping characteristics.

The virtual work associated with the springs and viscous dampers can be expressed as:

(20)
δ W b = 0 L l = 1 n b - K y y b l U 2 δ U 2 - K z y b l U 2 δ U 3 - K y z b l U 3 δ U 2 - K z z b l U 3 δ U 3 - C y y b l U ˙ 2 δ U ˙ 2 - C z y b l U ˙ 2 δ U ˙ 3   -   C y z b l U ˙ 3 δ U ˙ 2 - C z z b l U ˙ 3 δ U ˙ 3 Δ ( x - x b l ) d x ,

where nb denotes the number of the bearings, xbl the location, Kbl and Cbl the equivalent stiffness and damping coefficients of the lth bearings.

The virtual work done by the centrifugal force of rigid disks can be given by:

δ W e = 0 L p x u δ U 1 + p y u δ U 2 + p z u δ U 3 + Γ x u δ ϕ + Γ y u δ θ y + Γ z u δ θ z d x ,

where pxu, pyu, pzu denote centrifugal force per unit length, Γxu, Γyu, Γzu denote centrifugal moment per unit length. They are given by:

(21)
p x u = 0 ,       p y u = Ω 2 c o s Ω t l = 1 N D b 1 l D e l Δ ( x - x D l ) ,       p z u = Ω 2 s i n Ω t l = 1 N D b 1 l D e l Δ x - x D l ,
Γ x u = 0 ,     Γ y u = 0 ,     Γ z u = 0 ,

where the mass unbalance associated with disks is assumed as the concentrated masses b1lD=A ρlDdA at points with small distances of eccentricity el .

2.4. Motion equations of systems

In view of Eqs. (9), (16), (28) and (30), and employing the following variational principle:

(22)
t 0 t 1 δ U s - δ T s + T D - δ W b - δ W e d t = 0 ,

one obtains:

(23)
- F x ' + I 1 + I 1 D = p x u ,
- Q y ' + I 2 + I 2 D + F U 2 b = p y u ,
- Q z ' + I 3 + I 3 D + F U 3 b = p z u ,
- M x ' + I 6 + I 6 D = Γ x u ,
- M z ' + I 4 + I 4 D = Γ y u ,
- M y ' + I 5 + I 5 D = Γ z u ,

where:

(24)
F U 2 b = l = 1 n b K y y b l U 2 + K y z b l U 3 + K y y b l U ˙ 2 + K y z b l U ˙ 3 Δ x - x b l ,
F U 3 b = l = 1 n b K z z b l U 3 + K z y b l U 2 + K z z b l U ˙ 3 + K z y b l U ˙ 2 Δ ( x - x b l ) .

The boundary conditions at x=0, L are:

(25)
F x = F - x   o r   U 1 = U - 1 ,     Q y = Q - y   o r   U 2 = U - 2 ,     Q z = Q - z   o r   U 3 = U - 3 ,
M x = M - x   o r   θ y = θ - y ,     M y = M - y   o r   θ z = θ - z ,     M z = M - z   o r   ϕ = ϕ - .

As a result of CUS configuration, the equations of motion involving variables in terms of displacements can be reduced and split into two independent equation systems associated with both bending-transverse shear and extension-twist motions.

The motion equations of bending-transverse shear coupling:

(26)
- k 35 θ z ' ' - k 55 U 2 ' ' - θ y ' + b 1 U ¨ 2 - 2 Ω U ˙ 3 - Ω 2 U 2
            + l = 1 N D b 1 l D U ¨ 2 - 2 Ω U ˙ 3 - Ω 2 U 2 x - x D l + F U 2 b = p y u ,
- k 46 θ y ' ' - k 66 U 3 ' ' - θ z ' + b 1 U ¨ 3 + 2 Ω U ˙ 2 - Ω 2 U 3 + l = 1 N D b 1 l D U ¨ 3 + 2 Ω U ˙ 2 - Ω 2 U 3 x - x D l + F U 3 b = p z u ,
- k 44 θ y ' ' - k 46 U 3 ' ' - θ z ' - k 35 θ z ' - k 55 U 2 ' - θ y   - b 4 θ ¨ y + l = 1 N D - b 4 l D θ ¨ y   Δ x - x D l = 0 ,
- k 33 θ z ' ' - k 35 U 2 ' ' - θ y ' - k 46 θ y ' - k 66 U 3 ' - θ z   - b 5 θ ¨ z + l = 1 N D - b 5 l D θ ¨ z   Δ x - x D l = 0 .

The motion equations of extension-twist coupling:

(27)
- k 11 U 1 ' ' - k 12   ϕ ' ' + b 1 U ¨ 1 = 0 ,
- k 12 U 1 ' ' - k 22 ϕ ' ' + b 4 + b 5 ϕ ¨ - Ω 2 ϕ = 0 .

By using the coordinate transformation Eq. (26) can be transformed to the inertia reference frame X,Y,Z. The results are not presented in this paper for the sake of simplicity.

In present study emphasis is placed on only the problem involving the bending-transverse shear coupling.

2.5. Approximate solution method

In order to find the approximate solution of the rotating composite shaft, the quantities U2x,t, U3x,t, θy(x,t) and θz(x,t) are assumed in the form:

(28)
U 2 x , t = j = 1 N U j t α j ( x ) ,     U 3 x , t = j = 1 N U j t α j x ,
θ y x , t = j = 1 N Θ j t ψ j x ,     θ z x , t = j = 1 N Θ j t ψ j ( x ) ,

where αjx and ψj(x) denote mode shape functions which fulfill all the boundary conditions of the composite shaft, Ujt and Θjt are the generalized coordinates.

Substituting Eq. (28) into the governing Eq. (26) and applying Galerkin’s procedure, the following governing equations in matrix form can be found:

(29)
M X ¨ + C X ˙ + K Ω + K 1 X = f ,

where M denotes the mass matrix, C the gyroscopic matrix, K1 the elastic stiffness matrix, KΩ stiffness matrix which depends on the rotational speed Ω, f the generalized force vector. X¨,X˙ and X denote the generalized acceleration, velocity and displacement vectors respectively.

These matrices can be written as follows:

(30)
M = b 1 E i j + l = 1 N D b 1 l D E i j D l 0 0 0 0 b 1 E i j + l = 1 N D b 1 l D E i j D l 0 0 0 0 - b 4 F i j - l = 1 N D b 4 l D F i j l 0 0 0 0 - b 5 F i j - l = 1 N D b 5 l D F i j l ,
(31)
C = l = 1 N D C y y b l E i j b l - 2 Ω b 1 E i j - 2 Ω l = 1 N D b 1 l D E i j D l + l = 1 N D C y z b l E i j b l 0 0 2 Ω b 1 E i j - 2 Ω l = 1 N D b 1 l D E i j D l + l = 1 N D C z y b l E i j b l l = 1 N D C z z b l E i j b l 0 0 0 0 0 0 0 0 0 0 ,
(32)
K Ω = - b 1 Ω 2 E i j - Ω 2 l = 1 N D b 1 l D E i j D l 0 0 0 0 - b 1 Ω 2 E i j - Ω 2 l = 1 N D b 1 l D E i j D l 0 0 0 0 0 0 0 0 0 0 ,
(33)
K 1 = - k 55 I i j + l = 1 N D K y y b l   E i j b l l = 1 N D K y z b l   E i j b l k 55 J i j - k 35 K i j l = 1 N D K z y b l   E i j b l - k 66 I i j + l = 1 N D K z z b l   E i j b l - k 46 K i j k 66 J i j - k 55 O i j - k 46 P i j - k 44 Q i j + k 55 F i j 0 - k 35 P i j - k 66 O i j 0 - k 33 Q i j + k 66 F i j ,
(34)
f T = Ω 2 l = 1 N D b 1 l D e l B i D l ( cos Ω t sin Ω t   0   0 ) 1 × 4 N ,

where:

(35)
E i j = 0 L α i α j d x ,     F i j = 0 L ψ i ψ j d x ,     E i j D l = 0 L α i α j Δ x - x D l d x ,
F i j D l = 0 L ψ i ψ j Δ x - x D l d x ,     E i j b l = 0 L α i α j Δ x - x b l 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 ,     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 ,     B i D l = 0 L α i Δ x - x D l d x ,       ( i ,   j = 1,2 , , N ) .

The above approximate equations of the motion (29) can be easily simplified as the generalized eigenvalue problem. So the free vibration, stability and response of the composite shaft system can be studied by the eigenvalue equation. The resulting equation is not presented for the sake of simplicity.

3. Numerical simulations

3.1. Model verifications

In 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 frequencies are shown in Table 1 for an increasing number of mode shape functions. From Table 1 it can be seen that to obtain accurate results of the first three natural frequencies, no more than six mode shape functions are required. This indicates clearly that the convergence of the present model is quite good.

The natural frequencies of a cantilever composite shaft obtained for the problem without shear deformation using the present model together with those obtained in [14] are shown in Table 2 for different rotating speeds. A perfect agreement of numerical results with those in [14] can be seen.

Table 1. Effect of model number N on natural frequencies (Ω*= 0 and θ= 30°)

N
ω 1 *
ω 2 *
ω 3 *
ω 4 *
ω 5 *
ω 6 *
2
2.99
11.87
-
-
-
-
4
2.98
11.68
25.17
42.41
-
-
6
2.97
11.65
25.12
41.90
60.29
80.02
8
2.97
11.64
25.10
41.84
60.21
79.22
10
2.97
11.64
25.09
41.82
60.18
79.14

Table 2. Comparison of the natural frequencies of a cantilever composite shaft without shear deformation

Ω *
ω 1 *
ω 2 *
ω 3 *
ω 4 *
ω 5 *
ω 6 *
0
3.5160
3.5160
22.0345
22.0345
61.6973
61.6973
Ref. [14]
3.5160
3.5160
22.0340
22.0340
61.6970
61.6970
2
1.5160
5.5160
20.0345
24.0345
59.6973
69.6973
Ref. [14]
1.5160
5.5160
20.0340
24.0340
59.6970
63.6970
3.5
0.0160
7.0160
18.5345
25.5345
58.1973
65.1973
Ref. [14]
0.0000
7.0160
18.5340
25.5340
58.1970
65.1970
4
-
7.5160
18.0345
26.0345
57.6973
65.6973
Ref. [14]
-
7.5160
18.0340
26.0340
57.6970
65.6970
8
-
11.5160
14.0345
3.30345
53.6973
69.6973
Ref. [14]
-
11.5160
14.0340
30.0340
53.6970
69.6970

The variation of natural frequencies vs. rotating speed for simply supported shear shaft is shown in Fig. 2 from the present model with the ones from [7]. The numerical results are given in terms of the normalized natural frequencies and rotating rate that are defined by ω*=ω/ω0, Ω*=Ω/ω0, where the normalizing factor ω0= 138.85 rad/s is the fundamental frequency of the non-rotating shaft with θ=0°.

From this figure it clearly appears that when  Ω*= 0 a single zero-speed mode natural frequency is obtained. As the rotating speed Ω* is increased, due to the gyroscopic effect the natural frequencies ‘split’ into the upper and lower frequency branches. The upper frequency branch is associated with the forward whirling frequency (FW or F) motion. Similarly lower frequency branch is associated with the backward whirling frequency (BW or B) motion. The minimum rotating rate at which the BW frequency becomes zero is referred to as the critical rotating speed, which will cause the dynamically unstable motion of the rotating shaft. In addition to this, it can obviously be seen that the present numerical results are in good agreement with the results presented in [7].

Fig. 2. The natural frequency of a simply supported composite shaft versus rotating speed for different ply angles

 The natural frequency of a simply supported composite shaft versus rotating speed  for different ply angles

3.2. Results and discussion

In the following application the natural frequencies, critical rotating speeds and the unbalance response of a composite shaft system are computed using presented model. The composite rotor system is a composite shaft with one rigid disk supported by two bearings at the ends, the disk is at the mid-point of the shaft as shown in Fig. 3. The shaft geometrical properties are: L= 2.47 m, r= 0.0635 m, h= 0.001321 m [15]. Ratio of length over radius of the shaft is 38.89764 and radius over thickness 48.10606, which are suitable for features of the slender thin-walled shaft and consistent with assumption as used in the present model. The shaft is made by boron-epoxy composite material with parameters as follows: E11=2.11 GPa,E22=24.1 GPa,G12=G13=G23= 6.9 GPa, ν12= 0.36, ρ= 1967.0 kg/m3 [3]. The stacking sequence of the composite shaft is [±θ]5. The properties of the disk and the two bearings are given in Table 3.

Table 3. The properties of disk and bearings [5]

Disk
Bearings
b 1 l D (kg)
2.4364
-
e (10-5 m)
5.0000
-
b 4 l D , b 5 l D (kg m2)
0.1901
-
K y y ,   K z z (107 N/m)
-
1.75
C y y , C z z (107 Ns/m)
-
5.00

Fig. 3. The composite shaft system with bearings and disk

The composite shaft system with bearings and disk

Fig. 4. The Campbell diagram of a composite shaft system (θ= 0°)

The Campbell diagram of a composite shaft system (θ= 0°)

Fig. 5. The Campbell diagram of a composite shaft system (θ= 60°)

The Campbell diagram of a composite shaft system (θ= 60°)

Fig. 6. The Campbell diagram of a composite shaft system (θ= 90°)

 The Campbell diagram of a composite shaft system (θ= 90°)

Based on mode convergence examination it is found that N= 5 gives suitably converged eigenvalues. So for all results given in this paper N= 5 unless otherwise noted.

Figs. 4-6 present the Campbell diagrams which display the variation of the first three frequencies with respect to rotating speed for various fiber ply-angles. The crossover point of the straight line with BW curves yields the critical rotating speed of the composite shaft system. The predicted critical rotating speeds for the fiber ply-angles θ= 0°, 60°, 90° are Ωcr= 1414 rpm, 2477 rpm and 4023 rpm respectively. The results show that non-rotating natural frequencies and the critical rotating speeds increase as fiber ply-angle increases. The reason is that the bending siffnesses k33 and k44 increase when fiber ply-angle increases.

Fig. 7. The Campbell diagram of a composite shaft system (— with shear deformation; - - - without shear deformation)

The Campbell diagram of a composite shaft system  (— with shear deformation; - - - without shear deformation)

Fig. 8. The first two critical speeds of a composite shaft system versus ratio of length over radius

 The first two critical speeds of a composite shaft system versus ratio of length over radius

Fig. 9. The first two critical speeds of a composite shaft system versus ratio of radius over thickness

 The first two critical speeds of a composite shaft system versus ratio of radius over thickness

Fig. 7 presents the Campbell diagrams both with and without transversal shear effects. The critical rotating speeds predicted by the model with and without transversal shear deformation are found to be Ωcr= 2477 rpm and 3081 rpm for the first critical rotating speed and Ωcr= 23542 rpm and 24750 rpm for the second critical rotating speed respectively. The error of 24.38 % in the first critical rotating speed and 5.13 % in the second critical rotating speed can be observed for problems with and without shear effects. As a result, the classical composite thin-walled beams over-estimate the stability of composite shaft system.

Figs. 8 and 9 present the variation of critical rotating speed with respect to the shaft parameters L/r (ratio of length over radius) and r/h (ratio of radius over thickness) for different fiber ply-angles respectively. The results show that critical speeds decrease as L/r increases, whereas they increase as r/h increases.

Fig. 10. The first two natural frequencies of a composite shaft system versus ratio of length over radius (Ω= 20000 rpm)

 The first two natural frequencies of a composite shaft system versus  ratio of length over radius (Ω= 20000 rpm)

Fig. 11. The first two natural frequencies of a composite shaft system versus ratio of radius over thickness (Ω= 20000 rpm)

 The first two natural frequencies of a composite shaft system versus  ratio of radius over thickness (Ω= 20000 rpm)

Figs. 10 and 11 present the variation of the first two natural frequencies with respect to the shaft parameters L/r and r/h for different fiber ply-angles respectively. Because of the influences of rotating speed (Ω= 20000 rpm) there are the upper and lower frequency branches for each whirling mode as shown in these figures.

From Figs. 10 and 11 it can be seen that the variations of natural frequencies with the parameters L/r and r/h are similar to that previously presented for the critical rotating speed. However, from these figures it is evident that rotating speed has much more influence on the higher-mode ones than the lower-mode ones.

Fig. 12. Trajectory of a composite shaft system for a transient unbalance load: a) undercritical rotating speed Ω= 3000 rpm; b) critical rotating speed Ω=Ωcr= 4023 rpm; c) supercritical rotating speed Ω= 5000 rpm

Trajectory of a composite shaft system for a transient unbalance load:  a) undercritical rotating speed Ω= 3000 rpm; b) critical rotating speed Ω=Ωcr= 4023 rpm;  c) supercritical rotating speed Ω= 5000 rpm

a)

Trajectory of a composite shaft system for a transient unbalance load:  a) undercritical rotating speed Ω= 3000 rpm; b) critical rotating speed Ω=Ωcr= 4023 rpm;  c) supercritical rotating speed Ω= 5000 rpm

b)

Trajectory of a composite shaft system for a transient unbalance load:  a) undercritical rotating speed Ω= 3000 rpm; b) critical rotating speed Ω=Ωcr= 4023 rpm;  c) supercritical rotating speed Ω= 5000 rpm

c)

Fig. 12 presents the trajectories of the geometric center of the rigid disk for various rotating speeds (3000, 4023 and 5000 rpm). It clearly shows that when rotating speed is equal to the critical rotating speed Ω=Ωcr= 4023 rpm, a violent whirling response occurs due to resonant vibration. It can also be noted that the forward mode of a supercritical composite shaft system is stable since the internal damping of the composite shaft is not taken into account in the presented model.

Fig. 13. The time response of the displacement at the disk center: (a) U2(L/2,t); (b) U3(L/2,t)

The time response of the displacement at the disk center: (a) U2(L/2,t); (b) U3(L/2,t)

a)

The time response of the displacement at the disk center: (a) U2(L/2,t); (b) U3(L/2,t)

b)

Fig. 13 presents the time response of the displacement at the disk center which corresponds to that previously presented in Fig. 12. It can be seen from Fig. 12 and 13 that the unbalance responses arrive rapidly to the steady state vibration due to the external damping from bearings.

4. Conclusions

This paper deals with the structural modeling and dynamical analysis of rotating composite shaft. An analytical model capable of predicting of the natural frequency, the critical rotating speed and the unbalanced response of composite shaft system has been proposed based on the thin-walled composite beam theory referred to as VAM. Numerical simulations of the effects of various parameters including fiber ply-angle, geometric dimension and the transverse shear deformation on frequencies and critical rotating speeds are performed. The comparative study and the convergence examination of the approximate solution methodology show that the analytical model and numerical method developed in this paper can be used to highlight the dynamical behaviors of rotating composite shaft system. The present model can further be extended to incorporate the effects of internal damping of composite shaft for evaluating the dynamical instabilities due to internal damping. This improvement will be reported in the forthcoming paper.

Acknowledgements

The research is funded by the National Natural Science Foundation of China (Grant Nos. 11272190) and Shandong Provincial Natural Science Foundation of China (Grant Nos. ZR2011EEM031).

References

  1. Zinberg H., Symonds M. F. The development of an advanced composite tail rotor drive shaft. Presented at the 26th Annual National Forum of the American Helicopter Society, Washington, DC, 1970. [Search CrossRef]
  2. Dos Reis H. L. M., Goldman R. B., Verstrate P. H. Thin-walled laminated composite cylindrical tubes; Part III – Critical speed analysis. J. Compos. Technol. Res., Vol. 9, 1987, p. 58-62. [Search CrossRef]
  3. Kim C. D., Bert C. W. Critical speed analysis of laminated composite, hollow drive shafts. Composites Engineering, Vol. 3, 1993, p. 633-643. [Search CrossRef]
  4. Singh S. P., Gupta K. Composite shaft rotordynamic analysis using a layerwise theory. Journal of Sound and Vibration, Vol. 191, Issue 5, 1996, p. 739-756. [Search CrossRef]
  5. Chang M. Y., Chen J. K., Chang C. Y. A simple spinning laminated composite shaft model. International Journal of Solids and Structures, Vol. 41, 2004, p. 637-662. [Search CrossRef]
  6. Gubran H. B. H., Gupta K. The effect of stacking sequence and coupling mechanisms on the natural frequencies of composite shafts. Journal of Sound and Vibration, Vol. 282, 2005, p. 231-248. [Search CrossRef]
  7. Song O., Jeong N.-H., Librescu L. Implication of conservative and gyroscopic forces on vibration and stability of an elastically tailored rotating shaft modeled as a composite thin-walled beam. Journal of Acoustical Society of America, Vol. 109, Issue 3, 2001, p. 972-981. [Search CrossRef]
  8. Rehfield L. W. Design analysis methodology for composite rotor blades. Proc. Seventh DoD/NASA Conf. on Fibrous Composites in Structural Design, AFWAL-TR-85-3094, 1985. [Search CrossRef]
  9. Berdichevsky V. L., Armanios E., Badir A. M. Theory of anisotropic thin-walled closed-cross-section beams. Composites Engineering, Vol. 2, Issue 5-7, 1992, p. 411-432. [Search CrossRef]
  10. Park J. S., Kim J. H. Design and aeroelastic analysis of active twist rotor blades incorporating single crystal macro fiber composite actuators. Composites: Part B, Vol. 39, 2008, p. 1011-1025. [Search CrossRef]
  11. Cesnik C. E. S., Shin S. J. On the modeling of integrally actuated helicopter blades. International Journal of Solids and Structures, Vol. 38, 2001, p. 1765-1789. [Search CrossRef]
  12. Armanios E. A., Badir A. M. Free vibration analysis of anisotropic thin-walled closed-section beams. AIAA J., Vol. 33, Issue 10, 1995, p. 1905-1910. [Search CrossRef]
  13. Song O., Librescu L. Structural modeling and free vibration analysis of rotating composite thin-walled beams. Journal of the American Helicopter Society, 1997, p. 358-369. [Search CrossRef]
  14. Banerjee J. R., Su H. Development of a dynamic stiffness matrix for free vibration analysis of spinning beams. Comput. Struct., Vol. 82, 2004, p. 2189-2197. [Search CrossRef]
  15. Sino R., Baranger T. N., Chatelet E., Jacquet G. Dynamic analysis of a rotating composite shaft. Composites Science and Technology, Vol. 68, 2008, p. 337-345. [Search CrossRef]