# Vibration and Stability of Variable Cross Section Thin-Walled Composite Shafts with Transverse Shear

Authors: Ma Jing-min; Ren Yong-sheng
Journal: Shock and Vibration
 (2015)

---

## Abstract

A dynamic model of composite shaft with variable cross section is presented. Free vibration equations of the variable cross section thin-walled composite shaft considering the effect of shear deformation are established based on a refined variational asymptotic method and Hamilton’s principle. The numerical results calculated by Galerkin method are analyzed to indicate the effects of ply angle, taper ratio, and transverse shear deformation on the first natural frequency and critical rotating speed. The results are compared with those obtained by using finite element package ANSYS and available in the literature using other models.


---

## Body

## 1. Introduction

The shaft is an extremely important part in various mechanical transmissions. The shafts suffer not only static load such as bending moment, torque, and axial force which usually leads to its combined deformation, but also dynamic load caused by effect of inertia force and gyroscope which produces alternating stress and vibrations due to its rotation. Composite material has been widely utilized in the shaft structural design due to its various advantages such as lightweight, high strength, resistance to fatigue, and designability. Compared with metal shafts, composite shafts have outstanding advantages such as llighter, anti-vibration and high temperature resistant, which can reduce the friction loss, improve the fatigue life, and save running costs. Research on dynamic characteristics of anisotropic composite shaft is a focus of composite material structural dynamics. One of the studies of the laminated composite beams considering cross section warping is done by Librescu et al. [1–6]. In their works, some theories are proposed and applied to analyze the dynamic characteristics of uniform cross section, variable cross section, and containing intelligent material beams. The effects of transverse shear deformation, torsion warping, and the wall thickness have been considered in their theory. However, the warping caused by stretch-bending and torsion-bending coupled deformations is not included. A dynamic model of variable cross section laminated composite shafts based on the theory of Timoshenko beam is presented by Kim et al. [7–10]. The proposed model is of high precision for thick wall shaft due to using the shear factor calculated by the ratio of inner and outer diameter. Another high precision model of thin long beams considering cross section warping caused by extension-bending and torsion-bending deformation is presented by Armanios et al. [11–13]. A displacement function of that beam is proposed based on variational asymptotic method and vibration differential equations are derived and solved by using exact method. However, the propose function does not consider the effect of shear deformation. The further works are done by Yongsheng et al. by using a refined variational asymptotic method. The theory and model involve the free vibration and stability of conventional composite shaft [14] and the composite shaft with internal damping. In these studies, the effect of transverse shear deformation has been incorporated in theoretical formulation. In this paper, a new dynamic model for the analysis of rotating variable cross section composite shaft is presented. The equations of motion are derived using a refined variational asymptotic method and Hamilton’s principle. The free vibration characteristics and stability of variable cross section composite shaft are analyzed. The numerical results are compared with the results obtained by other models.


## 2. Spinning Shaft Model

The rotating composite shaft model is shown in Figure1. The shaft rotates about its axis with constant angular velocity (Ω). d, h, rT, rR, and L represent the outer diameter, wall thickness, the radius of tip middle cross section contour line, the radius of root middle cross section contour line, and the length of the shaft, respectively. (X,Y,Z) represents the inertial coordinate system, and the corresponding unit vector is (I,J,K); (x,y,z) represents the coordinate systems spinning with shaft and the corresponding unit vector is (i,j,k); (ξ,s,x) is a local coordinate system and s represents the axis along the tangent of the cross section middle contour line, while ξ represents the axis along the normal. θ is ply angle.Figure 1
Coordinate systems and geometry of thin-walled rotating shaft.


## 3. Governing Equations

### 3.1. Strain Energy of the Rotating Shaft

A displacement function based on variational asymptotic method is derived and utilized to analyze the dynamic characteristics of composite thin-walled closed section structures in [11]; however, transverse shear deformation is neglected in the displacement function. Considering transverse shear deformation, displacement function can be described as the following form [14]:(1)u1x,s,t=U1x,t-ysθyx,t-zθzx,t+gs,x,t,u2x,s,t=U2x,t-zφx,t,u3x,s,t=U3x,t+yφx,t,(2)θyx,t=U2′-2γzx,θzx,t=U3′-2γyx,where u1, u2, u3 are displacements of an arbitrary point on the shaft along x, y, z directions, respectively. U1(x,t), U2(x,t), and U3(x,t) denote the cross mean displacement along x-, y-, and z-axis. φ(x,t), θy(x,t), and θz(x,t) denote the twist angle and bending angle about y-, and z-axis, respectively. γzx and γyx are the transverse shear strains in planes xz and xy, respectively. g(s,x,t) is a warping displacement function, which is modified as [15] (3)gs,x,t=Gsφ′x,t+g1sU1′x,t+g2sθy′x,t+g3sθz′x,t,where g1s, g2s, g3s, and G(s) are related to physical behavior of the axial strain, bending curvatures, and the torsion twist rate, respectively. All of the ·′ denote the differentiations with respect to x.According to (1)~(3), the strains of the composite shaft are obtained as(4)γxx=U1′x,t-yθy′x,t-zθz′x,t,2γxs=dgds+rnφ′+U2′x,t-θyx,tdyds+U3′x,t-θzx,tdzds,2γxξ=U2′x,t-θyx,tdzds-U3′x,t-θzx,tdyds,where rn is the normal projection ofr which is the position vector of an arbitrary point on the cross section of the deformed shaft in the normal direction:(5)rn=ydzds-zdyds.The strain energy of the composite shaft U is(6)U=12∫0L∬Aσxxεxx+σxsεxs+σxξεxξdAdx,where σxx, σxs, and σxξ are the engineering stresses. Andεxx, εxs, and εxξ are associated engineering strains.Taking the variation of (6), one obtains(7)δUs=∫0L∫Aσxxδεxx+σxsδεxs+σxξδεxξdAdx=∫0LFx′δU1+Qy′δU2+Qz′δU3+Mx′δφ+Mz′+Qyδθy+My′+Qzδθzdx.Stress resultantsFx, Qy, Qz and stress couples Mx, My, Mz are as follows:(8)Fx=∮ΓNxxds,Mx=∮ΓNxxrnds,My=-∮ΓNxxzds,Mz=-∮ΓNxxyds,Qy=-∮ΓNxsdyds+Nxξdydsds,Qz=-∮ΓNxsdyds-Nxξdydsds,where Nxx, Nxs, and Nxξ are shell stress resultants and are of the following form:(9)
N
x
x
N
x
s
N
x
ξ
=
A
s
B
s
2
0
B
s
2
C
s
4
0
0
0
D
s
γ
x
x
2
γ
x
s
γ
x
ξ
,where(10)As=A-11-A-122A-22,Bs=2A-16-A-12A-26A-22,Cs=4A-66-A-262A-22,Ds=κA¯44-A¯452A¯55,A-ij=∑k=1NQ-ijzk-zk-1,i,j=1,2,6;i,j=4,5.Parameters A(s) and B(s) denote the reduced axial and coupling stiffness, while parameters C(s) and D(s) are the reduced shear stiffness, A-ij are the in-plane stiffness components, Q-ij are transformed stiffness of each layer, N is layer numbers, and (zk-zk-1) is thickness of the kth layer. κ is shear factor of the cross section which changes with the cross section and material properties and is of the following form [16]:(11)κ=6Exx1-m¯41+m¯2Gxyυxy2m¯6+18m¯4-18m¯2-2-Exx7m¯6+27m¯4-27m¯2-7,where m¯=Rin/Rout and Rin is the inner radius, while Rout is the outer radius. Exx, υxy, and Gxy are elastic modulus, Poisson’s ratio, and shear modulus.Combining (4), (9), and (10), (8) can be altered as(12)Fx=C11U1′+C12φ′+C13θz′+C14θy′+C15U2′-θy+C16U3′-θz,Mx=C21U1′+C22φ′+C23θz′+C24θy′+C25U2′-θy+C26U3′-θz,My=C31U1′+C32φ′+C33θz′+C34θy′+C35U2′-θy+C36U3′-θz,Mz=C41U1′+C42φ′+C43θz′+C44θy′+C45U2′-θy+C46U3′-θz,Qy=C51U1′+C52φ′+C53θz′+C54θy′+C55U2′-θy+C56U3′-θz,Qz=C61U1′+C62φ′+C63θz′+C64θy′+C65U2′-θy+C66U3′-θz,where Cij=Cji(i,j=1,…,6) are the equivalent cross section stiffness coefficients of the composite shaft, which show the cross section geometry and material properties as(13)C11=∮ΓA-B2Cds+∮ΓB/Cds2∮Γ1/Cds,C12=∮ΓB/Cds∮Γ1/CdsAe,C13=-∮ΓA-B2Czds-∮ΓB/Cds∮ΓB/Czds∮Γ1/Cds,C14=-∮ΓA-B2Cyds-∮ΓB/Cds∮ΓB/Cyds∮Γ1/Cds,C15=12∮ΓBdydsds,C16=12∮ΓBdzdsds,C22=1∮Γ1/CdsAe2,C23=-∮ΓB/Czds∮Γ1/CdsAe,C24=-∮ΓB/Cyds∮Γ1/CdsAe,C25=14∮ΓrnCdydsds,C26=14∮ΓrnCdzdsds,C33=∮ΓA-B2Cz2ds+∮ΓB/Czds2∮Γ1/Cds,C34=∮ΓA-B2Cyzds+∮ΓB/Cyds∮ΓB/Czds∮Γ1/Cds,C35=-12∮ΓBzdydsds,C36=-12∮ΓBzdzdsds,C44=κ∮ΓA-B2Cy2ds+∮ΓB/Cyds2∮Γ1/Cds,C45=-12κ∮ΓBydydsds,C46=-12κ∮ΓBydzdsds,C55=κ∮Γ14Cdyds2+Ddzds2ds,C56=κ∮Γ14C-Ddydsdzdsds,C66=κ∮Γ14Cdzds2+Ddyds2ds,Ae=12∮Γydzds-zdydsds,where ∮Γ(·)ds denotes the integral around the loop of the midline cross section. Ae is the surrounded area of section’s midline. Compared with [12], 16 of the 36 stiffness coefficients are consistent with the ones without shear deformation, while the remaining 20 are new, caused by shear deformation.


### 3.2. Kinetic Energy of Spinning Shaft

The position vector of an arbitrary point on the spinning shaft is written as(14)r=y+u2i+z+u3j+x+u1k.From the above equation, the velocity of an arbitrary point in the fixed coordinate system can be given as(15)V=u˙2-Ωz+u3I+u˙3+Ωy+u2J+u˙1K.The kinetic energy can be written as(16)T=12∫0L∬AρV·VdAdx.Substituting (15) into (16) and omitting the ·′ items due to the negligible effect of warping displacement, the expression of kinematic energy is obtained.Taking variation of the kinematic energy yields(17)δT=-∫0LI1δU1+I2δU2+I3δU3+I4δθy+I5δθz+I6δφ,where(18)I1=mcU¨1-Szθ¨y-Syθ¨z,I2=mcU¨2-2ΩU˙3-Ω2U2-Sz2Ωφ+Ω2-Syφ¨-Ω2φ,I3=mcU¨3+2ΩU˙2-Ω2U3+Szφ¨-Ω2φ-Sy2Ωφ+Ω2,I4=SzU¨1-Izθ¨y-Iyzθ¨z,I5=SyU¨1-Iyzθ¨y-Iyθ¨z,I6=SzU¨3+2ΩU˙2-Ω2U3+Iy+Izφ¨-Ω2φ-SyU¨2-2ΩU˙3-Ω2U2,mc=∬AρdA,Sz=∬AρydA,Sy=∬AρzdA,Iz=∬Aρy2dA,Iy=∬Aρz2dA,Iyz=∫AρyzdA.And Ii(i=1,…,6) is equivalent cross section mass coefficients.


### 3.3. Governing Equations

Hamilton principle is of the following form:(19)∫t0t1δU-δTdt=0.Combining (7), (17), and (19) and considering the variable cross section of the shaft, the governing equations are obtained as(20)-Fx′+I1=0,-Qy′+I2=0,-Qz′+I3=0,-Mx′+I6=0,-Mz′-Qy+I4=0,-My′-Qz+I5=0.Substituting internal forces into (20), the free vibration equations are derived, which can be used to analyze the thin-walled shaft with arbitrary cross shape and variable cross section. The special properties of circular cross section shafts with circumferentially uniform stiffness configuration (CUS) [12] lead various components of stiffness coefficient to be zero. So the vibration equations can be simplified as(21)L1U1,φ=-C11U1′′-C12φ′′+mcU¨1=0,(22)L2U1,φ=-C12U1′′-C22φ′′+Iy+Izφ¨-Ω2φ=0,(23)L3U2,U3,θy,θz=-C35θz′′-C55U2′-θy′+mcU¨2-2ΩU˙3-Ω2U2=0,(24)L4U2,U3,θy,θz=-C46θy′′-C66U3′-θz′+mcU¨3+2ΩU˙2-Ω2U3=0,(25)L5U2,U3,θy,θz=-C44θy′′-C46U3′-θz′-C35θz′-C55U2′-θy-Izθ¨y=0,(26)L6U2,U3,θy,θz=-C33θz′′-C35U2′-θy′-C46θy′-C66U3′-θz-Iyθ¨z=0.Equations (23)~(26) are similar to (10a–d) in [6], and the vibration equations of circular cross section shafts with CUS configuration are decoupled as tensile-torsion coupling system ((21) and (22)) and bending-transverse shear coupling system ((23)~(26)). The latter one is emphasis in this paper. Neglecting the effect of transverse shear, the bending-transverse shear coupling system can be simplified due to θy=U2′, θz=U3′. In this case, the vibration equations become(27)C44U2′′′′+IzU¨2′′+mcU¨2-2ΩU˙3-Ω2U2=0,C33U3′′′′+IyU¨3′′+mcU¨3-2ΩU˙2-Ω2U3=0.The displacement boundary conditions of free vibration cantilever shaft are as(28)U20,t=U30,t=θy0,t=θz0,t=0.The corresponding force boundary conditions for circular cross section shafts with circumferentially uniform stiffness configuration (CUS) are as(29)QyL=0⟹B3U2,U3,θy,θz=C35Lθz′L+C55LU2′L-θyL=0,(30)QzL=0⟹B4U2,U3,θy,θz=C46Lθy′L+C66LU3′L-θzL=0,(31)MzL=0⟹B5U2,U3,θy,θz=C44Lθy′L+C46LU3′L-θzL=0,(32)MyL=0⟹B6U2,U3,θy,θz=C33Lθz′L+C35LU2′L-θyL=0.


### 3.4. Approximate Solution

Because the classical vibration mode function of cantilever shaft satisfies the displacement boundary condition, it does not satisfy the force condition, so General Galerkin method is used to solve vibration and stability of continuous system. Considering the force boundary conditions, and eliminating the shaft’s spatial variables by using General Galerkin method based on the assuming mode shape, the free vibration equations are reduced to be ordinary differential equations with respect to time.Assume that the axial, torsion, and bending deformation is of the following form:(33)U2=∑j=1NU1jtψjx,U3=∑j=1NU3jtψjx,θy=∑j=1Nθyjtφjx,θz=∑j=1Nθzjtφjx,where ψj(x) and φj(x) denote the mode shape of bending deflection and bending angle, respectively, and are taken to be the mode shapes of a uniform, nonrotating, isotropic, fixed-free Euler-Bernoulli beam; namely,(34)φj=sin⁡2j-1π2Lx,cos⁡βjcosh⁡βj=-1,λj=-cos⁡βj+cosh⁡βjsin⁡βj+sinh⁡βj,j=1,2,…N,ψj=cos⁡βjxL-cosh⁡βjxL+λjsin⁡βjxL-sinh⁡βjxL.Substituting (33) into (23), we get(35)L3∑j=1NU2jtψjx,∑j=1NU3jtψjx,∑j=1Nθyjtφjx,∑j=1Nθzjtφjx=0.Multiply (35) by the mode shape, and integrate it on the whole length of the shaft. The differential equation residual is (36)R1m=∫0LL3ψidx.Similar to (35), (29) is multiplied by ψ(L), and the boundary residual is (37)R1Bm=B3∑j=1NU2jtψjL,∑j=1NU3jtψjL,∑j=1NθyjtφjL,∑j=1NθzjtφjLψiL.Based on General Galerkin method, it follows that(38)R1m+R1Bm=0,m=1,…,N.Similarly, the remaining equations give(39)R2m+R2Bm=0,R3m+R3Bm=0,R4m+R4Bm=0,m=1,…,N,where (40)R2m=∫0LL4ψidx,R3m=∫0LL5φidx,R4m=∫0LL6φidx,R2Bm=B4∑j=1NU2jtψjL,∑j=1NU3jtψjL,∑j=1NθyjtφjL,∑j=1NθzjtφjLψiL,R3Bm=B5∑j=1NU2jtψjL,∑j=1NU3jtψjL,∑j=1NθyjtφjL,∑j=1NθzjtφjLφiL,R4Bm=B6∑j=1NU2jtψjL,∑j=1NU3jtψjL,∑j=1NθyjtφjL,∑j=1NθzjtφjLφiL.The final end product of the procedure is a set of equations and is of the form(41)MU¨2jU¨3jθ¨yjθ¨zj+CU˙2jU˙3jθ˙yjθ˙zj+KU2jU3jθyjθzj=0,where (42)M=mcψiψj0000mcψiψj0000-Izφiφj0000-Izφiφj,C=0-2Ωmcψiψj00-2Ωmcψiψj00000000000,K=-C55ψj′ψi′-mcΩ2ψjψi0C55φjψi′-C35φj′ψi′0-C66ψj′ψi′-mcΩ2ψjψi-C46φj′ψi′C66ψjψi′-C55ψj′φi-C46ψj′φi′-C44φj′φi′+C55φjφiC46φjφi′-C35φj′φi-C35ψj′φi′-C66ψj′φiC35φjφi′-C46φj′φi-C33φj′φi′+C66φjφi.The free vibration characteristics of the variable cross section shaft can be obtained by solving (41) with MATLAB.


## 3.1. Strain Energy of the Rotating Shaft

A displacement function based on variational asymptotic method is derived and utilized to analyze the dynamic characteristics of composite thin-walled closed section structures in [11]; however, transverse shear deformation is neglected in the displacement function. Considering transverse shear deformation, displacement function can be described as the following form [14]:(1)u1x,s,t=U1x,t-ysθyx,t-zθzx,t+gs,x,t,u2x,s,t=U2x,t-zφx,t,u3x,s,t=U3x,t+yφx,t,(2)θyx,t=U2′-2γzx,θzx,t=U3′-2γyx,where u1, u2, u3 are displacements of an arbitrary point on the shaft along x, y, z directions, respectively. U1(x,t), U2(x,t), and U3(x,t) denote the cross mean displacement along x-, y-, and z-axis. φ(x,t), θy(x,t), and θz(x,t) denote the twist angle and bending angle about y-, and z-axis, respectively. γzx and γyx are the transverse shear strains in planes xz and xy, respectively. g(s,x,t) is a warping displacement function, which is modified as [15] (3)gs,x,t=Gsφ′x,t+g1sU1′x,t+g2sθy′x,t+g3sθz′x,t,where g1s, g2s, g3s, and G(s) are related to physical behavior of the axial strain, bending curvatures, and the torsion twist rate, respectively. All of the ·′ denote the differentiations with respect to x.According to (1)~(3), the strains of the composite shaft are obtained as(4)γxx=U1′x,t-yθy′x,t-zθz′x,t,2γxs=dgds+rnφ′+U2′x,t-θyx,tdyds+U3′x,t-θzx,tdzds,2γxξ=U2′x,t-θyx,tdzds-U3′x,t-θzx,tdyds,where rn is the normal projection ofr which is the position vector of an arbitrary point on the cross section of the deformed shaft in the normal direction:(5)rn=ydzds-zdyds.The strain energy of the composite shaft U is(6)U=12∫0L∬Aσxxεxx+σxsεxs+σxξεxξdAdx,where σxx, σxs, and σxξ are the engineering stresses. Andεxx, εxs, and εxξ are associated engineering strains.Taking the variation of (6), one obtains(7)δUs=∫0L∫Aσxxδεxx+σxsδεxs+σxξδεxξdAdx=∫0LFx′δU1+Qy′δU2+Qz′δU3+Mx′δφ+Mz′+Qyδθy+My′+Qzδθzdx.Stress resultantsFx, Qy, Qz and stress couples Mx, My, Mz are as follows:(8)Fx=∮ΓNxxds,Mx=∮ΓNxxrnds,My=-∮ΓNxxzds,Mz=-∮ΓNxxyds,Qy=-∮ΓNxsdyds+Nxξdydsds,Qz=-∮ΓNxsdyds-Nxξdydsds,where Nxx, Nxs, and Nxξ are shell stress resultants and are of the following form:(9)
N
x
x
N
x
s
N
x
ξ
=
A
s
B
s
2
0
B
s
2
C
s
4
0
0
0
D
s
γ
x
x
2
γ
x
s
γ
x
ξ
,where(10)As=A-11-A-122A-22,Bs=2A-16-A-12A-26A-22,Cs=4A-66-A-262A-22,Ds=κA¯44-A¯452A¯55,A-ij=∑k=1NQ-ijzk-zk-1,i,j=1,2,6;i,j=4,5.Parameters A(s) and B(s) denote the reduced axial and coupling stiffness, while parameters C(s) and D(s) are the reduced shear stiffness, A-ij are the in-plane stiffness components, Q-ij are transformed stiffness of each layer, N is layer numbers, and (zk-zk-1) is thickness of the kth layer. κ is shear factor of the cross section which changes with the cross section and material properties and is of the following form [16]:(11)κ=6Exx1-m¯41+m¯2Gxyυxy2m¯6+18m¯4-18m¯2-2-Exx7m¯6+27m¯4-27m¯2-7,where m¯=Rin/Rout and Rin is the inner radius, while Rout is the outer radius. Exx, υxy, and Gxy are elastic modulus, Poisson’s ratio, and shear modulus.Combining (4), (9), and (10), (8) can be altered as(12)Fx=C11U1′+C12φ′+C13θz′+C14θy′+C15U2′-θy+C16U3′-θz,Mx=C21U1′+C22φ′+C23θz′+C24θy′+C25U2′-θy+C26U3′-θz,My=C31U1′+C32φ′+C33θz′+C34θy′+C35U2′-θy+C36U3′-θz,Mz=C41U1′+C42φ′+C43θz′+C44θy′+C45U2′-θy+C46U3′-θz,Qy=C51U1′+C52φ′+C53θz′+C54θy′+C55U2′-θy+C56U3′-θz,Qz=C61U1′+C62φ′+C63θz′+C64θy′+C65U2′-θy+C66U3′-θz,where Cij=Cji(i,j=1,…,6) are the equivalent cross section stiffness coefficients of the composite shaft, which show the cross section geometry and material properties as(13)C11=∮ΓA-B2Cds+∮ΓB/Cds2∮Γ1/Cds,C12=∮ΓB/Cds∮Γ1/CdsAe,C13=-∮ΓA-B2Czds-∮ΓB/Cds∮ΓB/Czds∮Γ1/Cds,C14=-∮ΓA-B2Cyds-∮ΓB/Cds∮ΓB/Cyds∮Γ1/Cds,C15=12∮ΓBdydsds,C16=12∮ΓBdzdsds,C22=1∮Γ1/CdsAe2,C23=-∮ΓB/Czds∮Γ1/CdsAe,C24=-∮ΓB/Cyds∮Γ1/CdsAe,C25=14∮ΓrnCdydsds,C26=14∮ΓrnCdzdsds,C33=∮ΓA-B2Cz2ds+∮ΓB/Czds2∮Γ1/Cds,C34=∮ΓA-B2Cyzds+∮ΓB/Cyds∮ΓB/Czds∮Γ1/Cds,C35=-12∮ΓBzdydsds,C36=-12∮ΓBzdzdsds,C44=κ∮ΓA-B2Cy2ds+∮ΓB/Cyds2∮Γ1/Cds,C45=-12κ∮ΓBydydsds,C46=-12κ∮ΓBydzdsds,C55=κ∮Γ14Cdyds2+Ddzds2ds,C56=κ∮Γ14C-Ddydsdzdsds,C66=κ∮Γ14Cdzds2+Ddyds2ds,Ae=12∮Γydzds-zdydsds,where ∮Γ(·)ds denotes the integral around the loop of the midline cross section. Ae is the surrounded area of section’s midline. Compared with [12], 16 of the 36 stiffness coefficients are consistent with the ones without shear deformation, while the remaining 20 are new, caused by shear deformation.


## 3.2. Kinetic Energy of Spinning Shaft

The position vector of an arbitrary point on the spinning shaft is written as(14)r=y+u2i+z+u3j+x+u1k.From the above equation, the velocity of an arbitrary point in the fixed coordinate system can be given as(15)V=u˙2-Ωz+u3I+u˙3+Ωy+u2J+u˙1K.The kinetic energy can be written as(16)T=12∫0L∬AρV·VdAdx.Substituting (15) into (16) and omitting the ·′ items due to the negligible effect of warping displacement, the expression of kinematic energy is obtained.Taking variation of the kinematic energy yields(17)δT=-∫0LI1δU1+I2δU2+I3δU3+I4δθy+I5δθz+I6δφ,where(18)I1=mcU¨1-Szθ¨y-Syθ¨z,I2=mcU¨2-2ΩU˙3-Ω2U2-Sz2Ωφ+Ω2-Syφ¨-Ω2φ,I3=mcU¨3+2ΩU˙2-Ω2U3+Szφ¨-Ω2φ-Sy2Ωφ+Ω2,I4=SzU¨1-Izθ¨y-Iyzθ¨z,I5=SyU¨1-Iyzθ¨y-Iyθ¨z,I6=SzU¨3+2ΩU˙2-Ω2U3+Iy+Izφ¨-Ω2φ-SyU¨2-2ΩU˙3-Ω2U2,mc=∬AρdA,Sz=∬AρydA,Sy=∬AρzdA,Iz=∬Aρy2dA,Iy=∬Aρz2dA,Iyz=∫AρyzdA.And Ii(i=1,…,6) is equivalent cross section mass coefficients.


## 3.3. Governing Equations

Hamilton principle is of the following form:(19)∫t0t1δU-δTdt=0.Combining (7), (17), and (19) and considering the variable cross section of the shaft, the governing equations are obtained as(20)-Fx′+I1=0,-Qy′+I2=0,-Qz′+I3=0,-Mx′+I6=0,-Mz′-Qy+I4=0,-My′-Qz+I5=0.Substituting internal forces into (20), the free vibration equations are derived, which can be used to analyze the thin-walled shaft with arbitrary cross shape and variable cross section. The special properties of circular cross section shafts with circumferentially uniform stiffness configuration (CUS) [12] lead various components of stiffness coefficient to be zero. So the vibration equations can be simplified as(21)L1U1,φ=-C11U1′′-C12φ′′+mcU¨1=0,(22)L2U1,φ=-C12U1′′-C22φ′′+Iy+Izφ¨-Ω2φ=0,(23)L3U2,U3,θy,θz=-C35θz′′-C55U2′-θy′+mcU¨2-2ΩU˙3-Ω2U2=0,(24)L4U2,U3,θy,θz=-C46θy′′-C66U3′-θz′+mcU¨3+2ΩU˙2-Ω2U3=0,(25)L5U2,U3,θy,θz=-C44θy′′-C46U3′-θz′-C35θz′-C55U2′-θy-Izθ¨y=0,(26)L6U2,U3,θy,θz=-C33θz′′-C35U2′-θy′-C46θy′-C66U3′-θz-Iyθ¨z=0.Equations (23)~(26) are similar to (10a–d) in [6], and the vibration equations of circular cross section shafts with CUS configuration are decoupled as tensile-torsion coupling system ((21) and (22)) and bending-transverse shear coupling system ((23)~(26)). The latter one is emphasis in this paper. Neglecting the effect of transverse shear, the bending-transverse shear coupling system can be simplified due to θy=U2′, θz=U3′. In this case, the vibration equations become(27)C44U2′′′′+IzU¨2′′+mcU¨2-2ΩU˙3-Ω2U2=0,C33U3′′′′+IyU¨3′′+mcU¨3-2ΩU˙2-Ω2U3=0.The displacement boundary conditions of free vibration cantilever shaft are as(28)U20,t=U30,t=θy0,t=θz0,t=0.The corresponding force boundary conditions for circular cross section shafts with circumferentially uniform stiffness configuration (CUS) are as(29)QyL=0⟹B3U2,U3,θy,θz=C35Lθz′L+C55LU2′L-θyL=0,(30)QzL=0⟹B4U2,U3,θy,θz=C46Lθy′L+C66LU3′L-θzL=0,(31)MzL=0⟹B5U2,U3,θy,θz=C44Lθy′L+C46LU3′L-θzL=0,(32)MyL=0⟹B6U2,U3,θy,θz=C33Lθz′L+C35LU2′L-θyL=0.


## 3.4. Approximate Solution

Because the classical vibration mode function of cantilever shaft satisfies the displacement boundary condition, it does not satisfy the force condition, so General Galerkin method is used to solve vibration and stability of continuous system. Considering the force boundary conditions, and eliminating the shaft’s spatial variables by using General Galerkin method based on the assuming mode shape, the free vibration equations are reduced to be ordinary differential equations with respect to time.Assume that the axial, torsion, and bending deformation is of the following form:(33)U2=∑j=1NU1jtψjx,U3=∑j=1NU3jtψjx,θy=∑j=1Nθyjtφjx,θz=∑j=1Nθzjtφjx,where ψj(x) and φj(x) denote the mode shape of bending deflection and bending angle, respectively, and are taken to be the mode shapes of a uniform, nonrotating, isotropic, fixed-free Euler-Bernoulli beam; namely,(34)φj=sin⁡2j-1π2Lx,cos⁡βjcosh⁡βj=-1,λj=-cos⁡βj+cosh⁡βjsin⁡βj+sinh⁡βj,j=1,2,…N,ψj=cos⁡βjxL-cosh⁡βjxL+λjsin⁡βjxL-sinh⁡βjxL.Substituting (33) into (23), we get(35)L3∑j=1NU2jtψjx,∑j=1NU3jtψjx,∑j=1Nθyjtφjx,∑j=1Nθzjtφjx=0.Multiply (35) by the mode shape, and integrate it on the whole length of the shaft. The differential equation residual is (36)R1m=∫0LL3ψidx.Similar to (35), (29) is multiplied by ψ(L), and the boundary residual is (37)R1Bm=B3∑j=1NU2jtψjL,∑j=1NU3jtψjL,∑j=1NθyjtφjL,∑j=1NθzjtφjLψiL.Based on General Galerkin method, it follows that(38)R1m+R1Bm=0,m=1,…,N.Similarly, the remaining equations give(39)R2m+R2Bm=0,R3m+R3Bm=0,R4m+R4Bm=0,m=1,…,N,where (40)R2m=∫0LL4ψidx,R3m=∫0LL5φidx,R4m=∫0LL6φidx,R2Bm=B4∑j=1NU2jtψjL,∑j=1NU3jtψjL,∑j=1NθyjtφjL,∑j=1NθzjtφjLψiL,R3Bm=B5∑j=1NU2jtψjL,∑j=1NU3jtψjL,∑j=1NθyjtφjL,∑j=1NθzjtφjLφiL,R4Bm=B6∑j=1NU2jtψjL,∑j=1NU3jtψjL,∑j=1NθyjtφjL,∑j=1NθzjtφjLφiL.The final end product of the procedure is a set of equations and is of the form(41)MU¨2jU¨3jθ¨yjθ¨zj+CU˙2jU˙3jθ˙yjθ˙zj+KU2jU3jθyjθzj=0,where (42)M=mcψiψj0000mcψiψj0000-Izφiφj0000-Izφiφj,C=0-2Ωmcψiψj00-2Ωmcψiψj00000000000,K=-C55ψj′ψi′-mcΩ2ψjψi0C55φjψi′-C35φj′ψi′0-C66ψj′ψi′-mcΩ2ψjψi-C46φj′ψi′C66ψjψi′-C55ψj′φi-C46ψj′φi′-C44φj′φi′+C55φjφiC46φjφi′-C35φj′φi-C35ψj′φi′-C66ψj′φiC35φjφi′-C46φj′φi-C33φj′φi′+C66φjφi.The free vibration characteristics of the variable cross section shaft can be obtained by solving (41) with MATLAB.


## 4. Model Validation and Discussion and Free Vibration Analysis

With the application of the vibration differential equations derived in this paper, the variable cross section shaft presented in [10] is analyzed in our work, and the numerical results are compared with those obtained by [10] and finite element method, as presented in Figure 2. It can be seen that the variation trends of the first natural frequencies with the taper ratio (taper is defined in [10]) obtained by the three methods are consistent. Our model is developed based on variational asymptotic method [11–14] which is exact only for thin-walled and slender shafts, and yet it appears that Kim’s model [10] is available for both thin-walled and thick-walled shafts. So Figure 2 shows obvious difference between them using the geometries for thick walled shafts in [10]. As a result, the shaft in Table 1 is analyzed by the present model and compared with Kim’s model.Table 1
Properties of the analyzed shaft.
r
R
+
r
T
=
0.02 m
[θ/−θ]3 Layup
Length (L) = 0.5 m
Ply thickness = 0.15 mm
v
12
=
v
13
=
0.24
v
23
=
0.2067
E
11 = 192 Gpa
E
22
=
E
33
=
7.24 Gpa
G
12
=
G
13
=
4.07 GPa
G
23
=
3 GPaFigure 2
Changes comparison with TR (Ω=0 rad/s).Figures3-4 and Tables 2-3 illustrate the first natural frequencies of a rotating and a stationary shaft considering or nonconsidering the transverse shear deformation. It is easy to find that the transverse shear deformation has only a little effect when ply angle is in the range of 15°~40°. The largest effect occurs at the θ=0° and TR=2. At the same time, Figures 3 and 4 reveal that shear deformation has same effects on the rotating and stationary shafts.Table 2
Shear stress effects on the first natural frequency (Ω = 0 rad/s).
Taper  ratio
Plyangle (°)
Withshear (Hz)
Withoutshear (Hz)
Changerate (%)
TR = 0
0°
158.76
165
3.8
10°
154.57
157.56
1.9
20°
125.85
126.12
0.2
60°
36.20
34.95
−3.6
TR = 1
0°
209.63
221.98
5.6
10°
205.72
211.99
3.0
20°
168.47
169.68
0.71
60°
48.1
47.02
−2.3
TR = 2
0°
269.77
290.42
7.11
10°
267.90
277.34
3.4
20°
222.80
222.0
−0.36
60°
64.33
61.52
−4.57Table 3
Shear stress effects on the first natural frequency (Ω = 100 rad/s).
Taper  ratio
Plyangle (°)
Withshear (Hz)
Withoutshear (Hz)
Changerate (%)
TR = 0
0°
142.85
149.09
4.2
10°
138.67
141.67
2.1
20°
109.95
110.22
0.2
60°
20.29
19.05
−6.5
TR = 1
0°
193.72
206.08
6
10°
189.80
196.09
3.2
20°
152.55
153.78
0.8
60°
32.20
31.12
−3.4
TR = 2
0°
253.87
274.52
7.5
10°
255.03
271.74
6.1
20°
206.88
206.10
−0.37
60°
48.42
45.62
−6.13Figure 3
The influence of transverse shear deformation (Ω=0 rad/s).Figure 4
The influence of transverse shear deformation (Ω=100 rad/s).Tables4 and 5 list the effect of mode number on the first natural frequencies and also show the effects of ply method. When the ply method is θ/-θ3, the present model and Kim’s model both show good convergence as shown in Table 4. The convergence of Kim’s model is lost when the ply method is θ6, while the present model still shows good convergence.Table 4
Effect of mode numberN on the first natural frequencies ([θ/−θ]3, Ω = 0 rad/s).
N
TR = 0
TR = 1
TR = 2
0°
20°
60°
0°
20°
60°
0°
20°
60°
2
Kim
159.74
125.07
35.34
213.77
168.12
47.31
279.68
222.62
63.31
Present
158.76
125.85
36.29
209.63
168.47
48.10
269.77
222.80
64.33
3
Kim
159.48
124.28
34.93
213.63
167.30
46.92
278.54
219.02
61.54
Present
158.50
125.45
35.14
209.38
168.04
47.09
268.32
219.14
61.97
4
Kim
159.29
124.12
34.86
213.44
167.10
46.87
278.39
218.79
61.36
Present
158.29
125.29
35.02
209.07
167.99
46.99
267.99
219.02
61.59
5
Kim
159.25
124.06
34.82
213.43
167.07
46.85
278.37
218.55
61.27
Present
158.24
125.25
34.96
209.03
167.97
46.99
267.86
218.89
61.46Table 5
Effect of mode numberN on the first natural frequencies ([θ]6, Ω = 0 rad/s).
N
TR = 0
TR = 1
TR = 2
0°
20°
60°
0°
20°
60°
0°
20°
60°
2
Kim
159.74
816.93
36.66
213.77
127.38
52.46
279.68
220.04
70.05
Present
158.76
92.04
34.54
209.63
122.05
46.22
269.77
161.33
61.86
3
Kim
159.48
295.27
141.32
213.63
366.79
174.64
278.54
448.68
32.20
Present
158.50
92.71
34.25
209.38
124.39
45.94
268.32
162.97
60.22
4
Kim
159.29
799.21
49.16
213.44
184.04
70.50
278.39
142.48
94.92
Present
158.29
92.92
34.16
209.07
124.80
45.90
267.99
163.08
60.06
5
Kim
159.25
219.81
216.16
213.43
371.12
17.24
278.37
453.50
39.785
Present
158.24
93.35
34.15
209.03
125.54
45.90
267.86
164.05
60.01Figure6 shows the effect of ply angle on different taper ratio shafts with 0 rotation speed, and it also shows the difference among ANSYS, present model, and Kim’s model. The first natural frequency is reduced obviously with the increase of ply angle, while it increased with the increase of taper ratio. And the three kinds of results are completely consistent with each other.Figures7–9 show the effect of rotating speed on different taper ratio and ply angle shafts. And they also reveal the difference among the three models. When the rotating speed is more than 0, two natural frequencies occur for each mode. The higher value is called the forward precession mode, and the lower one represents a backward precession mode. So 1B and 1F denote the first backward precession frequency and first forward precession, respectively, in all of the figures. The same effects of rotating speed to the first natural frequency can be obtained. The first backward precession frequency (1B) is decreasing with increasing rotating speed for different taper ratio and different ply angle shafts. The first forward precession frequency (1F) is exactly the reverse. In Figures 7–9, the results obtained by the present model are close to the ones obtained by ANSYS software. However, the results obtained by Kim et al. [10] are different from that of ANSYS software obviously; the reason may be that the effect of spinning inertia force has not been considered by Kim et al. [10].Figure10 shows the effect of ply angle on 1B with 100 rad/s rotating speed and also displays the difference among the three methods. The effect of ply angle is similar to Figure 5.Figure 5
The first natural frequency versus ply angle (Ω=0 rad/s).Figure 6
Changes comparison with rotating speed (θ=0°).Figure 7
Changes comparison with rotating speed (θ=20°).Figure 8
Changes comparison with rotating speed (θ=60°).Figure 9
The first natural frequency versus ply angle (Ω=100 rad/s).Figure 10
The first natural frequency versus taper ratio.Figure11 shows the effect of taper ratio on natural frequencies of different taper ratio and ply angle rotating shaft or stationary shaft. With the increase of taper ratio, all of the natural frequencies of the shafts are increased. And the lower value ply angle is more obviously affected than high ones.Figure 11
Campbell diagrams of the shaft withθ=0°.A conclusion that lower ply angle and large taper ratio could increase the first natural frequency could be made from Figures6–11. This conclusion is important for rotating shafts. If rotating speed of the shafts is under the first natural frequency, the shaft can keep dynamic stability. The critical rotating speed could be increased with the increase of first natural frequency, which could improve the transmission efficiency in the fields of car, ship, and plane. If the milling bar is modeled as a rotating cantilever shaft, the machining efficiency could be improved as well.In a word, the analytical model presented in this work is verified and validated to a decisive extent in Figures6–10.


## 5. Stability Analysis

Figures11–13 show the Campbell diagrams of the shafts with different ply angle and different taper ratio. They illustrate that rotating shaft with 0° ply angle and 2 taper ratio is of maximal rotating critical speed for the shafts. Figures 11–13 show that the critical rotating speeds increase with taper ratio and decrease with ply angle.Figure 12
Campbell diagrams of the shaft withθ=20°.Figure 13
Campbell diagrams of the shaft withθ=60°.All the data and figures illustrate that the stability of rotating shaft can be dramatically improved by utilizing variable cross section and appropriate ply angle.


## 6. Conclusions

A refined variational asymptotic method and Hamilton’s principle are used to establish the vibration differential equation of rotating variable cross section composite shaft which is solved by applying Extended Galerkin method. And the effects of ply angle, taper ratio, and spinning speed on vibration characteristics are analyzed. Critical rotating speeds of the shafts with various ply angle and taper ratio are also calculated. The results illustrate that small ply angle and big taper ratio can improve the natural frequencies and critical rotating speeds effectively. The numerical results are compared with the results obtained by other models.



---

*Source: 697493-2015-09-03.xml*