Abstract
Structural modeling and dynamical analysis of rotating composite shaft are conducted in this paper. A thinwalled 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 thinwalled 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 timeintegration method.
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 thickshell theories of firstorder approximation to derive the motion equations of the rotating composite thinwalled 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 firstorder 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 thinwalled beam theory [8] that presented a composite thinwalled 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 thinwalled shaft model is proposed by means of the composite thinwalled 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 rotorcomposite 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 timeintegration 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 $\Omega $ is shown in Fig. 1. To describe the motion of the shaft the following coordinate systems are considered: (1) $\left(X,Y,Z\right)$ is inertial reference system whose origin $O$ is located in the geometric center, (2) $\left(x,y,z\right)$ is rotating reference system with the common origin $O\text{.}$$\left(I,J,K\right)\mathrm{}$and $\left(i,j,k\right)$ denote the unit vectors of the reference systems $\left(X,Y,Z\right)$ and $\left(x,y,z\right)\text{,}$ respectively. In addition a local coordinate system $\left(\xi ,s,x\right)$ is used, where $\xi $ 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 thinwalled elastic cylinder, satisfying $d\ll L$, $h\ll d$, $h\ll r$, 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 $\theta \left(y\right)=\theta \left(y\right)$, $\theta \left(z\right)=\theta \left(z\right)$, is considered (Fig. 1).
Fig. 1Composite thinwalled shaft of a circular cross section
It has been shown [1012] that VAM is an asymptotically correct theory which can be used effectively for the analysis of tubular composite thinwalled 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:
${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),$
in which:
${\theta}_{z}\left(x,t\right)={U}_{3}^{\text{'}}\left(x,t\right)2{\gamma}_{yx},$
where ${U}_{1}\left(x,t\right),{U}_{2}\left(x,t\right),{U}_{3}\left(x,t\right)$ denote the rigidbody translations along the $x$, $y$ and $z$axis, while $\varphi \left(x,t\right),{\theta}_{y}\left(x,t\right),{\theta}_{z}\left(x,t\right)$ denote the twist about the $z$axis and rotations about the $x$ and $y$axis respectively. ${\gamma}_{zx}$ and ${\gamma}_{yx}$ denote the transverse shear strains in the planes $xz$ and $xy$ respectively.
From the classical VAM the warping function $g\left(s,x,t\right)$ is modified as:
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)$ 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 straindisplacement relations [9], while referring to the [13], the strain of the composite shaft can be written as:
$2{\gamma}_{xs}=\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},$
where ${r}_{n}$ denotes the project of the position vector $\mathit{r}$ of an arbitrary point on the crosssection of the deformed shaft in the normal direction:
The position vector $\mathit{r}$ can be expressed as:
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:
The strain energy of the composite shaft ${U}_{s}$ can be written as:
where ${\sigma}_{xx}\text{,}\mathrm{}$${\sigma}_{xs}$ and ${\sigma}_{x\xi}$ represent the engineering stresses associated with the engineering strains$\mathrm{}{\epsilon}_{xx}$, ${\epsilon}_{xs}$ and ${\epsilon}_{x\xi}$.
Taking variation for the above strain energy expression one obtains:
Stress resultants ${F}_{x}$, ${Q}_{y}$, ${Q}_{z}$ and stress couples ${M}_{x}$, ${M}_{y}$, ${M}_{z}$can be defined as:
${M}_{x}=\oint {N}_{xs}{r}_{n}ds,{M}_{y}=\oint {N}_{xx}zds,{M}_{z}=\oint {N}_{xx}yds,$
${Q}_{y}=\oint \left({N}_{xs}\frac{dy}{ds}+{N}_{x\xi}\frac{dz}{ds}\right)ds,{Q}_{Z}=\oint \left({N}_{xs}\frac{dz}{ds}{N}_{x\xi}\frac{dy}{ds}\right)ds,$
where ${N}_{xx}$, ${N}_{xs}$ and ${N}_{x\xi}$ are shell stress resultants defined according to the following expressions:
${N}_{x\xi}\mathrm{}\mathrm{}=D\left(s\right){\gamma}_{x\xi},$
where:
$C\left(s\right)=4\left[{A}_{66}\frac{{{(A}_{26})}^{2}}{{A}_{22}}\right],D\left(s\right)=\left({A}_{44}\frac{{A}_{45}^{2}}{{A}_{55}}\right),$
${A}_{ij}=\sum _{k=1}^{N}{\overline{Q}}_{ij}^{\left(k\right)}\left({Z}_{k}{Z}_{k1}\right)\mathrm{},\left(i,j=\mathrm{1,2},6;i,j=\mathrm{4,5}\right).$
Parameters $A\left(s\right)$ and $B\left(s\right)$ denote the reduced axial and coupling stiffnesses respectively, while $C\left(s\right)$ and $D\left(s\right)$denote the reduced shear stiffnesses. ${A}_{ij}$ denote the inplane stiffnesses components and ${\stackrel{}{Q}}_{ij}^{\left(k\right)}$ 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:
${M}_{x}={k}_{12}{U}_{1}^{\text{'}}+{k}_{22}{\varphi}^{\text{'}}+{k}_{23}{\theta}_{z}^{\text{'}}+{k}_{24}{\theta}_{y}^{\text{'}}+{k}_{25}\left({U}_{2}^{\text{'}}{\theta}_{y}\right)+{k}_{26}\left({U}_{3}^{\text{'}}{\theta}_{z}\right),$
${M}_{y}={k}_{13}{U}_{1}^{\text{'}}+{k}_{23}{\varphi}^{\text{'}}+{k}_{33}{\theta}_{z}^{\text{'}}+{k}_{34}{\theta}_{y}^{\text{'}}+{k}_{35}\left({U}_{2}^{\text{'}}{\theta}_{y}\right)+{k}_{36}\left({U}_{3}^{\text{'}}{\theta}_{z}\right),$
${M}_{z}={k}_{14}{U}_{1}^{\text{'}}+{k}_{24}{\varphi}^{\text{'}}+{k}_{34}{\theta}_{z}^{\text{'}}+{k}_{44}{\theta}_{y}^{\text{'}}+{k}_{45}\left({U}_{2}^{\text{'}}{\theta}_{y}\right)+{k}_{46}\left({U}_{3}^{\text{'}}{\theta}_{z}\right),$
${Q}_{y}={k}_{51}{U}_{1}^{\text{'}}+{k}_{52}{\varphi}^{\text{'}}+{k}_{53}{\theta}_{z}^{\text{'}}+{k}_{54}{\theta}_{y}^{\text{'}}+{k}_{55}\left({U}_{2}^{\text{'}}{\theta}_{y}\right)+{k}_{56}\left({U}_{3}^{\text{'}}{\theta}_{z}\right),$
${Q}_{z}={k}_{61}{U}_{1}^{\text{'}}+{k}_{62}{\varphi}^{\text{'}}+{k}_{63}{\theta}_{z}^{\text{'}}+{k}_{64}{\theta}_{y}^{\text{'}}+{k}_{65}\left({U}_{2}^{\text{'}}{\theta}_{y}\right)+{k}_{66}\left({U}_{3}^{\text{'}}{\theta}_{z}\right),$
where ${k}_{ij}\mathrm{}(i,j=1,\dots ,6)$ are the stiffness coefficients of the composite shaft, which can be expressed in terms of the cross sectional geometry and material properties as:
${k}_{12}=\left[{\oint}_{\Gamma}^{}\left(\frac{B}{C}\right)ds/{\oint}_{\Gamma}^{}\left(\frac{1}{C}\right)ds\right]{A}_{e},$
${k}_{13}={\oint}_{\Gamma}^{}\left(A\frac{{B}^{2}}{C}\right)zds\left\{\left[{\oint}_{\Gamma}^{}\left(\frac{B}{C}\right)ds{\oint}_{\Gamma}^{}\left(\frac{B}{C}\right)zds\right]/{\oint}_{\Gamma}^{}\left(\frac{1}{C}\right)ds\right\},$
${k}_{14}={\oint}_{\Gamma}^{}\left(A\frac{{B}^{2}}{C}\right)yds\left\{\left[{\oint}_{\Gamma}^{}\left(\frac{B}{C}\right)ds{\oint}_{\Gamma}^{}\left(\frac{B}{C}\right)yds\right]/{\oint}_{\Gamma}^{}\left(\frac{1}{C}\right)ds\right\},$
${k}_{22}=\left[1/{\oint}_{\Gamma}^{}\left(\frac{1}{C}\right)ds\right]{A}_{e}^{2},$
${k}_{23}=\left[{\oint}_{\Gamma}^{}\left(\frac{B}{C}\right)zds/{\oint}_{\Gamma}^{}\left(\frac{1}{C}\right)ds\right]{A}_{e},$
${k}_{24}=\left[{\oint}_{\Gamma}^{}\left(\frac{B}{C}\right)yds/{\oint}_{\Gamma}^{}\left(\frac{1}{C}\right)ds\right]{A}_{e},$
${k}_{33}={\oint}_{\Gamma}^{}\left(A\frac{{B}^{2}}{C}\right){z}^{2}ds+\left\{{\left[{\oint}_{\Gamma}^{}\left(\frac{B}{C}\right)zds\right]}^{2}/{\oint}_{\Gamma}^{}\left(\frac{1}{C}\right)ds\right\},$
${k}_{34}={\oint}_{\Gamma}^{}\left(A\frac{{B}^{2}}{C}\right)yzds+\left\{\left[{\oint}_{\Gamma}^{}\left(\frac{B}{C}\right)yds{\oint}_{\Gamma}^{}\left(\frac{B}{C}\right)zds\right]/{\oint}_{\Gamma}^{}\left(\frac{1}{C}\right)ds\right\},$
${k}_{44}={\oint}_{\Gamma}^{}\left(A\frac{{B}^{2}}{C}\right){y}^{2}ds+\left\{{\left[{\oint}_{\Gamma}^{}\left(\frac{B}{C}\right)yds\right]}^{2}/{\oint}_{\Gamma}^{}\left(\frac{1}{C}\right)ds\right\},$
${k}_{15}=\frac{1}{2}{\oint}_{\Gamma}^{}B\frac{dy}{ds}ds,{k}_{16}=\frac{1}{2}{\oint}_{\Gamma}^{}B\frac{dz}{ds}ds,{k}_{25}=\frac{1}{4}{\oint}_{\Gamma}^{}{r}_{n}C\frac{dy}{ds}ds,$
${k}_{26}=\frac{1}{4}{\oint}_{\Gamma}^{}{r}_{n}C\frac{dz}{ds}ds,{k}_{35}=\frac{1}{2}{\oint}_{\Gamma}^{}Bz\frac{dy}{ds}ds,{k}_{36}=\frac{1}{2}{\oint}_{\Gamma}^{}Bz\frac{dz}{ds}ds,$
${k}_{45}=\frac{1}{2}{\oint}_{\Gamma}^{}By\frac{dy}{ds}ds,{k}_{46}=\frac{1}{2}{\oint}_{\Gamma}^{}By\frac{dz}{ds}ds,$
${k}_{55}={\oint}_{\Gamma}^{}\left[\frac{1}{4}{C\left(\frac{dy}{ds}\right)}^{2}+{D\left(\frac{dz}{ds}\right)}^{2}\right]ds,{k}_{56}={\oint}_{\Gamma}^{}\left(\frac{1}{4}CD\right)\frac{dy}{ds}\frac{dz}{ds}ds,$
${k}_{66}={\oint}_{\Gamma}^{}\left[\frac{1}{4}{C\left(\frac{dz}{ds}\right)}^{2}+{D\left(\frac{dy}{ds}\right)}^{2}\right]ds.$
The coefficients ${k}_{ij}(i,j=\mathrm{1,2},\dots ,4)$ represent the stiffnesses, as described in [9, 12] for the case of a nonshearable thinwalled beam, while the shear contribution to the stiffness coefficients is represented by the new terms ${k}_{ij}\left(i=\mathrm{1,2},\dots ,6;j=\mathrm{5,6}\right)$, ${k}_{ji}\mathrm{}(i=\mathrm{1,2},\dots ,4;j=\mathrm{5,6})$.
The kinetic energy of the composite shaft ${T}_{s}$ can be written as:
where $\rho $ 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 ${T}_{s}$ yields:
where:
${I}_{2}={b}_{1}\left({\ddot{U}}_{2}2\Omega {\dot{U}}_{3}{\Omega}^{2}{U}_{2}\right){b}_{2}\left(2\Omega \dot{\varphi}+{\Omega}^{2}\right){b}_{3}(\ddot{\varphi}{\Omega}^{2}\varphi ),$
${I}_{3}={b}_{1}\left({\ddot{U}}_{3}+2\Omega {\dot{U}}_{2}{\Omega}^{2}{U}_{3}\right)+{b}_{2}\left(\ddot{\varphi}{\Omega}^{2}\varphi \right){b}_{3}(2\Omega \dot{\varphi}+{\Omega}^{2}),$
${I}_{4}={b}_{2}{\ddot{U}}_{1}{b}_{4}{\ddot{\theta}}_{y}{b}_{6}{\ddot{\theta}}_{z}\mathrm{},$
${I}_{5}={b}_{3}{\ddot{U}}_{3}{b}_{6}{\ddot{\theta}}_{y}{b}_{5}{\ddot{\theta}}_{z}\mathrm{},$
${I}_{6}={b}_{2}\left({\ddot{U}}_{3}+2\Omega {\dot{U}}_{2}{\Omega}^{2}{U}_{3}\right)+\left({b}_{4}+{b}_{5}\right)\left(\ddot{\varphi}{\Omega}^{2}\varphi \right){b}_{3}({\ddot{U}}_{2}2\Omega {\dot{U}}_{3}{\Omega}^{2}{U}_{2}),$
${b}_{1}={\iint}_{A}^{}\rho dA,{b}_{2}={\iint}_{A}^{}\rho ydA,{b}_{3}={\iint}_{A}^{}\rho zdA,$
${b}_{4}={\iint}_{A}^{}\rho {y}^{2}dA,{b}_{5}={\iint}_{A}^{}\rho {z}^{2}dA,{b}_{6}={\iint}_{A}^{}\rho yzdA.$
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:
in which:
${I}_{2l}^{D}={b}_{1l}^{D}({\ddot{U}}_{2}2\Omega {\dot{U}}_{3}{\Omega}^{2}{U}_{2}){b}_{2l}^{D}\left(2\Omega \dot{\varphi}+{\Omega}^{2}\right){b}_{3l}^{D}(\ddot{\varphi}{\Omega}^{2}\varphi ),$
${I}_{3l}^{D}={b}_{1l}^{D}\left({\ddot{U}}_{3}+2\Omega {\dot{U}}_{2}{\Omega}^{2}{U}_{3}\right)+{b}_{2l}^{D}\left(\ddot{\varphi}{\Omega}^{2}\varphi \right){b}_{3l}^{D}(2\Omega \dot{\varphi}+{\Omega}^{2}),$
${I}_{4l}^{D}={b}_{2l}^{D}{\ddot{U}}_{1}{b}_{4l}^{D}{\ddot{\theta}}_{y}{b}_{6l}^{D}{\ddot{\theta}}_{Z},$
${I}_{5l}^{D}={b}_{3l}^{D}{\ddot{U}}_{3}{b}_{6l}^{D}{\ddot{\theta}}_{y}{b}_{5l}^{D}{\ddot{\theta}}_{Z},$
${I}_{6l}^{D}={b}_{2l}^{D}\left({\ddot{U}}_{3}+2\Omega {\dot{U}}_{2}{\Omega}^{2}{U}_{3}\right)+\left({b}_{4l}^{D}+{b}_{5l}^{D}\right)\left(\ddot{\varphi}{\Omega}^{2}\varphi \right){b}_{3l}^{D}\left({\ddot{U}}_{2}2\Omega {\dot{U}}_{3}{\Omega}^{2}{U}_{2}\right),$
${b}_{1l}^{D}={\iint}_{A}^{}{\rho}_{l}^{D}dA,{b}_{2l}^{D}={\iint}_{A}^{}{\rho}_{l}^{D}ydA,\mathrm{}{b}_{3l}^{D}={\iint}_{A}^{}{\rho}_{l}^{D}zdA,$
${b}_{4l}^{D}={\iint}_{A}^{}{\rho}_{l}^{D}{y}^{2}dA,{b}_{5l}^{D}={\iint}_{A}^{}{\rho}_{l}^{D}{z}^{2}dA,{b}_{6l}^{D}={\iint}_{A}^{}{\rho}_{l}^{D}yzdA,$
where ${\rho}_{l}^{D}$ denotes the mass density of the disk $l$. The symbol $\Delta (x{x}_{Dl})$ denotes Dirac delta function, ${N}_{D}$ the number of the disks mounted on the shaft and ${x}_{Dl}$ the location of the $l$th 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:
where ${n}_{b}$ denotes the number of the bearings, ${x}_{bl}$ the location, ${K}^{bl}$ and ${C}^{bl}$ the equivalent stiffness and damping coefficients of the $l$th bearings.
The virtual work done by the centrifugal force of rigid disks can be given by:
$\delta {W}_{e}={\int}_{0}^{L}\left({p}_{x}^{u}\delta {U}_{1}+{p}_{y}^{u}\delta {U}_{2}+{p}_{z}^{u}\delta {U}_{3}+{\Gamma}_{x}^{u}\delta \varphi +{\Gamma}_{y}^{u}\delta {\theta}_{y}+{\Gamma}_{z}^{u}\delta {\theta}_{z}\right)dx,$
where ${p}_{x}^{u}$, ${p}_{y}^{u}$, ${p}_{z}^{u}$ denote centrifugal force per unit length, ${\Gamma}_{x}^{u}$, ${\Gamma}_{y}^{u}$, ${\Gamma}_{z}^{u}$ denote centrifugal moment per unit length. They are given by:
${\Gamma}_{x}^{u}=0,{\Gamma}_{y}^{u}=0,{\Gamma}_{z}^{u}=0,$
where the mass unbalance associated with disks is assumed as the concentrated masses ${b}_{1l}^{D}={\iint}_{A}^{}{\rho}_{l}^{D}dA$ at points with small distances of eccentricity ${e}_{l}$.
2.4. Motion equations of systems
In view of Eqs. (9), (16), (28) and (30), and employing the following variational principle:
one obtains:
${Q}_{y}^{\mathrm{\text{'}}}+{I}_{2}+{I}_{2}^{D}+{F}_{{U}_{2}}^{b}={p}_{y}^{u},$
${Q}_{z}^{\mathrm{\text{'}}}+{I}_{3}+{I}_{3}^{D}+{F}_{{U}_{3}}^{b}={p}_{z}^{u},$
${M}_{x}^{\mathrm{\text{'}}}+{I}_{6}+{I}_{6}^{D}={\Gamma}_{x}^{u},$
${M}_{z}^{\mathrm{\text{'}}}+{I}_{4}+{I}_{4}^{D}={\Gamma}_{y}^{u},$
${M}_{y}^{\mathrm{\text{'}}}+{I}_{5}+{I}_{5}^{D}={\Gamma}_{z}^{u},$
where:
${F}_{{U}_{3}}^{b}=\sum _{l=1}^{{n}_{b}}\left({K}_{zz}^{bl}{U}_{3}+{K}_{zy}^{bl}{U}_{2}+{K}_{zz}^{bl}{\dot{U}}_{3}+{K}_{zy}^{bl}{\dot{U}}_{2}\right)\mathrm{\Delta}(x{x}_{bl}).$
The boundary conditions at $x=0$, $L$ are:
${M}_{x}={\stackrel{}{M}}_{x}\text{}\mathrm{o}\mathrm{r}\text{}{\theta}_{y}={\stackrel{}{\theta}}_{y},{M}_{y}={\stackrel{}{M}}_{y}\text{}\mathrm{o}\mathrm{r}\text{}{\theta}_{z}={\stackrel{}{\theta}}_{z},{M}_{z}={\stackrel{}{M}}_{z}\text{}\mathrm{o}\mathrm{r}\text{}\varphi =\stackrel{}{\varphi}.$
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 bendingtransverse shear and extensiontwist motions.
The motion equations of bendingtransverse shear coupling:
$+\sum _{l=1}^{{N}_{D}}\left[{b}_{1l}^{D}\left({\ddot{U}}_{2}2\Omega {\dot{U}}_{3}{\Omega}^{2}{U}_{2}\right)\right]\left(x{x}_{Dl}\right)+{F}_{{U}_{2}}^{b}={p}_{y}^{u},$
${k}_{46}{\theta}_{y}^{\text{'}\text{'}}{k}_{66}\left({U}_{3}^{\text{'}\text{'}}{\theta}_{z}^{\text{'}}\right)+{b}_{1}\left({\ddot{U}}_{3}+2\Omega {\dot{U}}_{2}{\Omega}^{2}{U}_{3}\right)+\sum _{l=1}^{{N}_{D}}\left[{b}_{1l}^{D}\left({\ddot{U}}_{3}+2\Omega {\dot{U}}_{2}{\Omega}^{2}{U}_{3}\right)\right]\left(x{x}_{Dl}\right)+{F}_{{U}_{3}}^{b}={p}_{z}^{u},$
${k}_{44}{\theta}_{y}^{\text{'}\text{'}}{k}_{46}\left({U}_{3}^{\text{'}\text{'}}{\theta}_{z}^{\text{'}}\right){k}_{35}{\theta}_{z}^{\text{'}}{k}_{55}\left({U}_{2}^{\text{'}}{\theta}_{y}^{}\right){b}_{4}{\ddot{\theta}}_{y}+\sum _{l=1}^{{N}_{D}}\left({b}_{4l}^{D}{\ddot{\theta}}_{y}^{}\right)\Delta \left(x{x}_{Dl}\right)=0,$
${k}_{33}{\theta}_{z}^{\text{'}\text{'}}{k}_{35}\left({U}_{2}^{\text{'}\text{'}}{\theta}_{y}^{\text{'}}\right){k}_{46}{\theta}_{y}^{\text{'}}{k}_{66}\left({U}_{3}^{\text{'}}{\theta}_{z}^{}\right){b}_{5}{\ddot{\theta}}_{z}+\sum _{l=1}^{{N}_{D}}\left({b}_{5l}^{D}{\ddot{\theta}}_{z}^{}\right)\Delta \left(x{x}_{Dl}\right)=0.$
The motion equations of extensiontwist coupling:
${k}_{12}{U}_{1}^{\text{'}\text{'}}{k}_{22}{\varphi}^{\text{'}\text{'}}+\left({b}_{4}+{b}_{5}\right)\left(\ddot{\varphi}{\Omega}^{2}\varphi \right)=0.$
By using the coordinate transformation Eq. (26) can be transformed to the inertia reference frame $\left(X,Y,Z\right)$. 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 bendingtransverse shear coupling.
2.5. Approximate solution method
In order to find the approximate solution of the rotating composite shaft, the quantities ${U}_{2}\left(x,t\right)$, ${U}_{3}\left(x,t\right)$, ${\theta}_{y}(x,t)$ and ${\theta}_{z}(x,t)$ are assumed in the form:
${\theta}_{y}\left(x,t\right)=\sum _{j=1}^{N}{\Theta}_{j}\left(t\right){\psi}_{j}\left(x\right),\mathrm{}\mathrm{}{\theta}_{z}\left(x,t\right)=\sum _{j=1}^{N}{\Theta}_{j}\left(t\right){\psi}_{j}\left(x\right),$
where ${\alpha}_{j}\left(x\right)$ and ${\psi}_{j}\left(x\right)$ denote mode shape functions which fulfill all the boundary conditions of the composite shaft, ${U}_{j}\left(t\right)$ and ${\Theta}_{j}\left(t\right)$ 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:
where $\left[M\right]$ denotes the mass matrix, $\left[C\right]$ the gyroscopic matrix, $\left[{K}_{1}\right]$ the elastic stiffness matrix, $\left[{K}_{\Omega}\right]$ stiffness matrix which depends on the rotational speed $\Omega $, $\left\{f\right\}$ the generalized force vector. $\left\{\ddot{X}\right\}\text{,}$$\left\{\dot{X}\right\}$ and $\left\{X\right\}$ denote the generalized acceleration, velocity and displacement vectors respectively.
These matrices can be written as follows:
where:
${F}_{ij}^{Dl}={\int}_{0}^{L}{\psi}_{i}{\psi}_{j}\Delta \left(x{x}_{Dl}\right)dx,{E}_{ij}^{bl}={\int}_{0}^{L}{\alpha}_{i}{\alpha}_{j}\Delta \left(x{x}_{bl}\right)dx,{I}_{ij}={\int}_{0}^{L}{\alpha}_{i}{\alpha}_{j}^{\text{'}\text{'}}dx,$
${J}_{ij}={\int}_{0}^{L}{\alpha}_{i}{\psi}_{j}^{\text{'}}dx,{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}_{i}^{Dl}={\int}_{0}^{L}{\alpha}_{i}\Delta \left(x{x}_{Dl}\right)dx,(i,j=\mathrm{1,2},\dots ,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 1Effect of model number N on natural frequencies (Ω*= 0 and θ= 30°)
$N$  ${\omega}_{1}^{*}$  ${\omega}_{2}^{*}$  ${\omega}_{3}^{*}$  ${\omega}_{4}^{*}$  ${\omega}_{5}^{*}$  ${\omega}_{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 2Comparison of the natural frequencies of a cantilever composite shaft without shear deformation
${\Omega}^{*}$  ${\omega}_{1}^{*}$  ${\omega}_{2}^{*}$  ${\omega}_{3}^{*}$  ${\omega}_{4}^{*}$  ${\omega}_{5}^{*}$  ${\omega}_{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 ${\omega}^{*}=\omega /{\omega}_{0}$, ${\Omega}^{*}=\Omega /{\omega}_{0}$, where the normalizing factor ${\omega}_{0}=$138.85 rad/s is the fundamental frequency of the nonrotating shaft with $\theta =0\xb0$.
From this figure it clearly appears that when ${\Omega}^{*}=$0 a single zerospeed mode natural frequency is obtained. As the rotating speed ${\Omega}^{*}$ 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. 2The 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 midpoint 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 thinwalled shaft and consistent with assumption as used in the present model. The shaft is made by boronepoxy composite material with parameters as follows: ${E}_{11}=\text{2.11 GPa,}$${E}_{22}=\text{24.1}\text{}\text{GPa}\text{,}$${G}_{12}={G}_{13}={G}_{23}=$6.9 GPa, ${\nu}_{12}=$0.36, $\rho =$1967.0 kg/m^{3} [3]. The stacking sequence of the composite shaft is ${[\pm \theta ]}_{\text{5}}$. The properties of the disk and the two bearings are given in Table 3.
Table 3The properties of disk and bearings [5]
Disk  Bearings  
${b}_{1l}^{D}$ (kg)  2.4364   
$e$ (10^{5} m)  5.0000   
${b}_{4l}^{D},{b}_{5l}^{D}$ (kg m^{2})  0.1901   
${K}_{yy},{K}_{zz}$ (10^{7} N/m)    1.75 
${C}_{yy},{C}_{zz}$ (10^{7} Ns/m)    5.00 
Fig. 3The composite shaft system with bearings and disk
Fig. 4The Campbell diagram of a composite shaft system (θ= 0°)
Fig. 5The Campbell diagram of a composite shaft system (θ= 60°)
Fig. 6The 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. 46 present the Campbell diagrams which display the variation of the first three frequencies with respect to rotating speed for various fiber plyangles. 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 plyangles $\theta =\mathrm{}$0°, 60°, 90° are ${\Omega}_{cr}=$1414 rpm, 2477 rpm and 4023 rpm respectively. The results show that nonrotating natural frequencies and the critical rotating speeds increase as fiber plyangle increases. The reason is that the bending siffnesses ${k}_{33}$ and ${k}_{44}$ increase when fiber plyangle increases.
Fig. 7The Campbell diagram of a composite shaft system (— with shear deformation;    without shear deformation)
Fig. 8The first two critical speeds of a composite shaft system versus ratio of length over radius
Fig. 9The 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 ${\Omega}_{cr}=$ 2477 rpm and 3081 rpm for the first critical rotating speed and ${\Omega}_{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 thinwalled beams overestimate 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 plyangles respectively. The results show that critical speeds decrease as $L/r$ increases, whereas they increase as $r/h$ increases.
Fig. 10The first two natural frequencies of a composite shaft system versus ratio of length over radius (Ω= 20000 rpm)
Fig. 11The 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 plyangles respectively. Because of the influences of rotating speed ($\Omega =$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 highermode ones than the lowermode ones.
Fig. 12Trajectory 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)
b)
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 $\Omega ={\Omega}_{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. 13The time response of the displacement at the disk center: (a) U2(L/2,t); (b) U3(L/2,t)
a)
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 thinwalled composite beam theory referred to as VAM. Numerical simulations of the effects of various parameters including fiber plyangle, 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.
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 III – Critical speed analysis. J. Compos. Technol. Res., 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. 231248.

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. Proc. Seventh DoD/NASA Conf. on Fibrous Composites in Structural Design, AFWALTR853094, 1985.

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

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

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

Armanios E. A., Badir A. M. Free vibration analysis of anisotropic thinwalled closedsection beams. AIAA J., Vol. 33, Issue 10, 1995, p. 19051910.

Song O., Librescu L. Structural modeling and free vibration analysis of rotating composite thinwalled beams. Journal of the American Helicopter Society, 1997, p. 358369.

Banerjee J. R., Su H. Development of a dynamic stiffness matrix for free vibration analysis of spinning beams. Comput. Struct., Vol. 82, 2004, p. 21892197.

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.
About this article
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).