Generalized Differential Quadrature Method for Free Vibration Analysis of a Rotating Composite Thin-Walled Shaft

Abstract

A refined variational asymptotic method (VAM) and Hamilton’s principle were used to establish the free vibration differential equations of a rotating composite thin-walled shaft with circumferential uniform stiffness (CUS) configuration. The generalized differential quadrature method (GDQM) was adopted to discretize and solve the governing equations. The accuracy and efficiency of the GDQM were validated in analyzing the frequency of a rotating composite shaft. Compared to the available results in literature, the computational results by the GDQM are accurate. In addition, effects of boundary conditions, rotating speed, ply angle, ratio of radius over thickness, and ratio of length over radius on the frequency characteristics were also investigated.

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

1. Introduction

A rotating shaft, as an important structure, is commonly used in many applications. Composite materials have been widely applied in shaft structural designs owing to their advantages of strength, specific modulus, corrosion resistance, wear resistance, fatigue life, light weight, high specific stiffness, and design ability.

In the past decades, great effort has been made in developing more accurate and appropriate models and methods to analyze dynamic characteristics of shaft structures. Zinberg and Symonds [1] determined the critical speeds of rotating anisotropic cylindrical shafts simply supported at the ends using equivalent modulus beam theory (EMBT). Dos Reis et al. [2] evaluated the rotordynamic performance of shaft configuration of Zinberg and Symmonds [1] using the finite element method. Kim and Bert [3] derived the motion equations of the rotating composite thin-walled shafts by means of the thin- and thick-shell theories of first-order approximation. They investigated the critical speeds of composite shafts including the effect of bending-twisting coupling. Bert and Kim [4] studied the dynamic instability of drive shafts subjected to fluctuating torque and/or rotational speed by means of various thin shell theories. The theories included the combined effects of torsion and rotation. Singh and Gupta [5] developed two composite spinning shaft models by using EMBT and layerwise beam theory (LBT), respectively. It was indicated that a discrepancy existed between the critical speeds obtained from the two theories for the unsymmetric laminated composite shaft. Chen and Peng [6, 7] analyzed the stability of a rotating shaft subjected to axial periodic forces using the finite element method. In their study, each node contains 5 degrees of freedom. Based on a first-order shear deformable beam theory, Chang et al. [8] presented a simple spinning composite shaft model. They found the approximate solution of the system using the finite element model with three-node one-dimensional line elements that have 6 degrees of freedom in each node. Lagrangian interpolation functions were used to approximate the displacement fields of the shaft. Chang et al. [9] studied the vibration behaviors of the rotating composite shafts containing randomly oriented reinforcements. The finite element model of the composite shaft [8] was extended to the case that contains the fiber inclusions by taking effective elastic moduli into account. Gubran and Gupta [10] modified EMBT with shear deformation, rotary inertia, and gyroscopic effects and analysed the natural frequencies of composite tubular shafts. Song et al. [11] investigated the vibrational and stability behavior of spinning circular shafts modeled as thin-walled beams considering the effects of transverse shear, rotatory inertias, the axial compressive load, and boundary conditions. Sino et al. [12] developed a simplified homogenized beam theory (SHBT) to analyse the sensitivity of the frequencies and instability thresholds regarding shear effect, stacking order, and fiber orientation. Boukhalfa et al. [13] performed free vibration analysis of rotating composite shafts using the p-version hierarchical finite element method with trigonometric shape functions. Boukhalfa and Hadjoui [14] also studied the vibratory behavior of rotating composite shafts. The study incorporated the rotary inertia, gyroscopic effects, transverse shear deformation, and the coupling effect due to the lamination of composite layers. Alwan et al. [15] analyzed solid shafts and composite tube shafts using ANSYS. Boukhalfa [16] considered the dynamic behavior of the spinning Functionally Graded Material (FGM) shaft on rigid bearings. A p-version hierarchical finite element was employed to define the model. Considering the effects of stacking sequences and shear-normal coupling on rotating composite shafts, Arab et al. [17] developed a finite element based on Equivalent Single Layer Theory (ESLT) to model the rotating composite shaft using the Timoshenko beam theory. More recently, Layerwise Shaft Theory (LST) is developed based on shaft finite element theory by Arab et al. [18]. They investigated the dynamic analysis of rotating laminated shafts including t the influences of stacking sequence, fiber orientation, and shear-normal coupling.

Many studies have been done on the vibration analysis of composite shafts using Galerkin methods. Kim et al. [19] investigated the free vibration of a rotating tapered composite Timoshenko shaft. The spatial solutions to the equations of motion were obtained by the general Galerkin method. Oh et al. [20] studied vibration and instability of functionally graded cylindrical spinning thin-walled beams also using the Galerkin method. Na et al. [21] evaluated the vibration and stability of a cylindrical shaft modeled as a tapered thin-walled composite beam and adopted the extended Galerkin method to solve the eigenvalue problem. Ren et al. [22] conducted structural modeling and dynamic analysis of a rotating composite shaft and also used the Galerkin method to discretize and solve the governing equations. Furthermore, Ren et al. [23] studied the dynamic behavior of an internally damped rotating composite shaft considered as an anisotropic thin-walled Timoshenko beam. The Galerkin method was also employed to discretize and solve the equations of motion.

The generalized differential quadrature method (GDQM) is a well-organized numerical technique originating from the differential quadrature method (DQM) proposed by Bellman and Casti [24]. The DQM has been promoted as a potential alternative to conventional numerical solution techniques such as the Galerkin method and finite element method.

So far, no publications have reported on using the GDQM to conduct the free vibration analysis of rotating composite thin-walled shaft. In this paper, the equations of motion of the rotating shaft were established by a refined variational asymptotic method and Hamilton’s principle. The motion governing equations are partial differential equations. With the aid of the GDQM, the partial differential governing equations were transformed approximately into a set of linear algebraic governing equations. Imposing the given boundary conditions, the numerical eigenvalue equations for the free vibration of the shaft were derived and then solved. The dynamic characteristics of the rotating shaft were then analyzed. The computational results using the GDQM were compared to those available in literature, and it was validated that the GDQM is accurate and efficient for the frequency analysis of a rotating composite shaft. Effects of boundary conditions, rotating speed, ply angle, ratio of radius over thickness, and ratio of length over radius on the frequency characteristics were also investigated.

2. Model of Rotating Composite Thin-Walled Shaft

The model of the composite thin-walled rotating shaft is shown in Figure 1, where $L$, $r$, and h denote the length, the radius of curvature of the middle surface, and the thickness, respectively, satisfying $r\ll L$ and $h\ll r$. The following coordinate systems are defined to describe the motion of the rotating shaft: $(1)$  $(X,Y,Z)$ represents the inertial coordinate system whose origin is at the geometric center, and the corresponding unit vector is (I, J,K); $(2)$ ($x,y,z)$ is a rotating reference system and the corresponding unit vector is (i, j, k); $(3)$  $(x,s,\xi )$ is a local coordinate system, where s and $\xi$ represent the directions, tangent, and normal to the middle surface of the shaft, respectively. Inertial $(X,Y,Z)$ and rotating $(x,y,z)$ coordinate systems are shown in Figure 2. The shaft rotates along its longitudinal $x$-axis at a constant rate $\Omega$. $\theta$ represents the ply angle measured by the s-axis positive coordinate.

[figure omitted; refer to PDF] [figure omitted; refer to PDF]

2.1. Stress-Strain Formulae

The stress-strain relations of a lamina shown in Figure 3 [13] in the principle material directions are as follows [25]:$$\begin{array}{}\text{(1)}& \left\{\begin{array}{c}{\sigma }_{11}\\ {\sigma }_{22}\\ {\sigma }_{33}\\ {\tau }_{12}\\ {\tau }_{13}\\ {\tau }_{23}\end{array}\right\}=\left[\begin{array}{cccccc}{Q}_{11}& Q_{12}& Q_{13}& 0& 0& 0\\ Q_{12}& Q_{22}& Q_{23}& 0& 0& 0\\ Q_{13}& Q_{23}& Q_{33}& 0& 0& 0\\ 0& 0& 0& Q_{66}& 0& 0\\ 0& 0& 0& 0& Q_{55}& 0\\ 0& 0& 0& 0& 0& Q_{44}\end{array}\right]\left\{\begin{array}{c}{\epsilon }_{11}\\ {\epsilon }_{22}\\ {\epsilon }_{33}\\ {\gamma }_{12}\\ {\gamma }_{13}\\ {\gamma }_{23}\end{array}\right\}\end{array}$$

[figure omitted; refer to PDF]

The above formula can be abbreviated as$$\begin{array}{}\text{(2)}& \left\{\sigma \right\}=\left[Q\right]\left\{\epsilon \right\}\end{array}$$where [Q] is the stiffness matrix and the material constants $Q_{ij}$ can be determined according to the material properties of orthogonal anisotropic laminates [25].$$\begin{array}{c}\text{(3)}& Q_{\mathrm{11}}=E_{\mathrm{11}}\frac{\mathrm{1}-v_{\mathrm{23}}{v}_{\mathrm{32}}}{\Delta },Q_{\mathrm{22}}=E_{\mathrm{22}}\frac{\mathrm{1}-v_{\mathrm{31}}{v}_{\mathrm{13}}}{\Delta }\\ Q_{\mathrm{33}}=E_{\mathrm{33}}\frac{\mathrm{1}-v_{\mathrm{12}}{v}_{\mathrm{21}}}{\Delta },\\ Q_{\mathrm{44}}={\kappa G}_{\mathrm{23}},\\ Q_{\mathrm{55}}={\kappa G}_{\mathrm{13}},\\ Q_{\mathrm{66}}={\kappa G}_{\mathrm{12}}\\ Q_{\mathrm{12}}=E_{\mathrm{11}}\frac{v_{\mathrm{21}}+v_{\mathrm{31}}{v}_{\mathrm{23}}}{\Delta }=E_{\mathrm{22}}\frac{v_{\mathrm{12}}+v_{\mathrm{32}}{v}_{\mathrm{13}}}{\Delta }\\ Q_{\mathrm{13}}=E_{\mathrm{11}}\frac{v_{\mathrm{31}}+v_{\mathrm{21}}{v}_{\mathrm{32}}}{\Delta }=E_{\mathrm{22}}\frac{v_{\mathrm{13}}+v_{\mathrm{12}}{v}_{\mathrm{23}}}{\Delta }\\ Q_{\mathrm{23}}=E_{\mathrm{22}}\frac{v_{\mathrm{32}}+v_{\mathrm{12}}{v}_{\mathrm{31}}}{\Delta }=E_{\mathrm{33}}\frac{v_{\mathrm{23}}+v_{\mathrm{21}}{v}_{\mathrm{13}}}{\Delta }\\ \Delta =\mathrm{1}-v_{\mathrm{12}}{v}_{\mathrm{21}}-v_{\mathrm{23}}{v}_{\mathrm{32}}-v_{\mathrm{31}}{v}_{\mathrm{13}}-\mathrm{2}{v}_{\mathrm{21}}{v}_{\mathrm{32}}{v}_{\mathrm{13}}\end{array}$$where $E_{ij}$ are the elastic moduli, $G_{ij}$ are the shear moduli, $v_{ij}$ are Poisson’s ratios, and $\kappa$ is the Timoshenko shear coefficient [26].

Consider an arbitrary layer of the laminate whose fiber direction makes an angle $\eta$ with respect to the x-axis of the cylindrical coordinate system $(x,s,\xi )$ as shown in Figure 4 [13]. The stress-strain relationship can be written as [25]$$\begin{array}{}\text{(4)}& \left\{\begin{array}{c}\begin{array}{c}{\sigma }_{xx}\\ {\sigma }_{ss}\\ {\sigma }_{\xi \xi }\\ {\tau }_{\xi s}\\ {\tau }_{x\xi }\\ {\tau }_{xs}\end{array}\end{array}\right\}=\left[\begin{array}{cccccc}{\overline{Q}}_{\mathrm{11}}& {\overline{Q}}_{\mathrm{12}}& {\overline{Q}}_{\mathrm{13}}& \mathrm{0}& \mathrm{0}& {\overline{Q}}_{\mathrm{16}}\\ {\overline{Q}}_{\mathrm{12}}& {\overline{Q}}_{\mathrm{22}}& {\overline{Q}}_{\mathrm{23}}& \mathrm{0}& \mathrm{0}& {\overline{Q}}_{\mathrm{26}}\\ {\overline{Q}}_{\mathrm{13}}& {\overline{Q}}_{\mathrm{23}}& {\overline{Q}}_{\mathrm{33}}& \mathrm{0}& \mathrm{0}& {\overline{Q}}_{\mathrm{36}}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& {\overline{Q}}_{\mathrm{44}}& {\overline{Q}}_{\mathrm{45}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& {\overline{Q}}_{\mathrm{45}}& {\overline{Q}}_{\mathrm{55}}& \mathrm{0}\\ {\overline{Q}}_{\mathrm{16}}& {\overline{Q}}_{\mathrm{26}}& {\overline{Q}}_{\mathrm{36}}& \mathrm{0}& \mathrm{0}& {\overline{Q}}_{\mathrm{66}}\end{array}\right]\left\{\begin{array}{c}{\epsilon }_{xx}\\ {\epsilon }_{ss}\\ {\epsilon }_{\xi \xi }\\ {\gamma }_{\xi s}\\ {\gamma }_{x\xi }\\ {\gamma }_{xs}\end{array}\right\}\end{array}$$

[figure omitted; refer to PDF]

The above equation can be abbreviated as$$\begin{array}{}\text{(5)}& \left\{\sigma \right\}=\left[\overline{Q}\right]\left\{\epsilon \right\}\end{array}$$where $\left[\overline{Q}\right]$ is the transformed stiffness matrix of the layer, which can be obtained by transformation matrix [T]:$$\begin{array}{}\text{(6)}& \left[\overline{Q}\right]={\left[\mathrm{T}\right]}^{-\mathrm{1}}\left[Q\right]{\left[\mathrm{T}\right]}^{-\mathrm{T}}\end{array}$$in which$$\begin{array}{}\text{(7)}& \left[T\right]=\left[\begin{array}{cccccc}{c}^{\mathrm{2}}& s^{\mathrm{2}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{2}cs\\ s^{\mathrm{2}}& c^{\mathrm{2}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& -\mathrm{2}cs\\ \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& c& -s& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& s& c& \mathrm{0}\\ -cs& cs& \mathrm{0}& \mathrm{0}& \mathrm{0}& \left(c^{\mathrm{2}}-s^{\mathrm{2}}\right)\end{array}\right],\hspace{1em}c=\cos \eta ,s=\sin \eta \end{array}$$

The displacement field of the shaft based on the refined VAM thin-walled beam theory is assumed in the form [22]:$$\begin{array}{c}\text{(8)}& u\left(x,y,z,t\right)=U\left(x,t\right)-y\left(s\right){\psi }_y\left(x,t\right)-z\left(s\right){\psi }_z\left(x,t\right)+g\left(s,x,t\right)v\left(x,y,z,t\right)=V\left(x,t\right)-z\left(s\right)\varphi \left(x,t\right)\\ w\left(x,y,z,t\right)=W\left(x,t\right)+y\left(s\right)\varphi \left(x,t\right)\end{array}$$where $U\left(x,t\right)$, $V\left(x,t\right)$, and $W\left(x,t\right)$ denote the rigid body displacements in the $x,y,z$ directions, while $\varphi \left(x,t\right)$, ${\psi }_y\left(x,t\right)$, and ${\psi }_z\left(x,t\right)$ are the rotation angles of the cross-section about x-, y-, and z-axis, respectively. $g\left(s,x,t\right)$ is the warping function, which is modified as$$\begin{array}{}\text{(9)}& g\left(s,x,t\right)=G\left(s\right){\varphi }^{\prime }\left(x,t\right)+g_{\mathrm{1}}\left(s\right)U^{\prime }\left(x,t\right)+g_{\mathrm{2}}\left(s\right){\psi }_y^{\prime }\left(x,t\right)+g_{\mathrm{3}}\left(s\right){\psi }_z^{\prime }\left(x,t\right)\end{array}$$where $G\left(s\right),g_1\left(s\right),g_2\left(s\right)$, and $g_3(s)$ are related to physical behavior of the torsion twist rate, the axial strain, and the bending curvatures, respectively.

In the ($x,y,z$) frame,${\psi }_y\left(x,t\right)$ and ${\psi }_z\left(x,t\right)$can be expressed as$$\begin{array}{c}\text{(10)}& {\psi }_y\left(x,t\right)=V^{\prime }\left(x,t\right)-\mathrm{2}{\gamma }_{zx}{\psi }_z\left(x,t\right)=W^{\prime }\left(x,t\right)-\mathrm{2}{\gamma }_{yx}\end{array}$$where ${\gamma }_{zx}$ and ${\gamma }_{yx}$ are the transverse shear strains in the planes of $xz$ and $xy$, respectively.

According to the displacement representations (8), (9), and (10), the strains of the composite shaft are obtained by$$\begin{array}{c}\text{(11)}& {\epsilon }_{xx}=U^{\prime }\left(x,t\right)-y\left(s\right){\psi }_y^{\prime }\left(x,t\right)-z\left(s\right){\psi }_z^{\prime }\left(x,t\right)\mathrm{2}{\gamma }_{xs}=\frac{dg}{ds}+r_n{\varphi }^{\prime }\left(x,t\right)+\left(V^{\prime }\left(x,t\right)-{\psi }_y\left(x,t\right)\frac{dy}{ds}+\left(W^{\prime }\left(x,t\right)-{\psi }_z\left(x,t\right)\right)\frac{dz}{ds}\right)\\ \mathrm{2}{\gamma }_{x\xi }=\left(V^{\prime }\left(x,t\right)-{\psi }_y\left(x,t\right)\frac{dz}{ds}+\left(W^{\prime }\left(x,t\right)-{\psi }_z\left(x,t\right)\right)\frac{dy}{ds}\right)\end{array}$$where $r_n$ denotes the normal projection of $\mathit{r}$ of an arbitrary point on the cross section of the deformed shaft in the normal direction:$$\begin{array}{}\text{(12)}& r_n=y\left(s\right)\frac{dz\left(s\right)}{ds}-z\left(s\right)\frac{dy\left(s\right)}{ds}\end{array}$$

2.2. Strain Energy of Composite Shaft

The strain energy of the composite shaft is$$\begin{array}{}\text{(13)}& U=\frac{\mathrm{1}}{\mathrm{2}}{{\displaystyle \int }}_{\mathrm{0}}^L{{\displaystyle \iint }}_A^{\hspace{0.17em}}\left({\sigma }_{xx}{\epsilon }_{xx}+{\tau }_{xs}{\gamma }_{xs}+{\tau }_{x\xi }{\gamma }_{x\xi }\right)dAdx\end{array}$$where ${\sigma }_{xx}$, ${\tau }_{xs}$, and ${\tau }_{x\xi }$ are the engineering stresses associated with the engineering strains ${\epsilon }_{xx}$, ${\gamma }_{xs}$, and ${\gamma }_{x\xi }$.

Taking variation of the above strain energy expression, (14) is obtained:$$\begin{array}{}\text{(14)}& \delta U={{\displaystyle \int }}_{\mathrm{0}}^L{{\displaystyle \iint }}_A^{\hspace{0.17em}}\left({\sigma }_{xx}{\delta \epsilon }_{xx}+{\tau }_{xs}{\delta \gamma }_{xs}+{\tau }_{x\xi }\delta {\gamma }_{x\xi }\right)dAdx={{\displaystyle \int }}_{\mathrm{0}}^L\left[F_x^{\prime }\delta U+Q_y^{\prime }\delta V+Q_Z^{\prime }\delta W+M_x^{\prime }\delta \varphi +\left(M_z^{\prime }+Q_y\right)\delta {\psi }_y+\left(M_y^{\prime }+Q_z\right)\delta {\psi }_z\right]dx\end{array}$$

Define the stress resultants $F_x$,$Q_y,Q_z$ and stress couples $M_x,M_y,M_z$ as$$\begin{array}{c}\text{(15)}& F_x=\oint N_{xx}ds{M}_x=\oint N_{xs}{r}_{n}ds\\ M_y=-\oint N_{xx}z\left(s\right)ds\\ M_z=-\oint N_{xx}y\left(s\right)ds\\ Q_y=\oint \left[N_{xs}\frac{dy\left(s\right)}{ds}+N_{x\xi }\frac{dz\left(s\right)}{ds}\right]ds\\ Q_z=\oint \left[N_{xs}\frac{dz\left(s\right)}{ds}+N_{x\xi }\frac{dy\left(s\right)}{ds}\right]ds\end{array}$$where $N_{xx}$, $N_{xs}$, and $N_{x\xi }$ are shell stress resultants, which are defined according to the following expressions:$$\begin{array}{}\text{(16)}& \left[\begin{array}{c}{N}_{xx}\\ N_{xs}\\ N_{x\xi }\end{array}\right]=\left[\begin{array}{ccc}A\left(s\right)& \frac{B\left(s\right)}{\mathrm{2}}& \mathrm{0}\\ \frac{B\left(s\right)}{\mathrm{2}}& \frac{C\left(s\right)}{\mathrm{4}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \kappa D\left(s\right)\end{array}\right]\left[\begin{array}{c}{\epsilon }_{xx}\\ \mathrm{2}{\gamma }_{xs}\\ {\gamma }_{x\xi }\end{array}\right]\end{array}$$in which$$\begin{array}{c}\text{(17)}& A\left(s\right)=A_{\mathrm{11}}-\frac{A_{\mathrm{12}}^{\mathrm{2}}}{A_{\mathrm{22}}},B\left(s\right)=\mathrm{2}\left[A_{\mathrm{16}}-\frac{A_{\mathrm{12}}{A}_{\mathrm{26}}}{A_{\mathrm{22}}}\right],\\ C\left(s\right)=\mathrm{4}\left[A_{\mathrm{66}}-\frac{{\left(A_{26}\right)}^2}{A_{\mathrm{22}}}\right],\\ D\left(s\right)=\left[{\mathrm{A}}_{\mathrm{44}}-\frac{{\left({\mathrm{A}}_{\mathrm{45}}\right)}^2}{{\mathrm{A}}_{\mathrm{55}}}\right],\\ A_{ij}={{\displaystyle \int }}_{-h/\mathrm{2}}^{h/\mathrm{2}}{\overline{Q}}_{ij}d\xi ={\displaystyle \sum }_{k=\mathrm{1}}^N{\overline{Q}}_{ij}^k\left(z_k-z_{k-\mathrm{1}}\right)\hspace{1em}\left(i,j=\mathrm{1,2},\mathrm{6};i,j=\mathrm{4,5}\right)\end{array}$$

Combining (11) and (16), (15) can be expressed in the following form:$$\begin{array}{c}\text{(18)}& F_x=C_{\mathrm{11}}{U}^{\prime }+C_{\mathrm{12}}{\varphi }^{\prime }+C_{\mathrm{13}}{\psi }_z^{\prime }+C_{\mathrm{14}}{\psi }_y^{\prime }+C_{\mathrm{15}}\left(V^{\prime }-{\psi }_y\right)+C_{\mathrm{16}}\left(W^{\prime }-{\psi }_Z\right)M_x=C_{\mathrm{12}}{U}^{\prime }+C_{\mathrm{22}}{\varphi }^{\prime }+C_{\mathrm{23}}{\psi }_z^{\prime }+C_{\mathrm{24}}{\psi }_y^{\prime }+C_{\mathrm{25}}\left(V^{\prime }-{\psi }_y\right)+C_{\mathrm{26}}\left(W^{\prime }-{\psi }_Z\right)\\ M_y=C_{\mathrm{13}}{U}^{\prime }+C_{\mathrm{23}}{\varphi }^{\prime }+C_{\mathrm{33}}{\psi }_z^{\prime }+C_{\mathrm{34}}{\psi }_y^{\prime }+C_{\mathrm{35}}\left(V^{\prime }-{\psi }_y\right)+C_{\mathrm{36}}\left(W^{\prime }-{\psi }_Z\right)\\ M_z=C_{\mathrm{14}}{U}^{\prime }+C_{\mathrm{24}}{\varphi }^{\prime }+C_{\mathrm{34}}{\psi }_z^{\prime }+C_{\mathrm{44}}{\psi }_y^{\prime }+C_{\mathrm{45}}\left(V^{\prime }-{\psi }_y\right)+C_{\mathrm{46}}\left(W^{\prime }-{\psi }_Z\right)\\ Q_y=C_{\mathrm{15}}{U}^{\prime }+C_{\mathrm{25}}{\varphi }^{\prime }+C_{\mathrm{35}}{\psi }_z^{\prime }+C_{\mathrm{45}}{\psi }_y^{\prime }+C_{\mathrm{55}}\left(V^{\prime }-{\psi }_y\right)+C_{\mathrm{56}}\left(W^{\prime }-{\psi }_Z\right)\\ Q_z=C_{\mathrm{16}}{U}^{\prime }+C_{\mathrm{26}}{\varphi }^{\prime }+C_{\mathrm{36}}{\psi }_z^{\prime }+C_{\mathrm{46}}{\psi }_y^{\prime }+C_{\mathrm{56}}\left(V^{\prime }-{\psi }_y\right)+C_{\mathrm{66}}\left(W^{\prime }-{\psi }_Z\right)\end{array}$$where $C_{ij}=C_{ji}\hspace{0.17em}\hspace{0.17em}(i,j=1,\cdots ,6)$ are the stiffness coefficients of the composite shaft, which can be expressed according to the cross section geometry and material properties as follows:$$\begin{array}{c}\text{(19)}& C_{\mathrm{11}}={{\displaystyle \oint }}_{\Gamma }^{\hspace{0.17em}}\left(A-\frac{B^{\mathrm{2}}}{C}\right)ds+\left\{\frac{{\left[{\oint }_{\Gamma }^{\hspace{0.17em}}\left(B/C\right)ds\right]}^{\mathrm{2}}}{{\oint }_{\Gamma }^{\hspace{0.17em}}\left(\mathrm{1}/C\right)ds}\right\},C_{\mathrm{12}}=\left[\frac{{\oint }_{\Gamma }^{\hspace{0.17em}}\left(B/C\right)ds}{{\oint }_{\Gamma }^{\hspace{0.17em}}\left(\mathrm{1}/C\right)ds}\right]A_e,\\ C_{\mathrm{13}}=-{{\displaystyle \oint }}_{\Gamma }^{\hspace{0.17em}}\left(A-\frac{B^{\mathrm{2}}}{C}\right)zds-\left\{\frac{\left[{\oint }_{\Gamma }^{\hspace{0.17em}}\left(B/C\right)ds\left[{\oint }_{\Gamma }^{\hspace{0.17em}}\left(B/C\right)zds\right]\right]}{{\oint }_{\Gamma }^{\hspace{0.17em}}\left(\mathrm{1}/C\right)ds}\right\},\\ C_{\mathrm{14}}=-{{\displaystyle \oint }}_{\Gamma }^{\hspace{0.17em}}\left(A-\frac{B^{\mathrm{2}}}{C}\right)yds-\left\{\frac{\left[{\oint }_{\Gamma }^{\hspace{0.17em}}\left(B/C\right)ds\left[{\oint }_{\Gamma }^{\hspace{0.17em}}\left(B/C\right)yds\right]\right]}{{\oint }_{\Gamma }^{\hspace{0.17em}}\left(\mathrm{1}/C\right)ds}\right\},\\ C_{\mathrm{15}}=\frac{\mathrm{1}}{\mathrm{2}}{{\displaystyle \oint }}_{\Gamma }^{\hspace{0.17em}}B\frac{dy}{ds}ds,\\ C_{\mathrm{16}}=\frac{\mathrm{1}}{\mathrm{2}}{{\displaystyle \oint }}_{\Gamma }^{\hspace{0.17em}}B\frac{dz}{ds}ds,\\ C_{\mathrm{22}}=\left[\frac{\mathrm{1}}{{\oint }_{\Gamma }^{\hspace{0.17em}}\left(\mathrm{1}/C\right)ds}\right]A_e^{\mathrm{2}},\\ C_{\mathrm{23}}=-\left[\frac{{\oint }_{\Gamma }^{\hspace{0.17em}}\left(B/C\right)zds}{{\oint }_{\Gamma }^{\hspace{0.17em}}\left(\mathrm{1}/C\right)ds}\right]A_e,\\ C_{\mathrm{24}}=-\left[\frac{{\oint }_{\Gamma }^{\hspace{0.17em}}\left(B/C\right)yds}{{\oint }_{\Gamma }^{\hspace{0.17em}}\left(\mathrm{1}/C\right)ds}\right]A_e,\\ C_{\mathrm{25}}=\frac{\mathrm{1}}{\mathrm{4}}{{\displaystyle \oint }}_{\Gamma }^{\hspace{0.17em}}{r}_{n}C\frac{dy}{ds}ds,\\ C_{\mathrm{26}}=\frac{\mathrm{1}}{\mathrm{4}}{{\displaystyle \oint }}_{\Gamma }^{\hspace{0.17em}}{r}_{n}C\frac{dz}{ds}ds,\\ C_{\mathrm{33}}={{\displaystyle \oint }}_{\Gamma }^{\hspace{0.17em}}\left(A-\frac{B^{\mathrm{2}}}{C}\right)z^{\mathrm{2}}ds+\left\{\frac{{\left[{\oint }_{\Gamma }^{\hspace{0.17em}}\left(B/C\right)zds\right]}^{\mathrm{2}}}{{\oint }_{\Gamma }^{\hspace{0.17em}}\left(\mathrm{1}/C\right)ds}\right\},\\ C_{\mathrm{34}}={{\displaystyle \oint }}_{\Gamma }^{\hspace{0.17em}}\left(A-\frac{B^{\mathrm{2}}}{C}\right)yzds+\left\{\frac{\left[{\oint }_{\Gamma }^{\hspace{0.17em}}\left(B/C\right)yds{\oint }_{\Gamma }^{\hspace{0.17em}}\left(B/C\right)zds\right]}{{\oint }_{\Gamma }^{\hspace{0.17em}}\left(\mathrm{1}/C\right)ds}\right\},\\ C_{\mathrm{35}}=-\frac{\mathrm{1}}{\mathrm{2}}{{\displaystyle \oint }}_{\Gamma }^{\hspace{0.17em}}Bz\frac{dy}{ds}ds,\\ C_{\mathrm{36}}=-\frac{\mathrm{1}}{\mathrm{2}}{{\displaystyle \oint }}_{\Gamma }^{\hspace{0.17em}}Bz\frac{dz}{ds}ds,\\ C_{\mathrm{44}}=\kappa {{\displaystyle \oint }}_{\Gamma }^{\hspace{0.17em}}\left(A-\frac{B^{\mathrm{2}}}{C}\right)y^{\mathrm{2}}ds+\left\{\frac{{\left[{\oint }_{\Gamma }^{\hspace{0.17em}}\left(B/C\right)yds\right]}^{\mathrm{2}}}{{\oint }_{\Gamma }^{\hspace{0.17em}}\left(\mathrm{1}/C\right)ds}\right\},\\ C_{\mathrm{45}}=-\frac{\mathrm{1}}{\mathrm{2}}\kappa {{\displaystyle \oint }}_{\Gamma }^{\hspace{0.17em}}By\frac{dy}{ds}ds,\\ C_{\mathrm{46}}=-\frac{\mathrm{1}}{\mathrm{2}}\kappa {{\displaystyle \oint }}_{\Gamma }^{\hspace{0.17em}}By\frac{dz}{ds}ds,\\ C_{\mathrm{55}}=\kappa {{\displaystyle \oint }}_{\Gamma }^{\hspace{0.17em}}\frac{\mathrm{1}}{\mathrm{4}}C{\left(\frac{dy}{ds}\right)}^{\mathrm{2}}+D{\left(\frac{dz}{ds}\right)}^{\mathrm{2}}ds,\\ C_{\mathrm{56}}=\kappa {{\displaystyle \oint }}_{\Gamma }^{\hspace{0.17em}}\left(\frac{\mathrm{1}}{\mathrm{4}}C-D\right)\frac{dy}{ds}\frac{dz}{ds}ds,\\ C_{\mathrm{66}}=\kappa {{\displaystyle \oint }}_{\Gamma }^{\hspace{0.17em}}\frac{\mathrm{1}}{\mathrm{4}}C{\left(\frac{dz}{ds}\right)}^{\mathrm{2}}+D{\left(\frac{dy}{ds}\right)}^{\mathrm{2}}ds,\\ A_e=\frac{\mathrm{1}}{\mathrm{2}}\kappa {{\displaystyle \oint }}_{\Gamma }^{\hspace{0.17em}}\left(y\frac{dz}{ds}-z\frac{dy}{ds}\right)ds\end{array}$$where ${\oint }_{\Gamma }^{\hspace{0.17em}}\left(\bullet \right)ds$ denotes the integral around the loop of the midline cross section; $A_e$ is the surrounded area of the midline of the section.

2.3. Kinetic Energy of Composite Shaft

The position vector of an arbitrary point on the composite shaft can be written as$$\begin{array}{}\text{(20)}& \mathit{r}=\left(y+v\right)\mathit{i}\mathbf{+}\left(z+w\right)\mathit{j}\mathbf{+}\left(x+u\right)\mathit{k}\end{array}$$

From the above equation, the velocity of an arbitrary point can be given by$$\begin{array}{}\text{(21)}& \mathit{V}=\dot{\mathit{r}}=\left(\dot{v}-\Omega \left(z+w\right)\right)\mathit{i}\mathbf{+}\left(\dot{w}+\Omega \left(y+v\right)\right)\mathit{j}\mathbf{+}\dot{u}\mathit{k}\end{array}$$

The kinetic energy of the composite spinning shaft can be written as$$\begin{array}{}\text{(22)}& T=\frac{\mathrm{1}}{\mathrm{2}}{{\displaystyle \int }}_{\mathrm{0}}^L{{\displaystyle \iint }}_A^{\hspace{0.17em}}\rho \left(\mathit{V}\bullet \mathit{V}\right)dAdx\end{array}$$

The variation of kinetic energy can be written as$$\begin{array}{}\text{(23)}& \delta T=-{{\displaystyle \int }}_{\mathrm{0}}^L\left(I_{\mathrm{1}}\delta U+I_{\mathrm{2}}\delta V+I_{\mathrm{3}}\delta W+I_{\mathrm{4}}{\psi }_y+I_{\mathrm{5}}{\psi }_z+I_{\mathrm{6}}\varphi \right)dx\end{array}$$in which$$\begin{array}{c}\text{(24)}& I_1={\mathrm{m}}_{\mathrm{c}}\ddot{\mathrm{U}}-{\mathrm{S}}_{\mathrm{z}}{\ddot{\psi }}_y-{\mathrm{S}}_{\mathrm{y}}{\ddot{\psi }}_z{I}_{\mathrm{2}}=m_c\left(\ddot{V}-\mathrm{2}\Omega \dot{W}-{\Omega }^{\mathrm{2}}V\right)-S_z\left(\mathrm{2}\Omega \dot{\varphi }+{\Omega }^{\mathrm{2}}\right)-S_y\left(\ddot{\varphi }-{\Omega }^{\mathrm{2}}\varphi \right)\\ I_{\mathrm{3}}=m_c\left(\ddot{W}-\mathrm{2}\Omega \dot{V}-{\Omega }^{\mathrm{2}}W\right)+S_z\left(\ddot{\varphi }+{\Omega }^{\mathrm{2}}\varphi \right)-S_y\left(\mathrm{2}\Omega \dot{\varphi }+{\Omega }^{\mathrm{2}}\right)\\ I_{\mathrm{4}}=S_z\ddot{U}-I_z{\ddot{\psi }}_y-I_{yz}{\ddot{\psi }}_z\\ I_{\mathrm{5}}=S_y\ddot{W}-I_{yz}{\ddot{\psi }}_y-I_y{\ddot{\psi }}_z\\ I_{\mathrm{6}}=S_z\left(\ddot{W}+\mathrm{2}\Omega \dot{U}-{\Omega }^{\mathrm{2}}W\right)+\left(I_y+I_z\right)\left(\ddot{\varphi }-{\Omega }^{\mathrm{2}}\varphi \right)-S_y\left(\ddot{V}-\mathrm{2}\Omega \dot{W}-{\Omega }^{\mathrm{2}}V\right)\\ m_c={{\displaystyle \iint }}_A^{\hspace{0.17em}}\rho dA,\\ S_z={{\displaystyle \iint }}_A^{\hspace{0.17em}}\rho ydA,\\ S_y={{\displaystyle \iint }}_A^{\hspace{0.17em}}\rho zdA,\\ I_z={{\displaystyle \iint }}_A^{\hspace{0.17em}}\rho y^{\mathrm{2}}dA,\\ I_y={{\displaystyle \iint }}_A^{\hspace{0.17em}}\rho z^{\mathrm{2}}dA,\\ I_{yz}={{\displaystyle \iint }}_A^{\hspace{0.17em}}\rho yzdA\end{array}$$where $I_i\left(i=1,\cdots ,6\right)$ are equivalent cross section mass coefficients.

2.4. Governing Equations

To derive the vibration equation of the composite thin-walled shaft, Hamilton's principle was applied:$$\begin{array}{}\text{(25)}& {{\displaystyle \int }}_{t_{\mathrm{0}}}^{t_{\mathrm{1}}}\left(\delta U-\delta T\right)dt=\mathrm{0}\end{array}$$

Combining (14), (23), and (25), the governing equations of the composite shaft are obtained:$$\begin{array}{c}\text{(26)}& -F_x^{\prime }+I_{\mathrm{1}}=\mathrm{0}-Q_y^{\prime }+I_{\mathrm{2}}=\mathrm{0}\\ -Q_z^{\prime }+I_{\mathrm{3}}=\mathrm{0}\\ -M_z^{\prime }-Q_y+I_{\mathrm{4}}=\mathrm{0}\\ -M_y^{\prime }-Q_z+I_{\mathrm{5}}=\mathrm{0}\\ -M_x^{\prime }+I_{\mathrm{6}}=\mathrm{0}\end{array}$$

Substituting internal forces into (26), the free vibration equations can be derived. The special properties of the circular cross section shaft with CUS configuration cause some of stiffness coefficients to be zero. The equations of motion involving variables in terms of displacements can be reduced as$$\begin{array}{}\text{(27a)}& f_{\mathrm{1}}\left(U,\varphi \right)=-C_{\mathrm{11}}{U}^{\prime \prime }-C_{\mathrm{12}}{\varphi }^{\prime \prime }+m_c\ddot{U}=\mathrm{0}\text{(27b)}& f_{\mathrm{2}}\left(U,\varphi \right)=-C_{\mathrm{12}}{U}^{\prime \prime }-C_{\mathrm{22}}{\varphi }^{\prime \prime }+\left(I_y+I_z\right)\left(\ddot{\varphi }-{\Omega }^{\mathrm{2}}\varphi \right)=\mathrm{0}\text{(27c)}& f_{\mathrm{3}}\left(V,W,{\psi }_y,{\psi }_z\right)=-C_{\mathrm{35}}{\psi }_z^{\prime \prime }-C_{\mathrm{55}}\left(V^{\prime \prime }-{\psi }_y^{\prime }\right)+m_c\left(\ddot{V}-\mathrm{2}\Omega \dot{W}-{\Omega }^{\mathrm{2}}V\right)=\mathrm{0}\text{(27d)}& f_{\mathrm{4}}\left(V,W,{\psi }_y,{\psi }_z\right)=-C_{\mathrm{46}}{\psi }_y^{\prime \prime }-C_{\mathrm{66}}\left(W^{\prime \prime }-{\psi }_z^{\prime }\right)+m_c\left(\ddot{W}+\mathrm{2}\Omega \dot{V}-{\Omega }^{\mathrm{2}}W\right)=\mathrm{0}\text{(27e)}& f_{\mathrm{5}}\left(V,W,{\psi }_y,{\psi }_z\right)=-C_{\mathrm{44}}{\psi }_y^{\prime \prime }-C_{\mathrm{46}}\left(W^{\prime \prime }-{\psi }_z^{\prime }\right)-C_{\mathrm{35}}{\psi }_z^{\prime }-C_{\mathrm{55}}\left(V^{\prime }-{\psi }_y^{\hspace{0.17em}}\right)-I_z{\ddot{\psi }}_y=\mathrm{0}\text{(27f)}& f_{\mathrm{6}}\left(V,W,{\psi }_y,{\psi }_z\right)=-C_{\mathrm{33}}{\psi }_z^{\prime \prime }-C_{\mathrm{35}}\left(V^{\prime \prime }-{\psi }_y^{\prime }\right)-C_{\mathrm{46}}{\psi }_y^{\prime }-C_{\mathrm{66}}\left(W^{\prime }-{\psi }_z^{\hspace{0.17em}}\right)-I_y{\ddot{\psi }}_z=\mathrm{0}\end{array}$$

The vibration equations are decoupled as the tensile-torsion coupling system ((27a) and (27b)) and bending-transverse shear coupling system ( (27c)~(27f)). This paper only studies the latter.

3. GDQM Solution to Governing Equations

The GDQM was used to solve the governing equations. In the GDQM, derivative of any order of a function is approximated by a weighted linear sum of the function values at all the discrete points. Taking a function $f(x,t)$ as an example, the mathematic description of the GDQM is given by [27]$$\begin{array}{}\text{(28)}& {\left.\frac{{\partial }^{p}f\left(x,t\right)}{\partial x^p}\right|}_{x=x_j}={\displaystyle \sum }_{k=\mathrm{1}}^N{C}_{jk}^{p}f\left(x_k,t\right),\hspace{1em}j=\mathrm{1,2},\cdots ,N\end{array}$$where N is the number of total discrete grid points in the $x$ direction. $C_{jk}^p$ is the weighting coefficient related to the $p$th-order derivative, and the weighting coefficenct is obtained by the following.

If p= 1, then$$\begin{array}{c}\text{(29)}& A_{jk}^{\left(\mathrm{1}\right)}=\frac{M^{\left(\mathrm{1}\right)}\left(x_j\right)}{\left(x_j-x_k\right)M^{\left(\mathrm{1}\right)}\left(x_k\right)},\hspace{1em}j\ne k,j,k=\mathrm{1,2},\dots ,N{A}_{jj}^{\left(\mathrm{1}\right)}=-{\displaystyle \sum }_{k=\mathrm{1}\left(k\ne j\right)}^N{A}_{jk}^{\left(\mathrm{1}\right)},\hspace{1em}j=\mathrm{1,2},\cdots ,N\end{array}$$where $M^{\left(1\right)}\left(x\right)$ is the first derivative of $M(x)$ and they can be defined as$$\begin{array}{c}\text{(30)}& M\left(x\right)={\displaystyle \prod }_{k=\mathrm{1}}^N\left({x-x}_k\right),M^{\left(\mathrm{1}\right)}\left(x_k\right)={\displaystyle \prod }_{j=\mathrm{1}\left(j\ne k\right)}^N\left(x_k-x_j\right)\end{array}$$

If p>1, namely, for the second and higher order derivatives, the weighting coefficients are obtained by using the following simple recurrence relationship:$$\begin{array}{c}\text{(31)}& A_{jk}^{\left(n\right)}=p\left(A_{jk}^{\left(\mathrm{1}\right)}\bullet A_{jj}^{\left(n-\mathrm{1}\right)}-\frac{A_{jk}^{\left(n-\mathrm{1}\right)}}{x_j-x_k}\right),\hspace{1em}j\ne k,j,k=\mathrm{1,2},\dots ,N,n\ge \mathrm{2}{A}_{jj}^{\left(n\right)}=-{\displaystyle \sum }_{k=\mathrm{1}\left(k\ne j\right)}^N{A}_{jk}^{\left(n\right)},\hspace{1em}j=\mathrm{1,2},\cdots ,N\end{array}$$

To get denser population near boundaries, the sampling points are selected based on the Chebyshev–Gauss–Lobatto grid distribution.$$\begin{array}{}\text{(32)}& x_j=\frac{L}{\mathrm{2}}\left[\mathrm{1}-\cos \left(\frac{j-\mathrm{1}}{N-\mathrm{1}}\pi \right)\right],\hspace{1em}j=\mathrm{1,2},\cdots ,N\end{array}$$

To find the approximate solution of the composite shaft, assume that bending deformation and bending angle are of the following form:$$\begin{array}{c}\text{(33)}& v\left(x,t\right)=V\left(x\right)e^{i\omega t},w\left(x,t\right)=W\left(x\right)e^{i\omega t},\\ {\psi }_y\left(x,t\right)={\Psi }_y\left(x\right)e^{i\omega t},\\ {\psi }_z\left(x,t\right)={\Psi }_z\left(x\right)e^{i\omega t}\end{array}$$where $V\left(x\right),W\left(x\right),{\Psi }_y\left(x\right)$, and ${\Psi }_z\left(x\right)$ denote mode functions of the composite shaft; $\omega$ is the natural frequency; $i=\sqrt{-1}$.

Substituting (33) into the governing Equations (27c), (27d), (27e), and (27f), a group of ordinary differential equations with the coefficients varying toward the x-direction is produced:$$\begin{array}{}\text{(34)}& L^{\ast }{U}^{\ast }=\mathrm{0}\end{array}$$where ${U^{\ast }}^{\mathrm{T}}=\left\{V\left(x\right),W\left(x\right),{\Psi }_y\left(x\right),{\Psi }_z\left(x\right)\right\}$ and $L^{\ast }=\left[L_{ij}^{\ast }\right]\left(i,j=1,\cdots ,4\right)$ is a 4×4 matrix which includes the differential operative of $U^{\ast }$ defined as$$\begin{array}{}\text{(35)}& L^{\ast }=-{\omega }^{\mathrm{2}}\left[M\right]+i\omega \left[G\right]+\left[K\right]\end{array}$$where$$\begin{array}{c}\text{(36)}& \left[M\right]=\left[\begin{array}{cccc}{m}_c& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& m_c& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& -I_z& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& {-I}_y\end{array}\right]\left[G\right]=\left[\begin{array}{cccc}\mathrm{0}& -\mathrm{2}\Omega m_c& \mathrm{0}& \mathrm{0}\\ \mathrm{2}\Omega m_c& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\right]\\ \left[K\right]=\left[\begin{array}{cccc}-C_{\mathrm{55}}\frac{{\partial }^{\mathrm{2}}}{\partial x^{\mathrm{2}}}-m_c{\Omega }^{\mathrm{2}}& \mathrm{0}& C_{\mathrm{55}}\frac{\partial }{\partial x}& -C_{\mathrm{35}}\frac{{\partial }^{\mathrm{2}}}{\partial x^{\mathrm{2}}}\\ \mathrm{0}& -C_{\mathrm{66}}\frac{{\partial }^{\mathrm{2}}}{\partial x^{\mathrm{2}}}-m_c{\Omega }^{\mathrm{2}}& -C_{\mathrm{46}}\frac{{\partial }^{\mathrm{2}}}{\partial x^{\mathrm{2}}}& C_{\mathrm{66}}\frac{\partial }{\partial x}\\ {-C}_{\mathrm{55}}\frac{\partial }{\partial x}& -C_{\mathrm{46}}\frac{{\partial }^{\mathrm{2}}}{\partial x^{\mathrm{2}}}& -C_{\mathrm{44}}\frac{{\partial }^{\mathrm{2}}}{\partial x^{\mathrm{2}}}+C_{\mathrm{55}}& \left(C_{\mathrm{46}}-C_{\mathrm{35}}\right)\frac{\partial }{\partial x}\\ -C_{\mathrm{35}}\frac{{\partial }^{\mathrm{2}}}{\partial x^{\mathrm{2}}}& -C_{\mathrm{66}}\frac{\partial }{\partial x}& \left(C_{\mathrm{35}}-C_{\mathrm{46}}\right)\frac{\partial }{\partial x}& -C_{\mathrm{33}}\frac{{\partial }^{\mathrm{2}}}{\partial x^{\mathrm{2}}}+C_{\mathrm{66}}\end{array}\right]\end{array}$$

Considering (28), the differential quadrature discrete expression of the governing equations of the composite shaft can be obtained:$$\begin{array}{}\text{(37a)}& -C_{\mathrm{55}}{\displaystyle \sum }_{j=\mathrm{1}}^N{P}_{ij}^{\left(\mathrm{2}\right)}V\left(x_j\right)-m_c{\Omega }^{\mathrm{2}}V\left(x_i\right)+C_{\mathrm{55}}{\displaystyle \sum }_{j=\mathrm{1}}^N{P}_{ij}^{\left(\mathrm{1}\right)}{\Psi }_y\left(x_j\right)-C_{\mathrm{35}}{\displaystyle \sum }_{j=\mathrm{1}}^N{P}_{ij}^{\left(\mathrm{2}\right)}{\Psi }_z\left(x_j\right)-i\omega \mathrm{2}\Omega m_{c}W\left(x_i\right)-{\omega }^{\mathrm{2}}{m}_{c}V\left(x_i\right)=\mathrm{0}\text{(37b)}& -C_{\mathrm{66}}{\displaystyle \sum }_{j=\mathrm{1}}^N{P}_{ij}^{\left(\mathrm{2}\right)}W\left(x_j\right)-m_c{\Omega }^{\mathrm{2}}W\left(x_i\right)-C_{\mathrm{46}}{\displaystyle \sum }_{j=\mathrm{1}}^N{P}_{ij}^{\left(\mathrm{2}\right)}{\Psi }_y\left(x_j\right)+C_{\mathrm{66}}{\displaystyle \sum }_{j=\mathrm{1}}^N{P}_{ij}^{\left(\mathrm{1}\right)}{\Psi }_z\left(x_j\right)+i\omega \mathrm{2}\Omega m_{c}V\left(x_i\right)-{\omega }^{\mathrm{2}}{m}_{c}V\left(x_i\right)=\mathrm{0}\text{(37c)}& {-C}_{\mathrm{55}}{\displaystyle \sum }_{j=\mathrm{1}}^N{P}_{ij}^{\left(\mathrm{1}\right)}V\left(x_j\right)-C_{\mathrm{46}}{\displaystyle \sum }_{j=\mathrm{1}}^N{P}_{ij}^{\left(\mathrm{2}\right)}W\left(x_j\right)-C_{\mathrm{44}}{\displaystyle \sum }_{j=\mathrm{1}}^N{P}_{ij}^{\left(\mathrm{2}\right)}{\Psi }_y\left(x_j\right)+C_{\mathrm{55}}{\Psi }_y\left(x_i\right)+\left(C_{\mathrm{46}}-C_{\mathrm{35}}\right){\displaystyle \sum }_{j=\mathrm{1}}^N{P}_{ij}^{\left(\mathrm{1}\right)}{\Psi }_z\left(x_j\right)+{\omega }^{\mathrm{2}}{I_z\Psi }_y\left(x_i\right)=\mathrm{0}\text{(37d)}& -C_{\mathrm{35}}{\displaystyle \sum }_{j=\mathrm{1}}^N{P}_{ij}^{\left(\mathrm{2}\right)}V\left(x_j\right)-C_{\mathrm{66}}{\displaystyle \sum }_{j=\mathrm{1}}^N{P}_{ij}^{\left(\mathrm{1}\right)}W\left(x_j\right)+\left(C_{\mathrm{35}}-C_{\mathrm{46}}\right){\displaystyle \sum }_{j=\mathrm{1}}^N{P}_{ij}^{\left(\mathrm{1}\right)}{\Psi }_y\left(x_j\right)-C_{\mathrm{33}}{\displaystyle \sum }_{j=\mathrm{1}}^N{P}_{ij}^{\left(\mathrm{2}\right)}{\Psi }_z\left(x_j\right)+C_{\mathrm{66}}{\Psi }_z\left(x_i\right)+{\omega }^{\mathrm{2}}{I_y\Psi }_z\left(x_i\right)=\mathrm{0}\end{array}$$where $i=\mathrm{2},\cdots ,N-\mathrm{1}$.

The differential quadrature discrete expression can be simplified as$$\begin{array}{}\text{(38)}& \left(\left[M\right]{\omega }^{\mathrm{2}}\right)d+\left(\left[G\right]i\omega \right)d+\left[K_{dd}\right]d+\left[K_{db}\right]b=\mathrm{0}\end{array}$$

In the above equations, the dimensions of $\left[M\right]$, $\left[G\right]$, and $\left[K_{dd}\right]$ are $4\left(N-2\right)\times 4\left(N-2\right)$; the dimension of $\left[K_{db}\right]$ is $4\left(N-2\right)\times 8$; the dimensions of vectors $d$ and $b$ are $4\left(N-2\right)$ and 8, respectively, which can be expressed as$$\begin{array}{c}\text{(39)}& d^{\mathrm{T}}=\left\{V\left(x_{\mathrm{2}}\right),\cdots ,V\left(x_{N-\mathrm{1}}\right),W\left(x_{\mathrm{2}}\right),\cdots ,W\left(x_{N-\mathrm{1}}\right),{{\Psi }_y\left(x_{\mathrm{2}}\right),\cdots ,\Psi }_y\left(x_{N-\mathrm{1}}\right),{\Psi }_z\left(x_{\mathrm{2}}\right),\cdots ,{\Psi }_{\mathrm{z}}\left(x_{N-\mathrm{1}}\right)\right\}{b}^{\mathrm{T}}=\left\{V\left(x_{\mathrm{1}}\right),\cdots ,V\left(x_N\right),W\left(x_{\mathrm{1}}\right),\cdots ,W\left(x_N\right),{\Psi }_y\left(x_{\mathrm{1}}\right),\left.\cdots ,{\Psi }_y\left(x_N\right),{\Psi }_z\left(x_{\mathrm{1}}\right),\cdots ,{\Psi }_{\mathrm{z}}\left(x_N\right)\right\}\right\}\end{array}$$

Similarly, the discretized form of the boundary conditions becomes$$\begin{array}{}\text{(40)}& \left[K_{\mathit{b}\mathit{d}}\right]d+\left[K_{bb}\right]b=\mathrm{0}\end{array}$$

In the current application of the GDQM, three boundary conditions considered for shafts are

(a)

clamped at both edges (clamped-clamped, C-C)$$\begin{array}{c}\text{(41)}& V=\mathrm{0},W=\mathrm{0},\\ {\Psi }_y=\mathrm{0},\\ {\Psi }_z=\mathrm{0}\\ \left(x=\mathrm{0},L\right)\end{array}$$

(b)

simply supported at both edges (simply supported-simply supported, S-S)$$\begin{array}{c}\text{(42)}& V=\mathrm{0},W=\mathrm{0},\\ M_{\mathrm{y}}=\mathrm{0},\\ M_{\mathrm{z}}=\mathrm{0}\\ \left(x=\mathrm{0},L\right)\end{array}$$

(c)

cantilever shaft (clamped-free, C-F)$$\begin{array}{c}\text{(43)}& V=\mathrm{0},W=\mathrm{0},\\ {\Psi }_y=\mathrm{0},\\ {\Psi }_z=\mathrm{0}\\ \left(x=\mathrm{0}\right)\\ Q_{\mathrm{y}}=\mathrm{0},\\ Q_{\mathrm{z}}=\mathrm{0},\\ M_{\mathrm{y}}=\mathrm{0},\\ M_{\mathrm{z}}=\mathrm{0}\\ \left(x=L\right)\end{array}$$

Combine (40) and (38) to remove vector b, and then (44) is obtained:$$\begin{array}{}\text{(44)}& \left(\left[M\right]{\omega }^{\mathrm{2}}\right)d+\left(\left[G\right]i\omega \right)d+\left(\left[K_{dd}\right]\mathbf{-}\left[K_{db}\right]{\left[K_{bb}\right]}^{-\mathrm{1}}\left[K_{bd}\right]\right)d=\mathrm{0}\end{array}$$

Equation (44) can be abbreviated as$$\begin{array}{}\text{(45)}& \left(\left[M\right]{\omega }^{\mathrm{2}}\right)d+\left(\left[G\right]i\omega \right)d+\left[K\right]d=\mathrm{0}\end{array}$$

Equation (45) is a nonstandard eigenvalue equation that can be equivalently converted into a standard formula for a certain frequency [28]:$$\begin{array}{}\text{(46)}& \left(\left[\begin{array}{cc}\mathrm{0}& I\\ -K& -G\end{array}\right]-\omega \left[\begin{array}{cc}I& \mathrm{0}\\ \mathrm{0}& M\end{array}\right]\right)\left\{\begin{array}{c}d\\ \omega d\end{array}\right\}=\mathrm{0}\end{array}$$where I is an identity matrix and its dimension is 4(N − 2) × 4(N − 2).

By a mainstream method, the eigenvalue equation (46) can be solved, and the natural frequencies are achieved.

4. Numerical Results Analysis and Discussion

4.1. Convergence Analysis and Accuracy Verification

To validate the approximate calculation method in this paper, the composite hollow shafts made of the boron/epoxy composite material, which are considered by Bert and Kim [30], were investigated.

The properties of the material are $E_{11}=211$ GPa, $E_{22}=24.1$ GPa, $G_{12}=G_{13}=G_{23}=6.9$ GPa, $v_{12}=0.36$, and $\rho$ = 1967 kg/m³. The dimensions of the shaft are $L=2.47$ m, $r=0.06345$ m, and $h=1.321$ mm. The stacking sequence of the composite shaft is [90°/45°/45°/${0^{°}}_{\mathrm{6}}$/90°].

The effects of the number of grid points on the first five natural frequencies of the rotating composite shafts with various boundary conditions are shown in Figures 5–7, respectively. It can be seen from the figures that the first five natural frequencies converge quickly. When the number of grid points is larger than 12, the first five natural frequencies almost do not vary. This indicates that the convergence of the GDQM is good. Therefore, in later calculations, 12 uniform grid points are used.

[figure omitted; refer to PDF] [figure omitted; refer to PDF] [figure omitted; refer to PDF]

Table 1 presents the natural frequencies of cantilever composite shafts without shear deformation obtained from both present model and [29] at different rotating speeds. ${\Omega }^{\ast }$ and ${\omega }_i^{\ast }$ denote the nondimensional rotating speed and natural frequency, respectively. As shown in the Table, the numerical results agree well with those in [29].

Table 1

Natural frequencies of cantilever composite shafts.

$${\Omega }^{\ast }$$	$${\omega }_{\mathrm{1}}^{\ast }$$	$${\omega }_{\mathrm{2}}^{\ast }$$	$${\omega }_{\mathrm{3}}^{\ast }$$	$${\omega }_{\mathrm{4}}^{\ast }$$	$${\omega }_{\mathrm{5}}^{\ast }$$	$${\omega }_{\mathrm{6}}^{\ast }$$
0	3.5160	3.5160	22.0345	22.0345	61.6975	61.6975
Ref [29]	3.5160	3.5160	22.0340	22.0340	61.6970	61.6970
2	1.5160	5.5160	20.0345	24.0345	59.6975	63.6975
Ref [29]	1.5160	5.5160	20.0340	24.0340	59.6970	63.6970
3.5	0.0160	7.0160	18.5345	25.5345	58.1975	65.1975
Ref [29]	-	7.0160	18.5340	25.5340	58.1970	65.1970
4	0.4840	7.5160	18.0345	26.0345	57.6975	65.6975
Ref [29]	-	7.5160	18.0340	26.0340	57.6970	65.6970
8	4.4839	11.5160	14.0345	30.0345	53.6975	69.6975
Ref [29]	-	11.5160	14.0340	30.0340	53.6970	69.6970

4.2. Stability and Free Vibration Analysis

In the following example, the shaft is made of the graphite-epoxy composite material with parameters as follows:

$E_{11}$ = 206.8 GPa, $E_{22}$ = $E_{33}$ = 5.17 GPa, $G_{12}=3.1$ GPa, $G_{23}=G_{13}$ = 2.55 GPa, $v_{21}=v_{31}=0.00625$, $v_{32}=0.25$, and $\rho =1528.15$ kg/m³.

The dimensions of the shaft are $L=2.023$ m, $r=0.127$ m, and $h=0.381$ mm. The stacking sequence of the shaft is ${\left[\pm \theta \right]}_3$.

Figure 8 shows the first natural frequencies of the rotating composite shaft, which considers the shear deformation versus rotating speed for selected ply angles. When the rotating speed is zero, each ply angle corresponds to a single fundamental frequency because the bending mode frequencies in horizontal and vertical directions are the same for a circular cross section shaft. As soon as the rotation starts, a bifurcation of natural frequencies occurs because of the gyroscopic effect. Hence, a natural frequency curve splits into an upper one and a lower one. With the increase of rotational speed, the upper curve goes up corresponding to the forward feed (F), while the lower curve goes down corresponding to the backward feed (B).

[figures omitted; refer to PDF]

The shaft spinning speed at zero natural frequency is called the critical speed, which corresponds to the instability. When the spinning speed exceeds the critical speed, the natural frequency increases with the rotating speed. It also can be seen that both natural frequencies and critical spinning speed of the shaft with the clamped boundary condition are higher than those with simply supported and cantilever boundary conditions.

The numerical results of the Figure 8 are given in Table 2 to show the influence of the gyroscopic effect.

Table 2

The first natural frequencies of the shaft for different boundary conditions and various rotating speed ($\theta ={60}^{°}$).

Rotating  Speed $\Omega$ [rpm]	Frequency ω [Hz]
	1B(C-C)	1F(C-C)	1B(S-S)	1F(S-S)	1B(C-F)	1F(C-F)
0	289.3935	289.3935	207.0878	207.0878	49.0284	49.0284
100	287.7097	291.0775	205.4020	208.7737	47.3530	50.7041
200	286.0260	292.7615	203.7163	210.4598	45.6779	52.3801
300	284.3423	294.4457	202.0308	212.1461	44.0031	54.0564
400	282.6588	296.1299	200.3455	213.8326	42.3286	55.7330
500	280.9753	297.8142	198.6603	215.5191	40.6544	57.4099
600	279.2920	299.4987	196.9752	217.2059	38.9806	59.0871
700	277.6087	301.1832	195.2904	218.8928	37.3070	60.7646
800	275.9256	302.8678	193.6057	220.5798	35.6337	62.4424
900	274.2425	304.5526	191.9211	222.2670	33.9607	64.1205
1000	272.5596	306.2374	190.2367	223.9544	32.2881	65.7990

Figure 9 reveals the first three natural frequencies of the rotating composite thin-walled shaft with different boundary conditions versus the ply angle when rotation speed is 800 rpm. Figure 9(a) represents the clamped boundary condition. As shown in the figure, the maximum values of the first- and second-order natural frequencies occur at θ=70°, while the maximum values of the third order natural frequencies occur at θ=65°. Figure 9(b) represents the simply supported shaft. The figure shows the maximum values of the first order natural frequencies occur at θ=75°, while the second- and third-order natural frequencies occur at θ=70°. Figure 9(c) represents the cantilever composite shaft. The figure indicates that the maximum values of the first-order natural frequencies occur at θ=80°, and the maximum values of the second-order ones occur at θ=75°, while those of the third order occur at θ=70°. It also can be seen from Figure 9 that the higher the mode, the lager the variation amplitude for the three aforementioned boundary conditions. The numerical results of the figure are given in Table 3 to show the influence of the boundary conditions and ply angles.

Table 3

The first three natural frequencies of the shaft for different boundary conditions and ply angles ($\Omega$ = 800rpm).

Boundary conditions	Ply angle θ $[$°$]$	Frequency ω [Hz]
		1B	1F	2B	2F	3B	3F
C-C	0°	89.5956	116.5424	250.1030	277.5957	462.3860	490.3354
	30°	109.5489	136.5861	325.5890	353.6101	653.8300	683.2888
	60°	275.9256	302.8678	724.3444	751.8162	1314.0023	1341.9065
	75°	356.2460	382.9658	765.3613	792.1977	1233.6654	1260.4888
S-S	0°	60.1498	87.1267	211.7786	239.3596	423.3616	451.4875
	30°	71.8321	98.8569	263.8645	291.8418	569.8837	599.3474
	60°	193.2056	220.5798	618.5663	646.1292	1207.7548	1235.8342
	75°	288.0808	314.8729	730.6617	757.5166	1214.6430	1241.4759
C-F	0°	4.0188	30.8277	91.7179	119.2882	266.1023	294.4848
	30°	6.1303	32.9436	109.4554	137.1525	334.9323	364.0539
	60°	35.6337	62.4424	282.7588	310.3222	771.1312	799.4811
	75°	75.0446	101.8206	402.6232	429.7342	895.0152	922.1036

[figures omitted; refer to PDF]

Figure 10 shows the first three natural frequencies of the rotating composite thin-walled shaft versus the ratio of mean radius to thickness. The natural frequencies increase with the ratio of radius to thickness, and this tendency becomes more obvious when the mode is higher.

[figures omitted; refer to PDF]

Figure 11 shows the first three frequencies of the rotating composite thin-walled shaft versus the ratio of length to mean radius. The figure demonstrates that the natural frequencies decrease with the increase of the ratio of length to mean radius. In addition, the higher the mode, the faster the natural frequencies decrease.

[figures omitted; refer to PDF]

5. Conclusions

A refined VAM and Hamilton’s principle were used to establish the vibration differential equation of the rotating composite thin-walled shaft, and the equation was solved by the GDQM. The computational results by this method were compared with those available in literature, and it was validated that the GDQM is accurate and efficient for the frequency analysis of the rotating composite shaft. The influences of the boundary conditions, rotating speed, ply angle, ratio of radius over thickness, and ratio of length over radius on the frequency characteristics were discussed.

Compared with the numerical solution methods such as the Galerkin method and finite element method, the GDQM has the advantages of simple mathematical principle, fast convergence speed, high calculation accuracy, small calculation amount and less memory demand, etc. According to the results obtained in this work, the main conclusions are followed as below:

(1) When the composite shaft rotates, a bifurcation of natural frequencies appears due to the effects of gyroscopic forces. With the increase of spinning speed, the forward whirling frequencies increase, while the backward whirling frequencies decrease.

(2) For the same ply angle, with the clamped boundary condition, the composite thin-walled shaft has higher natural frequencies and critical spinning speed; on the other hand, with simply supported and cantilever boundary conditions, both natural frequencies and critical spinning speed are lower.

(3) The natural frequencies of the rotating composite thin-walled shaft increase with the ratio of diameter over thickness and decrease with the increasing ratio of length to diameter. The natural frequencies of higher-order mode vary faster than those of lower order mode.

Conflicts of Interest

The authors declare that there are no conflicts of interest.

Acknowledgments

Tis work was supported by the National Natural Science Foundation of China (Grants no. 11672166).

References

[1] H. Zinberg, M. F. Symonds, "The development of an advanced composite tail rotor driveshaft," Proceedings of the 26th Annual Forum of the American Helicopter Society, .

[2] H. L. M. dos Reis, R. B. Goldman, P. H. Verstrate, "Thin-walled laminated composite cylindrical tubes, part III: critical speed analysis," Journal of Composites, Technology and Research, vol. 9 no. 2, pp. 58-62, DOI: 10.1520/CTR10430J, 1987.

[3] C. Kim, C. W. Bert, "Critical speed analysis of laminated composite, hollow drive shafts," Composites Part B: Engineering, vol. 3 no. 7-8, pp. 633-643, DOI: 10.1016/0961-9526(93)90087-Z, 1993.

[4] C. W. Bert, C.-D. Kim, "Dynamic instability of composite-material drive shaft subject to fluctuating torque and/or rotational speed," Dynamics and Stability of Systems, vol. 10 no. 2, pp. 125-147, DOI: 10.1080/02681119508806199, 1995.

[5] S. P. Singh, K. Gupta, "Composite shaft rotordynamic analysis using a layerwise theory," Journal of Sound and Vibration, vol. 191 no. 5, pp. 739-756, DOI: 10.1006/jsvi.1996.0153, 1996.

[6] L.-W. Chen, W.-K. Peng, "Dynamic stability of rotating composite shafts under periodic axial compressive loads," Journal of Sound and Vibration, vol. 212 no. 2, pp. 215-230, DOI: 10.1006/jsvi.1997.1405, 1998.

[7] L.-W. Chen, W.-K. Peng, "The stability behavior of rotating composite shafts under axial compressive loads," Composite Structures, vol. 41 no. 3-4, pp. 253-263, DOI: 10.1016/S0263-8223(98)00022-1, 1998.

[8] M. Y. Chang, J. K. Chen, C. Y. Chang, "A simple spinning laminated composite shaft model," International Journal of Solids and Structures, vol. 41 no. 3-4, pp. 637-662, DOI: 10.1016/j.ijsolstr.2003.09.043, 2004.

[9] C.-Y. Chang, M.-Y. Chang, J. H. Huang, "Vibration analysis of rotating composite shafts containing randomly oriented reinforcements," Composite Structures, vol. 63 no. 1, pp. 21-32, DOI: 10.1016/S0263-8223(03)00121-1, 2004.

[10] H. B. H. Gubran, K. Gupta, "The effect of stacking sequence and coupling mechanisms on the natural frequencies of composite shafts," Journal of Sound and Vibration, vol. 282 no. 1-2, pp. 231-248, DOI: 10.1016/j.jsv.2004.02.022, 2005.

[11] O. Song, N. Jeong, L. Librescu, "Implication of conservative and gyroscopic forces on vibration and stability of an elastically tailored rotating shaft modeled as a composite thin-walled beam," The Journal of the Acoustical Society of America, vol. 109 no. 3, pp. 972-981, DOI: 10.1121/1.1348301, 2001.

[12] R. Sino, T. N. Baranger, E. Chatelet, G. Jacquet, "Dynamic analysis of a rotating composite shaft," Composites Science and Technology, vol. 68 no. 2, pp. 337-345, DOI: 10.1016/j.compscitech.2007.06.019, 2008.

[13] A. Boukhalfa, A. Hadjoui, S. M. Hamza Cherif, "Free vibration analysis of a rotating composite shaft using the p-version of the finite element method," International Journal of Rotating Machinery, vol. 2008, 2008.

[14] A. Boukhalfa, A. Hadjoui, "Free vibration analysis of an embarked rotating composite shaft using the hp-version of the FEM," Latin American Journal of Solids and Structures, vol. 7 no. 2, pp. 105-141, DOI: 10.1590/S1679-78252010000200002, 2010.

[15] V. Alwan, A. Gupta, A. S. Sekhar, R. Velmurugan, "Dynamic analysis of shafts of composite materials," Journal of Reinforced Plastics and Composites, vol. 29 no. 22, pp. 3364-3379, DOI: 10.1177/0731684410371404, 2010.

[16] A. Boukhalfa, "Dynamic analysis of a spinning functionally graded material shaft by the p - version of the finite element method," Latin American Journal of Solids and Structures, vol. 11 no. 11, pp. 2018-2038, DOI: 10.1590/S1679-78252014001100007, 2014.

[17] S. Ben Arab, J. Dias Rodrigues, S. Bouaziz, M. Haddar, "A finite element based on Equivalent Single Layer Theory for rotating composite shafts dynamic analysis," Composite Structures, vol. 178, pp. 135-144, DOI: 10.1016/j.compstruct.2017.06.052, 2017.

[18] S. Ben Arab, J. D. Rodrigues, S. Bouaziz, M. Haddar, "Dynamic analysis of laminated rotors using a layerwise theory," Composite Structures, vol. 182, pp. 335-345, DOI: 10.1016/j.compstruct.2017.09.033, 2017.

[19] W. Kim, A. Argento, R. A. Scott, "Free vibration of a rotating tapered composite Timoshenko shaft," Journal of Sound and Vibration, vol. 226 no. 1, pp. 125-147, DOI: 10.1006/jsvi.1999.2289, 1999.

[20] S.-Y. Oh, L. Librescu, O. Song, "Vibration and instability of functionally graded circular cylindrical spinning thin-walled beams," Journal of Sound and Vibration, vol. 285 no. 4-5, pp. 1071-1091, DOI: 10.1016/j.jsv.2004.09.018, 2005.

[21] S. Na, H. Yoon, L. Librescu, "Effect of taper ratio on vibration and stability of a composite thin-walled spinning shaft," Thin-Walled Structures, vol. 44 no. 3, pp. 362-371, DOI: 10.1016/j.tws.2006.02.007, 2006.

[22] Y. S. Ren, Q. Y. Dai, X. Q. Zhang, "Modeling and dynamic analysis of rotating composite shaft," Journal Vibroengineering, vol. 15 no. 4, pp. 1790-1806, 2013.

[23] Y. S. Ren, Y. H. Zhang, X. Q. Zhang, "Vibration and stability of internally damped rotating composite Timoshenko shaft," Journal of Vibroengineering, vol. 17 no. 8, pp. 4404-4420, 2015.

[24] R. Bellman, J. Casti, "Differential quadrature and long-term integration," Journal of Mathematical Analysis and Applications, vol. 34 no. 2, pp. 235-238, DOI: 10.1016/0022-247x(71)90110-7, 1971.

[25] J. R. Vinson, R. L. Sierakowski, The Behavior of Structures Composed of Composite Materials, 2004.

[26] S. Dharmarajan, H. Mccutchen, "Shear Coefficients for Orthotropic Beams," Journal of Composite Materials, vol. 7 no. 4, pp. 530-535, DOI: 10.1177/002199837300700411, 1973.

[27] L. Hua, K. Y. Lam, "Frequency characteristics of a thin rotating cylindrical shell using the generalized differential quadrature method," International Journal of Mechanical Sciences, vol. 40 no. 5, pp. 443-459, DOI: 10.1016/S0020-7403(97)00057-X, 1998.

[28] T. Y. Ng, H. Li, K. Y. Lam, "Generalized differential quadrature for free vibration of rotating composite laminated conical shell with various boundary conditions," International Journal of Mechanical Sciences, vol. 45 no. 3, pp. 567-587, DOI: 10.1016/S0020-7403(03)00042-0, 2003.

[29] J. R. Banerjee, H. Su, "Development of a dynamic stiffness matrix for free vibration analysis of spinning beams," Computers and Structures, vol. 82 no. 23, pp. 2189-2197, 2004.

[30] C. W. Bert, C. D. Kim, "Whirling of composite-material driveshafts including bending–twisting coupling and transverse shear deformation," Journal of Vibration and Acoustics, vol. 117 no. 1, pp. 17-21, DOI: 10.1115/1.2873861, 1995.