Abstract
The mechanical model for the dynamic behavior of an internally damped rotating composite shaft is derived using a refined variational asymptotic method and the principle of virtual work. The composite shaft is considered as an anisotropic thinwalled Timoshenko beam. The internal damping of composite shaft is modelled by adopting the multiscale damping analysis method. Galerkin’s method is employed to discretize and solve the equations of motion. The effect of design parameters including fiber orientation, length aspect ratio, stacking sequences and boundary conditions on the free vibration and stability of composite shaft is investigated.
1. Introduction
The rotating composite shafts are now used widely in many application areas involving the rotating machinery systems. In the last few years, there exists numerous researches related to predicting natural frequencies and critical speeds of composite shaft [19]. As rotating speed is above the first critical speed, the effect of internal damping on the dynamic stability of rotating shaft becomes increasingly important. The sources of internal damping in rotating shaft can be due to material damping or friction. Composite shafts have higher material damping than conventional metallic shafts. Therefore, accurate prediction of effects of composite internal damping on the dynamical behavior of rotating composite shafts are necessary for improving performance of the shaft system. So far, only less work seems to have been reported towards the stability analysis, particularly effects of internal damping on the stability of a rotating composite shaft [1013]. Singh and Gupta [10] 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 Kim et al. [11], Mazzei and Scott [12] as well as Montagnier and Hochard [13]. 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. [14] 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 modelled by finite element method. However, only the elastic coupling effects induced by symmetrical stacking are taken into account in their mode. Ren et al. [15] presented the dynamical analysis of rotating composite thinwalled shafts with internal damping based on variational asymptotically method (VAM) [16], which is able to represent the crosssectional warping and various elastic couplings by symmetrical or unsymmetrical stacking. Internal damping of the shaft is introduced via the multiscale damping mechanics [17] but transverse shear deformation has not been taken into account in the model.
In the present work, a refined VAM composite thinwalled beam theory which incorporates the effect of transverse shear deformation is used for the dynamical analysis of rotating composite thinwalled shafts with internal damping. The equations of motion of the rotating Timoshenko composite shafts are derived by the principle of virtual work. The internal damping of composite shaft is modelled based on the multiscale damping mechanics [17]. Galerkin’s method is applied to discretize and solve the governing equations. The critical rotating speeds and instability thresholds are obtained through numerical simulations. The effects of fiber orientation, length aspect ratio, stacking sequences and boundary conditions are assessed. The validity of the model is proved by comparing the results with those in literature and convergence examination.
2. Theoretical formulations and numerical method
2.1. The displacement and strain
The symbols $L$, $h$ and $r$ represent the length, wall thickness and radius of the shaft, respectively. The inertial reference system $(X,Y,Z)$ and rotating reference system $(x,y,z)$ will be used for describing the motion of the shaft. The ($I,J,K)$ and $(i,j,k)$ are the unit vectors with respect to the reference system $(X,Y,Z)$ and$(x,y,z)$, respectively. In addition, a curvilinear coordinate system $(\varsigma ,s,x)$ is also used, where and $s$ is the circumferential coordinate measured along the tangent to the middle surface, whereas $\varsigma $ is measured along the normal to the middle surface.
The displacement field incorporating shear deformation in the local frame denoted $(x,s,\xi )$ is assumed in the form [9]:
${u}_{2}\left(x,y,z,t\right)={U}_{2}\left(x,t\right)z\varphi \left(x,t\right),{u}_{3}\left(x,y,z,t\right)={U}_{3}\left(x,t\right)+y\varphi \left(x,t\right),$
where, ${U}_{1}\left(x,t\right)$, ${U}_{2}\left(x,t\right)$, ${U}_{3}(x,t)$ denote the average displacements along the $x$, $y$ and $z$axis, while $\varphi \left(x,t\right)$, ${\theta}_{y}\left(x,t\right)$, ${\theta}_{z}(x,t)$ denote the twist angles about the $x$axis and rotations about the $y$ and $z$axis, respectively.
The warping functions $g(s,x,t)$ can be assumed as follows:
In the above equation, the functions ${g}_{1}\left(s\right)$, ${g}_{2}\left(s\right)$, ${g}_{3}\left(s\right)$, $G\left(s\right)$ associated with physical behavior for the axial strain, the bending curvatures, and the torsion twist rate, respectively. The primes in Eq. (2) denote differentiation with respect to $x$.
In Eqs. (1) and (2), the expressions of ${\theta}_{y}\left(x,t\right)$, ${\theta}_{z}(x,t)$ are in the form as follows:
${\theta}_{z}\left(x,t\right)={U}_{3}^{\text{'}}\left(x,t\right)2{\gamma}_{y1}.$
The strain of the composite shaft can be expressed as [9]:
$2{\gamma}_{1s}=\frac{dg}{ds}+{r}_{n}{\varphi}^{\text{'}}+\left({U}_{2}^{\text{'}}\left(x,t\right){\theta}_{y}\left(x,t\right)\right)\frac{dy}{ds}+\left({U}_{3}^{\text{'}}\left(x,t\right){\theta}_{z}\left(x,t\right)\right)\frac{dz}{ds},$
$2{\gamma}_{x\xi}=\left({U}_{2}^{\text{'}}\left(x,t\right){\theta}_{y}\left(x,t\right)\right)\frac{dz}{ds}\left({U}_{3}^{\text{'}}\left(x,t\right){\theta}_{z}\left(x,t\right)\right)\frac{dy}{ds}.$
The position vector of an arbitrary point on the crosssection of the deformed shaft is $\mathbf{R}=\left(y+{u}_{2}\right)\mathbf{i}+\left(z+{u}_{3}\right)\mathbf{j}+\left(x+{u}_{1}\right)\mathbf{k}$. Here, $\mathbf{i}$, $\mathbf{j}$, $\mathbf{k}$ are unit vectors of the coordinate systems $(x,y,z)$. By taking the time derivatives of unit vectors and adopting the assumption of constant rate $\mathrm{\Omega}\text{,}$ the velocity and acceleration vector of an arbitrary point are:
and
respectively.
Fig. 1Geometry and coordinate systems of composite thinwalled shaft
2.2. Equations of motion
The equations of motion of the shaft can be described by a variational form:
where $\delta {U}_{s}$, $\delta {T}_{s}$ and $\delta {W}_{s}$ are the variation of the strain energy, kinetic energy and the dissipated energy of crosssection, respectively.
2.2.1. The strain energy and dissipated energy
The strain energy and dissipated energy of composite shaft can be described by multiscale modelling for each composite ply, for the laminate, for crosssection and eventually for the complete shaft.
2.2.1.1. Equivalent crosssection stiffness matrix
The strain energy density of a composite ply can be expressed as:
The stressstrain relation is written as:
where, $\left[{\stackrel{}{Q}}_{ij}\right]$ is equivalent offaxis stiffness matrix of a composite ply with respect to the system $osx\varsigma $, $\kappa =1/(2{G}_{12}{\nu}_{12}/{E}_{1})$ is a Timoshenko shear coefficient [18].
Substituting Eq.(7) into (6) and integrate over the thickness $\xi $ to obtain the strain energy of the laminate as following:
$+4{A}_{55}{\gamma}_{\xi s}^{2}+8{A}_{45}{\gamma}_{\xi 1}{\gamma}_{\xi s},$
where the laminate inplane stiffnesses are:
${A}_{ij}={\int}_{h/2}^{h/2}{\stackrel{}{Q}}_{ij}d\xi =2\sum _{k=1}^{N/2}{\stackrel{}{Q}}_{ij}^{k}\left({h}_{k}{h}_{k1}\right),\left(i,j=4,5\right).$
For the case of no internal pressure acting on the shaft, Eq. (8) can be simplified by using free hoop stress resultants assumption (${N}_{ss}=0$, ${N}_{\xi s}=0$) and then integrating the resulted expression around the crosssectional midline to get the strain energy of crosssection in terms of the generalized displacements and the equivalent crosssection stiffness matrix, as follows:
where ${\left\{\delta \right\}}^{T}=\left\{{U}_{1}^{\text{'}}{\varphi}^{\text{'}}{\theta}_{z}^{\text{'}}{\theta}_{y}^{\text{'}}{U}_{2}^{\text{'}}{\theta}_{y}{U}_{3}^{\text{'}}{\theta}_{z}\right\}\text{,}$$\mathbf{K}$ is 6×6 symmetric stiffness matrix, the coefficients ${k}_{ij}$($i$, $j=$ 1, 2, …, 4) represent the stiffness, as described in [16] for the case of an nonshearable thinwalled beam, while the transverse shear deformation contribution to the stiffness coefficients is represented by the new terms ${k}_{ij}$($i=$1, 2,…, 6; $j=$ 5, 6), ${k}_{ji}$ ($i=$1, 2, 3, 4; $j=$5, 6) [9].
Injecting a Timoshenko shear coefficient $\kappa $, the new stiffness coefficient terms ${k}_{ij}$ ($i=$ 1, 2,…, 6; $j=$ 5, 6), ${k}_{ji}$ ($i=$ 1, 2,…,4; $j=$ 5, 6) can be expressed as:
${k}_{26}=\frac{1}{4}\kappa {\oint}_{\mathrm{\Gamma}}^{}{r}_{n}C\frac{dz}{ds}ds,{k}_{35}=\frac{1}{2}\kappa {\oint}_{\mathrm{\Gamma}}^{}Bz\frac{dy}{ds}ds,{k}_{36}=\frac{1}{2}\kappa {\oint}_{\mathrm{\Gamma}}^{}Bz\frac{dz}{ds}ds,$
${k}_{45}=\frac{1}{2}\kappa {\oint}_{\mathrm{\Gamma}}^{}By\frac{dy}{ds}ds,{k}_{46}=\frac{1}{2}\kappa {\oint}_{\mathrm{\Gamma}}^{}By\frac{dz}{ds}ds,$
${k}_{55}=\kappa {\oint}_{\mathrm{\Gamma}}^{}\left[\frac{1}{4}C{\left(\frac{dy}{ds}\right)}^{2}+D{\left(\frac{dz}{ds}\right)}^{2}\right]ds,{k}_{56}=\kappa {\oint}_{\mathrm{\Gamma}}^{}\left(\frac{1}{4}CD\right)\frac{dy}{ds}\frac{dz}{ds}ds,$
${k}_{66}=\kappa {\oint}_{\mathrm{\Gamma}}^{}\left[\frac{1}{4}C{\left(\frac{dz}{ds}\right)}^{2}+D{\left(\frac{dy}{ds}\right)}^{2}\right]ds,$
where:
$D\left(s\right)=\left({A}_{44}\frac{{A}_{45}^{2}}{{A}_{55}}\right),{A}_{ij}=\sum _{k=1}^{N}{\stackrel{}{Q}}_{ij}^{\left(k\right)}\left({z}_{k}{z}_{k1}\right),\left(i,j=1,2,6;ij=4,5\right).$
In above equations, the variables $A\left(s\right)$, $B\left(s\right)$, $C\left(s\right)$ and $D\left(s\right)$ represent the reduced axial, coupling, and shear stiffness of crosssection, respectively. The variables $A\left(s\right)$, $B\left(s\right)$, $C\left(s\right)$ correspond to the stiffness coefficients for nonshearable thinwalled beam [16]. The new variable $D\left(s\right)$ results from the transverse shear deformation.
By performing the integral along the length $L$, the strain energy of the shaft can be expressed as following:
2.2.1.2. Equivalent crosssection damping matrix
The dissipated energy density of a composite ply can be written as:
where ${\eta}_{ij}$ are offaxis damping coefficients of a composite ply.
The dissipated energy of the laminate has following form:
$+2\left({A}_{d26}+{A}_{d62}\right){\gamma}_{1s}{\gamma}_{ss}+4{A}_{d44}{\gamma}_{\xi 1}^{2}+4{A}_{d55}{\gamma}_{\xi s}^{2}+4\left({A}_{d45}+{A}_{d54}\right){\gamma}_{\xi 1}{\gamma}_{\xi s},$
where the laminate dampings inplane are:
${A}_{dij}={\int}_{h/2}^{h/2}{\eta}_{il}{\stackrel{}{Q}}_{lj}d\xi =2\sum _{k=1}^{N/2}{\eta}_{il}^{k}{\stackrel{}{Q}}_{lj}^{k}\left({h}_{k}{h}_{k1}\right),\left(i,j,l=4,5\right),$
where $\left[{\eta}_{ij}\right]={\left[R\right]}^{T}\left[{\eta}_{l}\right]\left[R\right]$, $\left[R\right]$ is transformation matrix. $\left[{\eta}_{l}\right]=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left({\eta}_{l1},{\eta}_{l2},{\eta}_{l4},{\eta}_{l5},{\eta}_{l6}\right)$ is the onaxis damping matrix within a ply.
Similar to the derivation of Eq. (10) for the strain energy of crosssection, the dissipated energy of crosssection can be obtained in terms of generalized displacements and the equivalent crosssection damping matrix as follows:
where $\left[\mathbf{C}\right]$ is 6×6 symmetric damping matrix. The formulation of its components ${c}_{ij}$ is analogous to stiffness components ${k}_{ij}$ as shown in Eqs. (10) and (11), but the terms $A$, $B$, $C$ and $D$ in Eqs. (10) and (11) should be replaced by the terms ${A}_{d}$, ${B}_{d}$, ${C}_{d}$ and ${D}_{d}$, respectively.
The variables ${A}_{d}$, ${B}_{d}$, ${C}_{d}$ and ${D}_{d}$ are given by:
${B}_{d}\left(s\right)={A}_{d16}+{A}_{d61}+2\left(\frac{{A}_{12}{A}_{26}}{{A}_{22}^{2}}\right){A}_{d22}\left(\frac{{A}_{26}}{{A}_{22}}\right)\left({A}_{d12}+{A}_{d21}\right)\left(\frac{{A}_{12}}{{A}_{22}}\right)\left({A}_{d26}+{A}_{d62}\right),$
${C}_{d}\left(s\right)=4\left[{A}_{d66}+\left(\frac{{A}_{26}^{2}}{{A}_{22}^{2}}\right){A}_{d22}\left(\frac{{A}_{26}}{{A}_{22}}\right)\left({A}_{d26}+{A}_{d62}\right)\right],$
${D}_{d}\left(s\right)=4\left[{A}_{d44}+\left(\frac{{A}_{45}^{2}}{{A}_{55}^{2}}\right){A}_{d55}\left(\frac{{A}_{45}}{{A}_{55}}\right)\left({A}_{d45}+{A}_{d54}\right)\right].$
The dissipated energy of the shaft can be expressed as following:
2.2.2. The kinetic energy
The kinetic energy of the crosssection of rotating shaft is:
where $\left[\mathbf{d}\mathbf{i}\mathbf{a}\mathbf{g}\left(\rho \right)\right]$ is a diagonal matrix with components equal to the mass density $\rho $ of a ply.
Substituting the expressions of $\mathbf{R}$ and $\ddot{\mathbf{R}}$ into Eq. (20), the kinetic energy of crosssection can be obtained as:
where:
${I}_{2}={b}_{1}\left({\ddot{U}}_{2}2\mathrm{\Omega}{\dot{U}}_{3}{\mathrm{\Omega}}^{2}{U}_{2}\right){b}_{2}\left(2\mathrm{\Omega}\dot{\varphi}+{\mathrm{\Omega}}^{2}\right){b}_{3}\left(\ddot{\varphi}{\mathrm{\Omega}}^{2}\varphi \right),$
${I}_{3}={b}_{1}\left({\ddot{U}}_{3}+2\mathrm{\Omega}{\dot{U}}_{2}{\mathrm{\Omega}}^{2}{U}_{3}\right)+{b}_{2}\left(\ddot{\varphi}{\mathrm{\Omega}}^{2}\varphi \right){b}_{3}\left(2\mathrm{\Omega}\dot{\varphi}+{\mathrm{\Omega}}^{2}\right),$
${I}_{4}={b}_{2}{\ddot{U}}_{1}{b}_{4}{\ddot{\theta}}_{y}{b}_{6}{\ddot{\theta}}_{z},{I}_{5}={b}_{3}{\ddot{U}}_{3}{b}_{6}{\ddot{\theta}}_{y}{b}_{5}{\ddot{\theta}}_{z},$
${I}_{6}={b}_{2}\left({\ddot{U}}_{3}+2\mathrm{\Omega}{\dot{U}}_{2}{\mathrm{\Omega}}^{2}{U}_{3}\right)+\left({b}_{4}+{b}_{5}\right)\left(\ddot{\varphi}{\mathrm{\Omega}}^{2}\varphi \right){b}_{3}\left({\ddot{U}}_{2}2\mathrm{\Omega}{\dot{U}}_{3}{\mathrm{\Omega}}^{2}{U}_{2}\right),$
${b}_{1}=\oint \rho h\left(s\right)ds,{b}_{2}=\oint \rho yh\left(s\right)ds,{b}_{3}=\oint \rho zh\left(s\right)ds,$
${b}_{4}=\oint \rho {y}^{2}h\left(s\right)ds,{b}_{5}=\oint \rho {z}^{2}h\left(s\right)ds,{b}_{6}=\oint \rho yzh\left(s\right)ds.$
The kinetic energy of the shaft can be expressed as following:
The approximate solutions are assumed for ${U}_{1}\left(x,t\right)$, ${U}_{2}\left(x,t\right)$, ${U}_{3}\left(x,t\right)$, ${\theta}_{y}\left(x,t\right)$,${\theta}_{z}\left(x,t\right)$ and $\varphi \left(x,t\right)$ in the form:
${U}_{3}\left(x,t\right)=\sum _{j=1}^{N}{U}_{3j}{\alpha}_{j}\left(x\right){e}^{i\lambda t},{\theta}_{y}\left(x,t\right)=\sum _{j=1}^{N}{\mathrm{\Theta}}_{yj}{\psi}_{j}\left(x\right){e}^{i\lambda t},$
${\theta}_{z}\left(x,t\right)=\sum _{j=1}^{N}{\mathrm{\Theta}}_{zj}{\psi}_{j}\left(x\right){e}^{i\lambda t},\varphi \left(x,t\right)=\sum _{j=1}^{N}{\mathrm{\Phi}}_{j}{\theta}_{j}\left(x\right){e}^{i\lambda t},$
where ${\beta}_{j}\left(x\right)$, ${\alpha}_{j}\left(x\right)$, ${\psi}_{j}\left(x\right)$ and ${\theta}_{j}\left(x\right)$ represent mode shape functions. These functions satisfy the geometric boundary conditions at $x=0$ and $x=L$, $\lambda $ is complex eigenvalues of the system, ${U}_{1j}$, ${U}_{2j}$, ${U}_{3j}$, ${\mathrm{\Theta}}_{yj}$, ${\mathrm{\Theta}}_{zj}$ and ${\mathrm{\Phi}}_{j}$ are undetermined constants, $i=\sqrt{1}$.
Substituting Eq. (24) into the energy expressions (13), (19), (23) and (5), and using Galerkin method, one obtains the motion equations in matrix form:
where ${\left\{U\right\}}^{T}=\left(\dots {U}_{1j}\dots {U}_{2j}\dots {U}_{3j}\dots {\mathrm{\Theta}}_{yj}\dots {\mathrm{\Theta}}_{zj}\dots {\mathrm{\Phi}}_{j}\dots \right)$ is a constant vector; $\left[\mathbf{M}\right]$ the mass matrix; $\left[\mathbf{G}\right]$ the gyroscopic matrix, $\left[\mathbf{C}\right]$ the damping matrix, $\left[\mathbf{K}\right]$ the stiffness matrix, $\left[{\mathbf{K}}_{\mathrm{\Omega}}\right]$ the stiffness matrix which depends on the rotational speed $\mathrm{\Omega}$. The detailed expressions of these matrices are follows:
where:
${K}_{ij}={\int}_{0}^{L}{\alpha}_{i}{\psi}_{j}^{\text{'}\text{'}}dx,{O}_{ij}={\int}_{0}^{L}{\psi}_{i}{\alpha}_{j}^{\text{'}}dx,{P}_{ij}={\int}_{0}^{L}{\psi}_{i}{\alpha}_{j}^{\text{'}\text{'}}dx,{Q}_{ij}={\int}_{0}^{L}{\psi}_{i}{\psi}_{j}^{\text{'}\text{'}}dx,$
${B}_{ij}={\int}_{0}^{L}{\beta}_{i}{\beta}_{j}dx,{N}_{ij}={\int}_{0}^{L}{\beta}_{i}{\beta}_{j}^{\text{'}\text{'}}dx,{D}_{ij}={\int}_{0}^{L}{\beta}_{i}{\theta}_{j}^{\text{'}\text{'}}dx,{G}_{ij}={\int}_{0}^{L}{\theta}_{i}{\theta}_{j}dx,$
${H}_{ij}={\int}_{0}^{L}{\theta}_{i}{\beta}_{j}^{\text{'}\text{'}}dx,{L}_{ij}={\int}_{0}^{L}{\theta}_{i}{\theta}_{j}^{\text{'}\text{'}}dx,\left(i,j=1,\dots ,N\right).$
From Eq. (25), one can obtain the following complex eigenvalue problem:
The complex eigenvalue $\lambda $ of Eq.(32) can be written as:
The imaginary part $\omega $ and the real part $\sigma $ of $\lambda $ give the damping natural frequency and the decay, respectively. If $\sigma $ is negative, the amplitude decays in time and the shaft has a stable behavior because the whirl motion tends to reduce its amplitude. On the other hand, if $\sigma $ is positive, the motion is unstable with the amplitude growing exponentially in time.
It should be emphasized that the equation (25) is obtained using the rotating reference system $(x,y,z)$. Besides, when the ply angle distribution is a CUS (Circumferentially Uniform Stiffness) configuration [20], Eq. (25) can be decoupled into the flapping bendingsweeping bending equation and extensiontwist equation. This paper will address the flapping bendingsweeping bending coupling and present all numerical results with reference to the rotating reference system. For simplicity, the resulting flapping bendingsweeping bending equation are not presented.
3. Numerical results
The natural frequency and modal damping of thinwalled beam with circular and box sections are predicted and shown in Table 1 and 2. The numerical results of present model are compared with those using shear beam finite element model [20] and nonshearable beam model [17]. As can been seen, the present model shows good agreement with those found in the literature.
The numerical calculations are performed by considering the composite shaft whose elastic characteristics are listed in Table 3. The shaft has circular crosssection of mean diameter $d=$0.352 m. Ply thickness is 0.000635 m, and the wall has 16 layers. The material properties of each ply are listed in Table 3.
Figs. 2 show the first three natural frequencies vs. rotating speed for simply supported shear shaft and for various ply angles. From these figures it can be seen that when rotating speed equals to zero a single nonrotating 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 motion. Similarly, lower frequency branch is associated with the backward whirling motion. The minimum rotating speed at which the lower frequency becomes zero is referred to as the critical rotating speed. It is also seen that with the decrease of ply angle the nonrotating natural frequencies increase. This is because the closer the fiber is oriented to the longitudinal direction of the shaft, the more they contribute to bending stiffnesses. For the rotating shaft with simply supported boundary condition, the first mode is symmetric while the second and third modes are the antisymmetric bending modes. As seen in these figures, rotating speed can cause significant changes to the natural frequencies. However, the effects of rotating speed on the mode shapes are not important [21]. The results show that the higher modes have much more critical rotating speed compared with the lower one. Moreover, the effect of ply angle appears to be more significant for the highermode ones.
Table 1Modal frequency and damping of various laminated tubular circular beams (L/d= 26)
Lamination  Mode  Natural frequency (Hz)  Loss factor (%)  
Ref. [20]  Present  Ref. [17]  Ref. [20]  Present  Ref. [17]  
[0_{2}/90_{2}/45_{2}/45_{2}]_{8}  First flapping  2.4  2.4888  2.5391  1.44  1.4576  1.5277 
Second flapping  15.0  15.3542  15.5312  1.44  1.4576  1.5277.  
First sweeping  2.4  2.4888  2.5391  1.44  1.4576  1.5277  
Second sweeping  15.0  15.3542  15.5312  1.44  1.4576  1.5277  
First torsional  49.2  48.1261  49.4021  1.60  1.6391  1.8203  
Second torsional  148.0  149.8562  152.0235  1.60  1.6391  1.8203  
[45/45]_{8}  First flapping  2.0  2.0704  2.3025  2.45  2.4712  2.6091 
Second flapping  13.0  13.7984  13.8541  2.44  2.4593  2.5125  
First sweeping  2.0  2.0704  2.3025  2.45  2.4712  2.6091  
Second sweeping  13.0  13.7984  13.8541  2.44  2.4593  2.5125  
First torsional  57.0  56.8652  58.9652  1.00  1.1047  1.3943  
Second torsional  172.0  173.5241  175.6352  1.00  1.1047  1.3943 
Table 2Modal frequency and damping of a boxsection beam for various skin laminate configurations (L/a= 14.36, a/b= 5)
Lamination  Mode  Natural frequency (Hz)  Loss factor (%)  
Ref. [20]  Present  Ref. [17]  Ref. [20]  Present  Ref. [17]  
[0_{2}/90_{2}/45_{2}/45_{2}]_{8 }  First flapping  2.4  2.2392  2.6532  1.44  1.4523  1.6134 
Second flapping  14.8  15.5809  16.0583  1.44  1.4523  1.6134  
First sweeping  8.3  8.3368  8.5981  1.44  1.4523  1.6134  
Second sweeping  50.9  49.0480  51.5689  1.45  1.4731  1.6243  
First torsional  46.9  47.3335  48.6523  1.61  1.6279  1.9756  
Second torsional  140.9  142.0065  145.2563  1.61  1.6279  1.9756  
[45/45]_{8}  First flapping  2.0  1.9058  2.2230  2.45  2.4715  2.7932 
Second flapping  12.7  13.4369  13.5246  2.45  2.4715  2.7932  
First sweeping  7.1  7.2860  7.2906  2.45  2.4715  2.7932  
Second sweeping  44.1  43.4169  45.8832  2.41  2.4437  2.6417  
First torsional  54.6  55.2003  56.5263  0.99  1.0752  1.2578  
Second torsional  164.2  165.6008  167.4320  0.99  1.0752  1.2578 
Table 3Mechanical properties of material
$\rho $ (kg/m^{3})  ${E}_{11}$ (GPa)  ${E}_{22}$ (GPa)  ${G}_{12}$ (GPa)  ${G}_{23}$ (GPa)  ${\nu}_{12}$  ${\eta}_{11}$ (%)  ${\eta}_{12}$ (%)  ${\eta}_{13}$ (%)  ${\eta}_{14}$ (%) 
1672  25.8  8.7  3.5  3.5  0.34  0.65  2.34  2.89  2.89 
Fig. 2The natural frequency of a simply supported composite shaft versus rotating speed for different ply angles (L/r= 52, θ/θ8)
a) First mode
b) Second mode
c) Third mode
Fig. 3The damping of a simply supported composite shaft versus rotating speed for different ply angles (L/r= 52, θ/θ8)
a) First mode
b) Second mode
c) Third mode
Fig. 4The natural frequency of a simply supported composite shaft versus rotating speed for different length aspect ratios (60/608)
a) First mode
b) Second mode
c) Third mode
Figs. 3(a)(c) present the first three flexural dampings with respect rotating speed for different ply angles. It is evident that as the rotating speed increases, the dampings of the backward modes decrease and are negative, this means the motion involving backward modes is stable. The dampings involving the forward modes are negative at the lower value of $\mathrm{\Omega}$ and increase with $\mathrm{\Omega}$ increasing, and at a certain value of $\mathrm{\Omega}$ the dampings become zero and then are positive. The value of $\mathrm{\Omega}$ corresponding to zero damping is called as the threshold of instability. Furthermore, it can also be seen that the threshold of instability decreases with ply angle, which is explained by the fact the closer the fiber is oriented to 90°, the higher the internal damping of composite (see specific damping capacity given in Table 3), thus, the lower threshold of instability is. Comparing with the lower mode, the higher modes have larger threshold of instability. This shows that the higher modes are more stable.
Figs. 4 and Figs. 5 show the first three natural frequencies and flexural dampings with respect rotating speed for various length aspect ratios. The results show that as the length aspect ratio increases, the nonrotating natural frequencies decrease and a corresponding decrease in the instability thresholds is also experienced. It can be concluded that for relatively large value of length aspect ratio, the shaft becomes less stability than for low length aspect ratio. A similar conclusion has been arrived by Sino, et al. [14].
As the results reveal, the higher mode is much stronger affected by length aspect ratio, than the lower mode.
Figs. 6 and 7 show the first three natural frequencies and and flexural dampings with respect rotating speed for two typical stacking sequences. The results show that stacking sequence can be used to increase the stability of the rotating composite shaft. From these figures it is evident that the effect of the stacking sequences on the higher modes are significant.
Fig. 5The damping of a simply supported composite shaft versus rotating speed for different length aspect ratios (60/608)
a) First mode
b) Second mode
c) Third mode
Fig. 6The natural frequency of a simply supported composite shaft versus rotating speed for different stacking sequences (L/r= 52)
a) First mode
b) Second mode
c) Third mode
Fig. 7The damping of a simply supported composite shaft versus rotating speed for different stacking sequences (L/r= 52)
a) First mode
b) Second mode
c) Third mode
Fig. 8The natural frequency of a composite shaft versus rotating speed for different boundary conditions (L/r= 52, 60/608)
a) First mode
b) Second mode
c) Third mode
Fig. 9The damping of a composite shaft versus rotating speed for different boundary conditions (L/r= 52, 60/608)
a) First mode
b) Second mode
c) Third mode
Figs. 8 and Figs. 9 show the first three natural frequencies and flexural dampings with respect rotating speed for three different boundary conditions. The results show that boundary condition has a significant effect on the stability of the rotating composite shaft. From Figs. 8, it is obviously that the boundary conditions can exert a great influence over the natural frequency of the higher modes.
Figs. 10 and 11 show the first critical speeds and instability thresholds with respect ply angle for various length aspect ratios, respectively. The results show that the first critical speed and instability threshold decrease with the ply angle and increase with a decrease in length aspect ratio. These trends are consistent with that presented in Figs. 25.
Figs. 12 and 13 show the effect of three different boundary conditions on the first critical speeds and instability thresholds of the shaft, respectively. The results show that the maximum critical speed and instability threshold occur for the clampedclamped (CC) shaft, the minimum ones for the clampedfree (CF) shaft, while the intermediate ones occur for the simply supported (SS) shaft.
Figs. 14 and 15 show the effect of the stacking sequence on the first critical speeds and instability thresholds of the shaft, respectively. The effect of the stacking sequence on the first critical speed and instability threshold is similar to those indicated in Figs. 6 and 7. It can be observed that the angleply ${\left[\pm \theta \right]}_{8}$ seems to exert a larger stabilizing effect, as compared to the stacking sequence ${\left[\theta \right]}_{16}$.
Convergence studies with our model indicate that five mode shape functions are sufficient to obtain the accurate results of the first critical speeds and instability thresholds for simply supported shaft as shown in Table 4. This shows that our model is sufficiently accurate even when only a small number of mode shape functions are used.
Fig. 10The variation of critical speeds with ply angle for different length aspect ratios (θ/θ8, simply supported)
Fig. 11The variation of instability thresholds with ply angle for different length aspect ratios (θ/θ8, simply supported)
Fig. 12The variation of critical speeds with ply angle for different boundary conditions (L/r= 52, θ/θ8)
Fig. 13The variation of instability thresholds with ply angle for different boundary conditions (L/r= 52, θ/θ8)
Fig. 14The variation of critical speeds with ply angle for different stacking sequences (L/r= 52, simply supported)
Fig. 15The variation of instability thresholds with ply angle for different stacking sequences (L/r= 52, simply supported)
Table 4Effect of model number N on first critical speed and instability threshold (L/r= 52, θ/θ8 simply supported)
Ply angle (degree)  First critical speed (rpm)  Instability threshold (rpm)  
$N=$1  $N=$3  $N=$5  $N=$1  $N=$3  $N=$5  
30  44.65  44.72  44.72  45.47  45.62  45.62 
60  40.44  40.64  40.64  42.03  42.35  42.35 
90  36.58  36.67  36.67  38.03  38.23  38.23 
4. Conclusions
A refined variational asymptotic approach and Galerkin’s method was used to model the free vibration equations and to predict the dynamical behavior the rotating composite Timoshenko shafts with internal damping. The present numerical results are compatible to those available in the literature. The critical speeds and instability thresholds are obtained using the present model. From the presented results, the following conclusions can be made:
1) The proper stacking sequence is useful for increasing critical rotating speeds and instability thresholds of composite shaft.
2) The clampedclamped composite shaft has the higher critical rotating speed and threshold of instability and consequently better dynamic stability than the shaft with other boundary conditions
3) The ply angles up to 50°^{}have a significant effect on dynamics of composite shaft and at the ply angle 0°, maximum critical speed and instability threshold are reached. However, for the ply angles in the range 50°90°, their effect appears to be marginal.
4) The instability threshold is always higher than the critical rotating speed.
References

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.

dos Reis H. L. M., Goldman R. B., Verstrate P. H. Thinwalled laminated composite cylindrical tubes; part 3 – critical speed analysis. Journal of Composites Technology and Research, Vol. 9, 1987, p. 5862.

Kim C. D., Bert C. W. Critical speed analysis of laminated composite, hollow drive shafts. Composites Engineering, Vol. 3, 1993, p. 633643.

Singh S. P., Gupta K. Composite shaft rotordynamic analysis using a layerwise theory. Journal of Sound and Vibration, Vol. 191, Issue 5, 1996, p. 739756.

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

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

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 thinwalled beam. Journal of Acoustical Society of America, Vol. 109, Issue 3, 2001, p. 972981.

Rehfield L. W. Design analysis methodology for composite rotor blades. Proceedings of the 7th DoD/NASA Conference on Fibrous Composites in Structural Design, 1985.

Ren Yongsheng, Dai Qiyi, Zhang Xingqi Modeling and dynamic analysis of rotating composite shaft. Journal of Vibroengineering, Vol. 15, Issue 4, 2013, p. 18161832.

Singh S. P., Gupta K. Free damped flexural vibration analysis of composite cylindrical tubes using beam and shell theories. Journal of Sound and Vibration, Vol. 172, Issue 2, 1994, p. 17190.

Kim W., Argento A., Scott R. A. Forced vibration and dynamic stability of a rotating tapered composite Timoshenko shaft: bending motions in endmilling operations. Journal of Sound and Vibration, Vol. 246, 2001, p. 563600.

Mazzei A. J., Scott R. A. Effects of internal viscous damping on the stability of a rotating shaft driven through a universal joint. Journal of Sound and Vibration, Vol. 265, 2003, p. 863885.

Montagnier O., Hochard C. H. Dynamics of a supercritical composite shaft mounted on viscoelastic supports. Journal of Sound and Vibration, Vol. 333, 2014, p. 470484.

Sino R., Baranger T. N., Chatelet E., Jacquet G. Dynamic analysis of a rotating composite shaft. Composites Science and Technology, Vol. 68, 2008, p. 337345.

Ren Yongsheng, Zhang Xingqi, Liuyanghang, Chen Xiulong Vibration and instability of rotating composite thinwalled shafts with internal damping. Shock and Vibration, 2014, p. 123271110.

Berdichevsky V. L., Armanios E., Badir A. M. Theory of anisotropic thinwalled closedcrosssection beams. Composites Engineering, Vol. 2, Issues 57, 1992, p. 411432.

Ren Yongsheng, Du Xianghong, Sun Shuangshuang, Teng Xiangmeng Structure damping of thinwalled composite onecell beams. Journal of Vibration and Shock, Vol. 31, Issue 3, 2012, p. 141152, (in Chinese).

Dharmarajan S., Mccutchen H. Shear coefficients for orthotropic beams. Journal of Composite Materials,Vol. 7, 1973, p. 530535.

Smith E. C., Chopra I. Formulation and evaluation of an analytical model for composite boxbeams. Journal of American Helicopter Society, Vol. 36, Issue 3, 1991, p. 2325.

Saravanos D. A., Varelis D., Plagianakos T. S., et al. A shear beam finite element for the damping analysis of tubular laminated composite beams. Journal of Sound and Vibration, Vol. 291, 2006, p. 802823.

Banerjee J. R., Su H. Development of a dynamic stiffness matrix for free vibration analysis of spinning beams. Computer and Structures, Vol. 82, 2004, p. 21892197.
About this article
The research is funded by the National Natural Science Foundation of China (Grant No. 11272190), Shandong Provincial Natural Science Foundation of China (Grant No. ZR2011EEM031).