Nonlinear free vibrations analysis of rotating shaft based on the Timoshenko beam theory with stretching nonlinearity

Abstract

The nonlinear axial–lateral free vibration behavior of the rotor is investigated by a new sight to the Timoshenko beam theory. In the modeling of the system based on this new nonlinear dynamic model in which the nonlinearities are originated from the stretching of beam centerline, the effects of rotary inertia, gyroscopic forces, and shear deformations are taken into account but the torsional deformation is neglected. According to this new methodology in which the curved geometry of the beam is represented by the Euler angles and the shear deformations, a nonlinear system of coupled variable-coefficient differential equations is derived which is solved by the method of multiple scales in order to determine the nonlinear natural frequencies. The accuracy of the solution method is inspected by comparing the numerical results with the results of two other beam theories and also by the comparison of the free vibration response obtained by the perturbation technique with the numerical integration of the governing equations. The effect of various parameters such as the rotational speed, slenderness ratio, and shear-to-bending stiffness ratio on the free vibrations of the system is inspected. The study demonstrates the nonlinear stretching effect on the free vibration behavior of the shear deformable spinning beam.

Introduction

Rotating machineries have many applications in various mechanical equipments such as compressors, pumps, fans, industrial turbines, aircraft gas turbine engines, spindles of machining tools, helicopter rotors, and motors. The correct analysis and prediction of the dynamic behavior of rotating machines are of great importance because any imperfection of the rotating components can lead to irreparable damages to the systems. Consequently, the rotordynamics is of great importance in engineering mechanics so that the results of the vibration analysis of rotor systems are essential to the design of industrial machineries. The stability and vibration behavior of rotating machineries have been extensively studied by many researchers.^1–6

Generally, several parameters such as the spin speed of the rotor, the mass imbalance of the shaft, the support characteristics of the system, and the properties of the foundation can affect the dynamic characteristics of the rotating systems. Consequently, various models are developed in the previous studies on this field in order to inspect the dynamic behavior of the rotor-bearing systems based on the linear or nonlinear and discrete or continuous mathematical modeling. Linear models which are extensively discussed in the rotordynamics^7–9 do not have the capability of the exact behavior prediction of the complicated systems. However, the nonlinear phenomena can be seen abundantly in rotordynamics. Although the nonlinear analysis of the rotating systems is far more difficult in comparison to its linear counterpart, the need for accurate results has increased the importance of considering the nonlinearity in the dynamic analysis of rotating equipment. In order to model nonlinearities in rotordynamics, various models have been considered in the aforementioned works, such as very simple lumped systems^10,11 and continuous systems in Shaw and Shaw ¹² and Leinonen. ¹³ Nonlinearities in rotor systems can be naturally due to the support conditions,^14–17 inertial or material reasons,^18–21 or geometrical reasons.^17,22–28 The geometric nonlinearity may be caused by nonlinear stretching or large curvatures or nonlinear relationship between the strain and displacement field. Nonlinear inertial effects are caused by the inertial effects of concentrated or distributed masses, and material nonlinearity occurs in the cases that the relationship between the stresses and the strains is nonlinear. Accordingly, lots of research works have been devoted to the nonlinear features of the rotordynamics and it still attracts the attention of many researchers. Ji and Zu ¹⁶ used the perturbation method to investigate the vibration behavior of rotating assemblies based on the linear Timoshenko beam under the effect of a nonlinear bearing pedestal model with nonlinear spring and linear damping characteristics. Cveticanin ²¹ studied the nonlinear vibrations of a rotating shaft which was originated from the elasticity and geometry using the Galerkin method. A geometrically nonlinear model was developed by Luczko ²² based on the Von-Karman nonlinearity assumptions in addition to shear deformations and gyroscopic effects. JianPing et al. ¹⁷ employing the finite element method analyzed the nonlinear dynamics of a rotating assembly based on a continuum model. They also examined the simple discrete model and found the significant differences between these two models. Rubel et al. ²³ presented a model for the dynamics of fast rotating, elastic beams supported in hydrodynamic bearings. Special focus was put on the influence of the nonlinear bearing reaction forces on the dynamics. Their continuous rotor was modeled using Euler–Bernoulli beam theory under the inclusion of rotatory inertia and gyroscopic effects. The finite element method was applied to discretize the beam equation, and the bearing forces were included into the model as point forces in the bearing nodes. Shad et al. ²⁴ investigated the dynamics of the rotor system analytically and numerically by considering nonlinearity due to higher order large deformations in bending. Ding and Zhang ²⁵ studied an isotropic flexible shaft, acted by nonlinear fluid-induced forces generated from oil-lubricated journal bearings and hydrodynamic seal using the Galerkin method. Hosseini and Khadem ²⁶ investigated the free vibrational behavior of a rotating shaft under the in-extensionality condition and nonlinear curvature and inertial effects. Neglecting the shear deformation, they considered the rotary inertia and gyroscopic effects and used the multiple scales method to study the free vibrations of the shaft. They found that the backward and forward frequencies appear simultaneously in the nonlinear natural frequencies. Khadem et al. ²⁷ studied large amplitude vibrations of an in-extensional simply supported rotor. They considered the effect of gyroscopic forces and rotary inertia and neglected the shear deformation. Using the multiple scales method, the effects of physical and geometrical parameters were investigated on the forced vibrations of the rotor. It was shown that in this case, only forward frequencies are excited and the stiffness of the system is of hardening type. Hosseini et al. ²⁸ studied primary resonances of a rotating shaft with nonlinearities in curvature and inertia and in-extensionality condition. They investigated large amplitude free vibrations of the spinning beam under six types of support conditions. They neglected the shear deformation, but took into account the gyroscopic effects and rotary inertia and showed that the forward modes are excited in the primary resonances and both forward and backward modes are excited in the free vibrations of the system. Ishida et al. ²⁹ studied the nonlinear forced vibrations of the spinning distributed mass rotor in which the nonlinear geometries originated from the stretching of the rotor centerline. Applying the method of multiple scales, Hosseini and Khadem ³⁰ studied the primary resonances of a rotor in which stretching nonlinearity, gyroscopic effects, and rotary inertia were considered but the shear deformation was ignored. Shahgholi and Khadem ³¹ investigated the large amplitude primary and parametric resonances of an asymmetrical rotating shaft with simply supported boundary conditions and stretching nonlinearities. They included the rotary inertia and gyroscopic effects in their model and neglect the shear deformation. The asymmetrical effect included the inequality of bending stiffness and mass moments of inertia in the principal axes directions, and the parametric excitation due to inequality between mass moments of inertia and bending stiffness in the direction of principal axes was considered. Hosseini and Zamanian ³² investigated the free vibration of rotor with simply supported boundary conditions under the effect of stretching nonlinearity. They included the gyroscopic force effects and rotary inertia in the study and neglected the shear deformations. It was shown that forward and backward whirling modes are excited simultaneously.

As known by author, none of the previous studies has dealt with the free vibrations analysis of the Timoshenko rotor undergoing nonlinearities originated from the extension of the rotor centerline. The problem is of importance because several previous works reported the occurrence of nonlinear phenomena in rotordynamics, and there are many rotating systems that the Euler–Bernoulli beam theory is not exact enough to predict their correct vibrational behavior.

The stretching nonlinear effects on the dynamic behavior of the rotating systems are of high importance in rotordynamics. This effect is more intensive for moderately thick rotating shafts in which the axial, bending, and shear stiffness is not large enough. Therefore, it is of high interest to develop a dynamic model that satisfactorily implemented the stretching nonlinear effects as well as the shear deformations.

In the past literature of the rotordynamics, two dependent variables are assumed as the degrees of freedom in the definition of the displacement field of the rotor based on the Timoshenko beam theory. These variables consist of the lateral displacement of the mid-surface and the total rotation angle cross section of the rotor. For non-rotating or stationary beams, the variables can be replaced by the total rotation angle and the shear angle of the cross section, while no difference is appeared in the governing equations of the Timoshenko beam and two models are equivalent. ³³ It was shown in the previous study of the authors ³⁴ that if the shear angle is employed as the dependent variable instead of the total rotation angle, the governing equations of the Timoshenko rotating shaft change and some variable coefficients appear in the gyroscopic and inertial terms of the equations. This difference is originated from the fact that in the above-mentioned paper, according to the first-order shear deformation or Timoshenko beam theory, it was assumed that the shear deformation is superposed on the beam cross section after elastic deformation of its centerline while the cross section remained normal to the beam centerline. This modification leads to some deviations in the natural frequencies and critical speeds of the rotor when the results are compared with the common model in the literature. The main discriminations were discussed carefully in Mirtalaie and Hajabasi. ³⁴

Following the previous study of the dynamic behavior of the Timoshenko rotating shaft based on this new insight on the degrees of freedom, the authors are interested in the development of the new dynamic model in which the stretching nonlinearities besides the shear deformation effects are implemented. So, this article follows our previous studies devoted to investigate the nonlinear free lateral vibrations of a flexible shear deformable rotor with stretching nonlinearities based on the new insight to the shear deformations of the beam cross section and was constructed as follows.

First, the concepts, assumptions, and geometric principles required for derivation of the mathematical model are represented. Then, employing Hamilton’s principle and based on the new representation of deformation of the cross section of the beam according to the Timoshenko beam theory, the nonlinear governing differential equations of motion as well as the boundary conditions are derived. Afterward, the method of multiple scales is used to examine the free vibration behavior of the rotor system, and the exact values of nonlinear natural frequencies in forward and backward whirling modes are determined for various parameters. The accuracy of the results is evaluated by the numerical integrations, and the effect of various system parameters on the free vibrational behaviors of the system such as rotating speed and slenderness ratio of the shaft is investigated.

Equations of motion

In this section, the nonlinear governing differential equations of motion of the rotor based on the Timoshenko beam theory with the effect of stretching are derived employing Hamilton’s principle. Here, a flexible horizontal rotating shaft in which the rotating speed is constant and the cross section is circular and uniform along the axial direction is assumed. The simply supported boundary conditions at the end-points of the shaft are considered such that both of the supports are fixed in axial displacements as demonstrated in Figure 1. Therefore, no shortening effect can be considered and the axial displacements appeared in the governing equations explicitly. Because the moderately thick shafts are assumed, the shear deformations should be taken into account. Furthermore, the effects of the rotary inertia are considered but the torsional vibrations, energy wasting mechanisms, and structural deficiencies are neglected. Furthermore, it is assumed that the amplitude of the vibration is large and stretching nonlinearity which is originated from the extension of the shaft centerline is considered.

Figure 1.

The configuration of rotor-bearing system and the inertial and moving coordinate system.

As shown in Figure 1, a fixed Cartesian coordinate system $(X,Y,Z)$ is used as the inertial reference frame which describes the original configuration of the shaft. The $(x,y,z)$ is the local orthogonal curvilinear moving coordinate system that rotates with the shaft angular velocity and describes the deformed geometry of the rotating shaft. This coordinate system is defined such that the x-axis is tangent to the deformed centerline of the shaft and the y and z coordinates show the normal and bi-normal vectors of the deformed centerline of the shaft.

The orientation of the inertial frame reference with respect to the local coordinate system is defined by a series of the Euler angles, employing their sequent represents the angles of the $(x,y,z)$ system regarding to $(X,Y,Z)$ . In this study, the sequence of Euler angles is selected as follows. First, the body rotates by the angle $\psi (x,t)$ about the Z-axis, and the updated reference frame is represented by $(x_1,y_1,z_1)$ , and then, it rotates by $\phi (x,t)$ about the y₂-axis so that the $(x_2,y_2,z_2)$ reference frame is obtained. Finally, the updated reference frame is rotated about the x₂-axis by the angle $\theta (x,t)$ in order to construct the local reference frame. This order of rotation is shown schematically in Figure 2. The shear and torsional warpings are not described up to this level of modeling because the plane which was normal to the mid-surface of the shaft in the undeformed geometry is mapped to the plane perpendicular to the curved mid-surface of the shaft in the deformed configuration. Based on the methodology developed in the previous study of the authors, ³⁴ first the curved geometry of the beam centerline after deformation is considered and the local coordinate system is constructed on it, and then, the shear deformations are represented in this local coordinate system. The main strategy of this methodology is the superposition of warpings on the planes defined by the tangent–normal–binormal local coordinate system. ³⁴ Consequently, the angular displacements of the normal plane to the mid-surface of the shaft modify as a result of shear deformation contribution.

Figure 2.

Euler angle sequences.

The displacement components of the elements on the middle surface of the beam are assumed as $u(x,t)$ for axial displacement in the x direction and $[v(x,t),w(x,t)]$ for lateral deflections in $(y,z)$ directions, respectively. Also, the shear deformation components are considered as $({\gamma }_y,{\gamma }_z)$ around the $(y,z)$ directions, respectively. In this study, the torsional warpings are neglected but the shear deformations due to lateral vibrations are taken into account, and it is assumed that the shear strain has uniform contribution through the cross section of the beam. Therefore, according to the Timoshenko beam theory, the displacement field is defined as ³³

$$\begin{array}{l}{u}_1=u_1(x,y,z,t)=u(x,t)+y[-{\theta }_3(x,t)+{\gamma }_z(x,t)]\\ +z[{\theta }_2(x,t)+{\gamma }_y(x,t)]\end{array}$$

(1-a)

$$u_2=u_2(x,y,z,t)=v(x,t)-z{\theta }_1(x,t)$$

(1-b)

$$u_3=u_3(x,y,z,t)=w(x,t)+y{\theta }_1(x,t)$$

(1-c)

in which $(u_1,u_2,u_3)$ stand for the displacements of a sample point on beam cross section in the $(x,y,z)$ directions, respectively. Moreover, $({\theta }_1,{\theta }_2,{\theta }_3)$ demonstrate the cross-sectional angles of rotation around the x-, y-, and z-axes, respectively. Also, $({\gamma }_y,{\gamma }_z)$ show the shear rotation angles of the cross section at the reference point around the (y, z)-axis, respectively. The geometrical representation of the cross-sectional lateral and angular displacements and the shear angles in y-x and z-x planes has been schematically shown in Figure 3. It can be seen that under the effect of transverse shear forces, the normal plane to the beam middle surface rotates more and the shape of the elements in the lateral view changes from dashed line rectangles to solid line trapezoids accordingly.

Figure 3.

Displacements and shear deformation components of the sample element.

By the definition, the non-zero components of linear strains are obtained as

$$\begin{array}{l}{\epsilon }_{11}=e+y\left[-\frac{\partial }{\partial x}{\theta }_3(x,t)+\frac{\partial }{\partial x}{\gamma }_z(x,t)\right]\\ +z\left[\frac{\partial }{\partial x}{\theta }_2(x,t)+\frac{\partial }{\partial x}{\gamma }_y(x,t)\right]\end{array}$$

(2-a)

$${\epsilon }_{12}=\frac{\partial }{\partial x}v(x,t)-z\frac{\partial }{\partial x}{\theta }_1(x,t)$$

(2-b)

$${\epsilon }_{13}=\frac{\partial }{\partial x}w(x,t)+y\frac{\partial }{\partial x}{\theta }_1(x,t)$$

(2-c)

where e is the extensional or axial strain along the reference line of the shaft which is defined as ³³

$$e=\sqrt{{[1+\frac{\partial u(x,t)}{\partial x}]}^2+{\left(\frac{\partial v(x,t)}{\partial x}\right)}^2+{\left(\frac{\partial w(x,t)}{\partial x}\right)}^2}-1$$

(3)

Keeping in mind the concept of local displacement and the geometrical representation of curvature which was completely explained in Nayfeh and Pai ³³ and Bower, ³⁵ the relation between the derivatives of rotation angles, curvature of the cross section, and Euler angles is represented as

$$\frac{\partial {\theta }_1(x,t)}{\partial x}={\rho }_1=-\frac{\partial }{\partial x}\psi (x,t)\sin \phi (x,t)+\frac{\partial }{\partial x}\theta (x,t)$$

(4-a)

$$\begin{array}{l}\frac{\partial {\theta }_2(x,t)}{\partial x}={\rho }_2=\frac{\partial }{\partial x}\psi (x,t)\sin \theta (x,t)\cos \phi (x,t)\\ +\frac{\partial }{\partial x}\phi (x,t)\cos \theta (x,t)\end{array}$$

(4-b)

$$\begin{array}{l}\frac{\partial {\theta }_3(x,t)}{\partial x}={\rho }_3=\frac{\partial }{\partial x}\psi (x,t)\cos \theta (x,t)\cos \phi (x,t)\\ -\frac{\partial }{\partial x}\phi (x,t)\sin \theta (x,t)\end{array}$$

(4-c)

As the torsional vibrations are neglected in this study, it can be concluded that $\theta (x,t)=\Omega t$ , in which $\Omega$ shows the spin speed of the shaft. Besides, the relations between the Euler angles and the derivatives of displacement components are represented as follows ³⁵

$$\begin{array}{l}\sin \varphi (x,t)=-\frac{1}{1+e}\frac{\partial w(x,t)}{\partial x},\\ \cos \varphi (x,t)=\frac{1}{1+e}\sqrt{{\left[1+\frac{\partial u(x,t)}{\partial x}\right]}^2+{\left[\frac{\partial v(x,t)}{\partial x}\right]}^2}\end{array}$$

(5-a)

$$\begin{array}{l}\sin \psi (x,t)=\frac{\frac{\partial v(x,t)}{\partial x}}{\sqrt{{\left[1+\frac{\partial u(x,t)}{\partial x}\right]}^2+{\left[\frac{\partial v(x,t)}{\partial x}\right]}^2}}\\ \cos \psi (x,t)=\frac{1+\frac{\partial u(x,t)}{\partial x}}{\sqrt{{\left[1+\frac{\partial u(x,t)}{\partial x}\right]}^2+{\left[\frac{\partial v(x,t)}{\partial x}\right]}^2}}\end{array}$$

(5-b)

which are employed in the following for the evaluation of potential and kinetic energy expressions because Hamilton’s principle which states that ${\displaystyle {\int }_{t_1}^{t_2}\delta (\mathrm{\Pi }-T)dt}=0$ is employed in this study to derive the governing equations. The potential or strain energy of the shaft is expressed as

$$\mathrm{\Pi }=\frac{1}{2}\underset{A}{\int }{\int }_0^L\{E({\epsilon }_{11}^2+{\epsilon }_{22}^2+{\epsilon }_{33}^2)+G({\epsilon }_{12}^2+{\epsilon }_{23}^2+{\epsilon }_{13}^2)\}\mathrm{dx}\mathrm{dA}$$

(6)

where E and G show Young’s modulus of elasticity and shear modulus, respectively. Substituting the strain components from equations (2), using equations (3)–(5), and expanding the results into a Taylor series by retaining terms up to the order of four, the strain energy of the shaft is simplified as

$$\mathrm{\Pi }=\underset{0}{\overset{L}{\int }}\left\{\begin{array}{c}\mathrm{EI}\left[\begin{array}{c}-\frac{\partial u}{\partial x}{\left(\frac{{\partial }^{2}v}{\partial x^2}\right)}^2+\frac{\partial u}{\partial x}\frac{{\partial }^{2}v}{\partial x^2}\frac{\partial {\gamma }_z}{\partial x}+\frac{\partial u}{\partial x}\frac{{\partial }^{2}w}{\partial x^2}\frac{\partial {\gamma }_y}{\partial x}-\frac{{\partial }^{2}w}{\partial x^2}\frac{\partial {\gamma }_y}{\partial x}-\frac{{\partial }^{2}v}{\partial x^2}\frac{\partial {\gamma }_z}{\partial x}\hfill \\ \multicolumn{1}{c}{+\frac{1}{2}{\left(\frac{{\partial }^{2}v}{\partial x^2}\right)}^2+\frac{1}{2}{\left(\frac{{\partial }^{2}w}{\partial x^2}\right)}^2-\frac{\partial u}{\partial x}{\left(\frac{{\partial }^{2}w}{\partial x^2}\right)}^2-\frac{{\partial }^{2}v}{\partial x^2}\frac{{\partial }^{2}u}{\partial x^2}\frac{\partial v}{\partial x}+\frac{{\partial }^{2}u}{\partial x^2}\frac{\partial {\gamma }_z}{\partial x}\frac{\partial v}{\partial x}-\frac{{\partial }^{2}w}{\partial x^2}\frac{{\partial }^{2}u}{\partial x^2}\frac{\partial w}{\partial x}}\\ \multicolumn{1}{c}{+\frac{{\partial }^{2}u}{\partial x^2}\frac{\partial {\gamma }_y}{\partial x}\frac{\partial w}{\partial x}+\frac{1}{2}{\left(\frac{\partial {\gamma }_y}{\partial x}\right)}^2+\frac{1}{2}{\left(\frac{\partial {\gamma }_z}{\partial x}\right)}^2}\end{array}\right]\hfill \\ \multicolumn{1}{c}{+\frac{1}{2}\mathrm{EA}[{\left(\frac{\partial u}{\partial x}\right)}^2+\frac{\partial u}{\partial x}{\left(\frac{\partial v}{\partial x}\right)}^2+\frac{\partial u}{\partial x}{\left(\frac{\partial w}{\partial x}\right)}^2]+\frac{1}{2}\kappa GA({\gamma }_y(x,t{)}^2+{\gamma }_z(x,t{)}^2)}\end{array}\right\}\mathrm{dx}$$

(7)

where A stands for the cross-sectional area and $\kappa$ shows the shear correction factor in the Timoshenko beam theory. In the above relation, because of symmetry of the cross section, the terms ${\int }_{A}yz\mathrm{dA},{\int }_{A}y\mathrm{dA}$ and ${\int }_{A}z\mathrm{dA}$ vanish. Moreover, ${\int }_A{y}^2\mathrm{dA}={\int }_A{z}^2\mathrm{dA}=I$ and ${\int }_A(y^2+z^2)\mathrm{dA}=J$ demonstrate the transverse and polar area moments of inertia of the cross section, respectively.

The kinetic energy of the shaft is defined by

$$T=\frac{1}{2}\underset{A}{\int }{\int }_0^L\rho \{{\left(\frac{\partial u_1}{\partial t}\right)}^2+{\left(\frac{\partial u_2}{\partial t}\right)}^2+{\left(\frac{\partial u_3}{\partial t}\right)}^2\}\mathrm{dx}\mathrm{dA}$$

(8)

in which $\rho$ is the mass density of the shaft. In the evaluation of the time derivatives of the displacement components using the displacement field (1), similar to the curvature components of the shaft, the components of the angular velocity vector of the moving frame are defined as ³³

$$\frac{\partial {\theta }_1(x,t)}{\partial t}={\omega }_1=\Omega -\frac{\partial }{\partial t}\psi (x,t)\sin \phi (x,t)+\frac{\partial }{\partial t}\theta (x,t)$$

(9-a)

$$\begin{array}{l}\frac{\partial {\theta }_2(x,t)}{\partial t}={\omega }_2=\frac{\partial }{\partial t}\psi (x,t)\sin (\theta (x,t))\cos \phi (x,t)\\ +\frac{\partial }{\partial t}\phi (x,t)\cos (\theta (x,t))\end{array}$$

(9-b)

$$\begin{array}{l}\frac{\partial {\theta }_3(x,t)}{\partial t}={\omega }_3=\frac{\partial }{\partial t}\psi (x,t)\cos (\theta (x,t))\cos \phi (x,t)\\ -\frac{\partial }{\partial t}\phi (x,t)\sin (\theta (x,t))\end{array}$$

(9-c)

Substituting these velocity components, using equations (3)–(5), and expanding the outcomes into Taylor series by retaining terms up to the order of four, the kinetic energy of the shaft is simplified as

$$\begin{array}{l}T={\displaystyle \underset{0}{\overset{L}{\int }}\{\frac{1}{2}{I}_m[\left(\frac{{\partial }^{2}v}{\partial x\partial t}\right)+\left(\frac{{\partial }^{2}w}{\partial x\partial t}\right)+{\left(\frac{\partial {\gamma }_y}{\partial t}\right)}^2+{\left(\frac{\partial {\gamma }_z}{\partial t}\right)}^2-2\frac{\partial u}{\partial x}{\left(\frac{{\partial }^{2}v}{\partial x\partial t}\right)}^2-2\frac{\partial u}{\partial x}{\left(\frac{{\partial }^{2}w}{\partial x\partial t}\right)}^2}\\ -2\frac{{\partial }^{2}u}{\partial x\partial t}\frac{{\partial }^{2}v}{\partial x\partial t}\frac{\partial v}{\partial x}-2\frac{{\partial }^{2}u}{\partial x\partial t}\frac{{\partial }^{2}w}{\partial x\partial t}\frac{\partial w}{\partial x}]+I_m[\sin (\Omega t)(-\frac{\partial {\gamma }_z}{\partial t}\frac{{\partial }^{2}w}{\partial x\partial t}+\frac{\partial {\gamma }_y}{\partial t}\frac{{\partial }^{2}v}{\partial x\partial t}\\ +\frac{\partial {\gamma }_z}{\partial t}\frac{{\partial }^{2}w}{\partial x\partial t}\frac{\partial u}{\partial x}-\frac{\partial {\gamma }_y}{\partial t}\frac{{\partial }^{2}v}{\partial x\partial t}\frac{\partial u}{\partial x}-\frac{\partial {\gamma }_y}{\partial t}\frac{{\partial }^{2}u}{\partial x\partial t}\frac{\partial v}{\partial x}+\frac{\partial {\gamma }_z}{\partial t}\frac{{\partial }^{2}u}{\partial x\partial t}\frac{\partial w}{\partial x})]\\ +I_m[\cos (\Omega t)(-\frac{\partial {\gamma }_y}{\partial t}\frac{{\partial }^{2}w}{\partial x\partial t}-\frac{\partial {\gamma }_z}{\partial t}\frac{{\partial }^{2}v}{\partial x\partial t}+\frac{\partial {\gamma }_z}{\partial t}\frac{{\partial }^{2}u}{\partial x\partial t}\frac{\partial v}{\partial x}+\frac{\partial {\gamma }_z}{\partial t}\frac{{\partial }^{2}v}{\partial x\partial t}\frac{\partial u}{\partial x}\\ +\frac{\partial {\gamma }_y}{\partial t}\frac{{\partial }^{2}w}{\partial x\partial t}\frac{\partial u}{\partial x}+\frac{\partial {\gamma }_y}{\partial t}\frac{{\partial }^{2}u}{\partial x\partial t}\frac{\partial w}{\partial x})]+\frac{1}{2}\rho A\left[{\left(\frac{\partial u}{\partial t}\right)}^2+{\left(\frac{\partial v}{\partial t}\right)}^2+{\left(\frac{\partial w}{\partial t}\right)}^2\right]+J_m\left[\Omega \frac{{\partial }^{2}v}{\partial x\partial t}\frac{\partial w}{\partial x}+\frac{1}{2}{\Omega }^2-2\Omega \frac{{\partial }^{2}v}{\partial x\partial t}\frac{\partial w}{\partial x}\frac{\partial u}{\partial x}-\Omega \frac{{\partial }^{2}u}{\partial x\partial t}\frac{\partial v}{\partial x}\frac{\partial w}{\partial x}\right]\}dx\end{array}$$

(10)

in which ${\int }_A\rho y\mathrm{dA}=0,{\int }_A\rho z\mathrm{dA}=0,{\int }_A\rho \mathrm{yz}\mathrm{dA}=0$ as a result of symmetry of the cross section. Also, ${\int }_A\rho y^2\mathrm{dA}={\int }_A\rho z^2\mathrm{dA}=I_m$ and ${\int }_A\rho (y^2+z^2)\mathrm{dA}=J_m$ are the diametrical and polar mass moments of inertia per length of the rotor, respectively. The variable (periodic) coefficient terms appeared in the kinetic energy expression lead to variable-coefficient terms in the governing differential equations which are the essential distinction of the new modeling methodology and the common models of the rotor based on Timoshenko beam theory. The origin of this distinction and other aspects of the modeling have been inspected carefully in the previous study of the authors. ³⁴

Inserting the potential and kinetic energy expressions in Hamilton’s principle and applying the variational operator, integrating by parts and doing some simplifications, the governing equations of motion are combined as

$$\begin{array}{c}2u″+U″\overline{U}\prime +U\prime \overline{U}″+\epsilon (-2\stackrel{··}{u}-U\prime \overline{\gamma }‴+U‴\overline{U}″+U″\overline{U}‴-\overline{U}″\gamma ″-U″\gamma ″+\overline{U}\prime U^{(\mathrm{IV})}+U\prime {\overline{U}}^{(\mathrm{IV})}\hfill \\ -\overline{U}\prime \gamma ‴)+{\epsilon }^2[-\overline{U}\prime \stackrel{··}{U}″-U\prime \stackrel{··}{\overline{U}}″-U″\stackrel{··}{\overline{U}}\prime -\overline{U}″\stackrel{··}{U}\prime +2i\Omega (\stackrel{·}{U}\prime \overline{U}″-\stackrel{·}{\overline{U}}″U\prime +\stackrel{·}{U}″\overline{U}\prime -\stackrel{·}{\overline{U}}\prime U″)\hfill \\ +\sin (\Omega \tau )(-iU\prime \stackrel{··}{\overline{\gamma }}\prime +i\overline{U}″\stackrel{··}{\gamma }+i\overline{U}\prime \stackrel{··}{\gamma }\prime -iU″\stackrel{··}{\overline{\gamma }}-\Omega U″\stackrel{·}{\overline{\gamma }}-\Omega \overline{U}\prime \stackrel{·}{\gamma }\prime -\Omega U\prime \stackrel{·}{\overline{\gamma }}\prime -\Omega \overline{U}″\stackrel{·}{\gamma })\hfill \\ +\cos (\Omega \tau )(\overline{U}″\stackrel{··}{\gamma }-i\Omega U\prime \stackrel{·}{\overline{\gamma }}\prime +U\prime \stackrel{··}{\overline{\gamma }}\prime +\overline{U}\prime \stackrel{··}{\gamma }\prime -i\Omega U″\stackrel{·}{\overline{\gamma }}+U″\stackrel{··}{\overline{\gamma }}+i\Omega \overline{U}″\stackrel{·}{\gamma }+i\Omega \overline{U}\prime \stackrel{·}{\gamma }\prime )]=0\hfill \end{array}$$

(11-a)

$$\begin{array}{c}{U}^{(\mathrm{IV})}+\stackrel{··}{U}-\gamma ‴-\epsilon \stackrel{··}{w}″+i\beta \stackrel{·}{U}″-3u‴U″-4u″U‴-2u\prime U^{(\mathrm{IV})}+u″\gamma ″+u\prime \gamma ‴-u^{(\mathrm{IV})}U\prime \hfill \\ -i\beta (2u″\stackrel{·}{U}\prime +2u\prime \stackrel{·}{U}″+\stackrel{·}{u}″U\prime +\stackrel{·}{u}\prime U″)+\epsilon (\stackrel{··}{u}\prime U″+2\stackrel{·}{u}″\stackrel{·}{U}\prime +2u″\stackrel{··}{U}\prime +2u\prime \stackrel{··}{U}″+2\stackrel{·}{u}\prime \stackrel{·}{U}″+\stackrel{··}{u}″U\prime )\hfill \\ -\epsilon \sin (\Omega \tau )[\Omega u\prime \stackrel{·}{\gamma }\prime +\Omega \stackrel{·}{\gamma }\prime +iu″\stackrel{··}{\gamma }+iu\prime \stackrel{··}{\gamma }\prime -i\stackrel{··}{\gamma }\prime -\Omega u″\stackrel{·}{\gamma }]\hfill \\ -\epsilon \cos (\Omega \tau )[u\prime \stackrel{··}{\gamma }\prime -\stackrel{··}{\gamma }\prime +i\Omega u\prime \stackrel{·}{\gamma }\prime +i\Omega u″\stackrel{·}{\gamma }+u″\stackrel{··}{\gamma }-i\Omega \stackrel{·}{\gamma }\prime ]-\frac{1}{\epsilon }(u\prime U″+u″U\prime )=0\hfill \end{array}$$

(11-b)

$$\begin{array}{c}U‴+\epsilon \stackrel{··}{\gamma }-\gamma ″+\lambda \gamma -u\prime U‴-u‴U\prime -2u″U″\hfill \\ -\epsilon \sin (\Omega \tau )[2i\stackrel{·}{u}\prime \stackrel{·}{U}\prime -\Omega \stackrel{·}{U}\prime +\Omega \stackrel{·}{u}\prime U\prime +\Omega u\prime \stackrel{·}{U}\prime +i\stackrel{··}{u}\prime U\prime -i\stackrel{··}{U}\prime +iu\prime \stackrel{··}{U}\prime ]\hfill \\ -\epsilon \cos (\Omega \tau )[-2\stackrel{·}{u}\prime \stackrel{·}{U}\prime -i\Omega \stackrel{·}{U}\prime +i\Omega \stackrel{·}{u}\prime U\prime +i\Omega u\prime \stackrel{·}{U}\prime -\stackrel{··}{u}\prime U\prime +\stackrel{··}{U}\prime -u\prime \stackrel{··}{U}\prime ]=0\hfill \end{array}$$

(11-c)

where the governing system of equations has been represented in nondimensional form employing the following dimensionless quantities

$$\begin{array}{ccc}\xi =\frac{x}{L}\hfill & \tau =\alpha t\hfill & {\Omega }^{*}=\frac{\Omega }{\alpha }\hfill \\ \multicolumn{1}{c}{U(\xi ,\tau )=\frac{U^{*}(x,t)}{L}}& u^{*}(\tau ,\tau )=\frac{u(x,t)}{L}\hfill & \gamma (\xi ,\tau )={\gamma }^{*}(x,t)\hfill \\ \multicolumn{1}{c}{\beta =\frac{\Omega J_m}{\sqrt{\mathrm{EI}\rho A}}}& \lambda =\frac{\kappa \mathrm{AG}{L}^2}{\mathrm{EI}}\hfill & \epsilon =\frac{I_m}{\rho A{L}^2}\hfill \\ \multicolumn{1}{c}{\alpha =\sqrt{\frac{\mathrm{EI}}{\rho A{L}^4}}}\end{array}$$

(12)

in which the following complex variables are used

$$U^{*}(x,t)=v(x,t)+iw(x,t)$$

(13-a)

$${\gamma }^{*}(x,t)={\gamma }_z(x,t)+i{\gamma }_y(x,t),i^2=-1$$

(13-b)

The asterisk for parameter $u^{*}(\tau ,\tau )$ has been dropped for simplicity. The over bar denotes the complex conjugate and the prime and dot stand for differentiation with respect to x and t, respectively. It is seen that the equations of motion are nonlinear, variable coefficients and coupled together. As the special cases, if the nonlinear terms are neglected, the system discussed in Mirtalaie and Hajabasi ³⁴ is obtained; furthermore, if the spin of the rotor is assumed to be zero, the equations of motion of the stationary Timoshenko beam ³³ are obtained which are decoupled. Also, if the nonlinear terms and shear deformations are eliminated, the equations of motion will be coincident with the common Rayleigh rotor model in the open literature. ³ If short bearings are considered as supports, the simply supported boundary conditions can be assumed which are expressed as

$$u^{*}(\xi ,\tau )=0,\xi =0,1$$

(14-a)

$$U(\xi ,\tau )=0,\xi =0,1$$

(14-b)

$$\frac{{\partial }^{2}U(\xi ,\tau )}{\partial {\xi }^2}-\frac{\partial \gamma (\xi ,\tau )}{\partial \xi }=0,\xi =0,1$$

(14-c)

In the following, the above-mentioned equations are solved in order to analyze the dynamic behavior of the system.

Multiple scales solution

Here, the method of multiple scales is employed to examine the nonlinear system of equations derived earlier. According to this method, the axial and lateral nondimensional displacement components and the shear deformation angles are perturbed as^36,37

$$u(\xi ,\tau )={\epsilon }^2{u}_2(\xi ,T_0,T_1)+{\epsilon }^3{u}_3(\xi ,T_0,T_1)$$

(15-a)

$$U(\xi ,\tau )=\epsilon U_1(\xi ,T_0,T_1)+{\epsilon }^2{U}_2(\xi ,T_0,T_1)$$

(15-b)

$$\gamma (\xi ,\tau )=\epsilon {\gamma }_1(\xi ,T_0,T_1)+{\epsilon }^2{\gamma }_2(\xi ,T_0,T_1)$$

(15-c)

where $T_i={\epsilon }^i\tau (i=0,1)$ demonstrate the time scales. The derivatives with respect to nondimensional time in terms of the time scales can be expressed as follows by employing the chain rule

$$\frac{d}{d\tau }=D_0+\epsilon D_1+\cdots$$

(16-a)

$$\begin{array}{l}\frac{d^2}{d{\tau }^2}={(D_0+\epsilon D_1+\mathrm{\dots })}^2\\ =D_0^2+2\epsilon D_0{D}_1+{\epsilon }^2{D}_1^2+\mathrm{\dots }\end{array}$$

(16-b)

in which $D_i=\partial /\partial T(i=0,1)$ . Using the above relations in the nondimensional equations of motion, simplifying the results, and equating the terms that are coefficients of like powers of $D_i=\partial /\partial T(i=0,1)$ read

$$\begin{array}{c}\multicolumn{1}{c}{O(\epsilon ):}\\ \multicolumn{1}{c}{D_0^2{U}_1+U_1^{(\mathrm{IV})}-\gamma ″{\prime }_1-\beta D_0{U″}_1=0}\end{array}$$

(17-a)

$${U‴}_1-{\gamma ″}_1+\lambda {\gamma }_1=0$$

(17-b)

$$\begin{array}{c}\multicolumn{1}{c}{O({\epsilon }^2):}\\ \multicolumn{1}{c}{{u″}_2=-\frac{1}{2}[\overline{U}{\prime }_1{U″}_1+U{\prime }_1{\overline{U}″}_1]}\end{array}$$

(18-a)

$$\begin{array}{c}\multicolumn{1}{c}{D_0^2{U}_2+U_2^{(\mathrm{IV})}-\gamma ″{\prime }_2-\beta D_0{U″}_2=D_0^2{U″}_1}\\ -i\beta D_1{U″}_1-2D_0{D}_1{U}_1+u{\prime }_2{U″}_1+{u″}_{2}U{\prime }_1\hfill \\ +[-D_0^2\gamma {\prime }_1-i\Omega D_0\gamma {\prime }_1]\cos (\Omega T_0)\hfill \\ +[\Omega D_0\gamma {\prime }_1-i{D}_0^2\gamma {\prime }_1]\sin (\Omega T_0)\hfill \end{array}$$

(18-b)

$${U^{‴}}_2-{{\gamma }^{″}}_2+\lambda {\gamma }_2=-D_0^2{\gamma }_1+[D_0^2{U^{\prime }}_1-i\Omega D_0{U^{\prime }}_1]\cos (\Omega T_0)+[-i{D}_0^2{U^{\prime }}_1-\Omega D_0{U^{\prime }}_1]\sin (\Omega T_0)$$

(18-c)

A similar procedure is performed for the boundary conditions and the results are obtained as follows

$$\begin{array}{c}O({\epsilon }^j):j=1,2\hfill \\ u_j(\xi ,\tau )=0,\xi =0,1\hfill \end{array}$$

(19-a)

$$U_j(\xi ,\tau )=0,\xi =0,1$$

(19-b)

$${U″}_j-\gamma {\prime }_j=0,\xi =0,1$$

(19-c)

After perturbation, first the linear equations (17) and the boundary conditions (19) for $j=1$ are attacked. As it has been explained extensively in Mirtalaie and Hajabasi, ³⁴ the solution of the system is obtained as

$$U_1(\xi ,T_0,T_1)=[F_f(T_1)e^{i{\omega }_f{T}_0}+F_b(T_1)e^{i{\omega }_b{T}_0}]\sin (n\pi \xi )$$

(20-a)

$$\begin{array}{l}{\gamma }_1(\xi ,T_0,T_1)=\frac{n^3{\pi }^3}{\lambda +n^2{\pi }^2}\\ [F_f(T_1)e^{i{\omega }_f{T}_0}+F_b(T_1)e^{i{\omega }_b{T}_0}]\cos (n\pi \xi )\end{array}$$

(20-b)

in which $i^2=-1$ , n is the mode number, and $F_f(T_1)$ and $F_b(T_1)$ are functions with complex values and will be determined at next step of perturbation method. The frequencies of forward ${\omega }_f$ and backward ${\omega }_b$ whirling modes are determined as follows ³⁴

$${\omega }_f=\frac{n^2{\pi }^2}{2}[\beta +\sqrt{\frac{{\beta }^2{n}^2{\pi }^2+\lambda {\beta }^2+4\lambda }{\lambda +n^2{\pi }^2}}]$$

(21-a)

$${\omega }_b=\frac{n^2{\pi }^2}{2}[\beta -\sqrt{\frac{{\beta }^2{n}^2{\pi }^2+\lambda {\beta }^2+4\lambda }{\lambda +n^2{\pi }^2}}]$$

(21-b)

Then, the equations obtained in the second scale, namely, $O({\epsilon }^2)$ are examined. Substituting equations (20) and (18-a) into equations (18-b) and (18-c) and using the expansions $\sin (\Omega T_0)=(1/2i)(e^{i\Omega T_0}-e^{-i\Omega T_0})$ and $\cos (\Omega T_0)=(1/2)(e^{i\Omega T_0}+e^{-i\Omega T_0})$ lead to

$$D_0^2{U}_2+U_2^{(IV)}-{{\gamma }^{‴}}_2-\beta D_0{U^{″}}_2=U_f(\xi ,T_1)e^{i{\omega }_f{T}_0}+U_b(\xi ,T_1)e^{i{\omega }_b{T}_0}+CC+NSTs$$

(22-a)

$${U^{‴}}_2-{{\gamma }^{″}}_2+\lambda {\gamma }_2={\mathrm{\Gamma }}_f(\xi ,T_1)e^{i{\omega }_f{T}_0}+{\mathrm{\Gamma }}_b(\xi ,T_1)e^{i{\omega }_b{T}_0}+CC+NSTs$$

(22-b)

in which the abbreviations CC and NSTs stand for complex conjugated and non-secular terms, respectively. The coefficients that produce secular terms are determined as

$$\begin{array}{c}\multicolumn{1}{c}{U_f(\xi ,T_1)=[+(\frac{3}{2}{n}^4{\pi }^4{F}_f{(T_1)}^2{\overline{F}}_f(T_1)+3n^4{\pi }^4{F}_f(T_1)F_b(T_1){\overline{F}}_b(T_1)){\cos }^2(n\pi \xi )}\\ +n^2{\pi }^2{\omega }_f^2{F}_f(T_1)+i(\beta n^2{\pi }^2-2{\omega }_f){\stackrel{·}{F}}_f(T_1)-\frac{1}{2}{n}^4{\pi }^4{F}_f(T_1)F_b(T_1){\overline{F}}_b(T_1)\hfill \\ -\frac{1}{4}{n}^4{\pi }^4{F}_f{(T_1)}^2{\overline{F}}_f(T_1)]\sin (n\pi \xi )\hfill \end{array}$$

(23-a)

$${\mathrm{\Gamma }}_f(\xi ,T_1)=\frac{n^3{\pi }^3{\omega }_f^2\cos (n\pi \xi )F_f(T_1)}{\lambda +n^2{\pi }^2}$$

(23-b)

$$\begin{array}{c}\multicolumn{1}{c}{U_b(\xi ,T_1)=[+(\frac{3}{2}{n}^4{\pi }^4{F}_b{(T_1)}^2{\overline{F}}_b(T_1)+3n^4{\pi }^4{F}_f(T_1)F_b(T_1){\overline{F}}_f(T_1)){\cos }^2(n\pi \xi )}\\ +n^2{\pi }^2{\omega }_b^2{F}_b(T_1)+i(\beta n^2{\pi }^2-2{\omega }_b){\stackrel{·}{F}}_b(T_1)-\frac{1}{2}{n}^4{\pi }^4{F}_f(T_1)F_b(T_1){\overline{F}}_f(T_1)\hfill \\ -\frac{1}{4}{n}^4{\pi }^4{F}_b{(T_1)}^2{\overline{F}}_b(T_1)]\sin (n\pi \xi )\hfill \end{array}$$

(23-c)

$${\mathrm{\Gamma }}_b(\xi ,T_1)=\frac{n^3{\pi }^3{\omega }_b^2\cos (n\pi \xi )F_b(T_1)}{\lambda +n^2{\pi }^2}$$

(23-d)

Following the Nayfeh, ³⁷ the inhomogeneous equations (22) have solution if and only if the solvability condition is fulfilled because the homogeneous part of equations (22) has non-trivial solutions. This condition represents the orthogonality of the right side of these equations to every solution of its adjoint problem. Because the homogeneous parts of equations (22) are self-adjoint, the solvability conditions are obtained as

$$\underset{0}{\overset{1}{\int }}\{U_f(\xi ,T_1)\sin (n\pi \xi )+{\mathrm{\Gamma }}_f(\xi ,T_1)\mathrm{\Gamma }\cos (n\pi \xi )\}d\xi =0$$

(24-a)

$$\underset{0}{\overset{1}{\int }}\{U_b(\xi ,T_1)\sin (n\pi \xi )+{\mathrm{\Gamma }}_b(\xi ,T_1)\mathrm{\Gamma }\cos (n\pi \xi )\}d\xi =0$$

(24-b)

Inserting from equations (23), the solvability conditions are finally determined as

$$\begin{array}{l}i{\mathrm{\Lambda }}_1^f{D}_1{F}_f(T_1)+{\mathrm{\Lambda }}_2^f{F}_f(T_1)\\ +{\mathrm{\Lambda }}_3[{\overline{F}}_f(T_1)F_f^2(T_1)+2{\overline{F}}_b(T_1)F_f(T_1)F_b(T_1)]=0\end{array}$$

(25-a)

$$\begin{array}{l}i{\mathrm{\Lambda }}_1^b{D}_1{F}_b(T_1)+{\mathrm{\Lambda }}_2^b{F}_b(T_1)\\ +{\mathrm{\Lambda }}_3[{\overline{F}}_b(T_1)F_b^2(T_1)+2{\overline{F}}_f(T_1)F_b(T_1)F_f(T_1)]=0\end{array}$$

(25-b)

wherein the coefficients are as follows

$${\mathrm{\Lambda }}_1^f=8(\lambda +n^2{\pi }^2{)}^2(-2{\omega }_f+\beta n^2{\pi }^2)$$

(26-a)

$${\mathrm{\Lambda }}_2^f=16{\omega }_f^2{n}^2{\pi }^2(\lambda n^2{\pi }^2+\frac{1}{2}{\lambda }^2+n^4{\pi }^4)$$

(26-b)

$${\mathrm{\Lambda }}_3=n^4{\pi }^4(\lambda +n^2{\pi }^2{)}^2$$

(26-c)

$${\mathrm{\Lambda }}_1^b=8(\lambda +n^2{\pi }^2{)}^2(-2{\omega }_b+\beta n^2{\pi }^2)$$

(26-d)

$${\mathrm{\Lambda }}_2^b=16{\omega }_b^2{n}^2{\pi }^2(\lambda n^2{\pi }^2+\frac{1}{2}{\lambda }^2+n^4{\pi }^4)$$

(26-e)

The complex valued functions $F_f(T_1)$ and $F_b(T_1)$ are written in polar form to solve the differential equations which govern the solvability conditions, namely

$$F_f(T_1)=f_f(T_1)e^{i{\theta }_f(T_1)}$$

(27-a)

$$F_b(T_1)=f_b(T_1)e^{i{\theta }_b(T_1)}$$

(27-b)

where $f(T_1)$ and $\theta (T_1)$ stand for the amplitudes and phase angles of the time-dependent part of the response, respectively. Inserting equations (27) in the solvability equations (25) and equating the real and imaginary parts to zero, these relations named modulation equations are acquired

$${\mathrm{\Lambda }}_1^f{D}_1{f}_f(T_1)=0$$

(28-a)

$${\mathrm{\Lambda }}_1^b{D}_1{f}_b(T_1)=0$$

(28-b)

$${\mathrm{\Lambda }}_1^f{D}_1{\theta }_f(T_1)-{\mathrm{\Lambda }}_2^f-{\mathrm{\Lambda }}_3{f}_f^2(T_1)-2{\mathrm{\Lambda }}_3{f}_b^2=0$$

(28-c)

$${\mathrm{\Lambda }}_1^b{D}_1{\theta }_b(T_1)-{\mathrm{\Lambda }}_2^b-{\mathrm{\Lambda }}_3{f}_b^2(T_1)-2{\mathrm{\Lambda }}_3{f}_f^2=0$$

(28-d)

From the first two equations, it can be concluded that the amplitude of the functions $F_f(T_1)$ and $F_b(T_1)$ is constant, namely, $f_f(T_1)=C_f$ and $f_b(T_1)=C_b$ . Employing this result in equations (28-c) and (28-d) and integrating the resultants, the phase angle of the response is obtained as

$${\theta }_f(T_1)=\frac{{\mathrm{\Lambda }}_3{C}_f^2+2{\mathrm{\Lambda }}_3{C}_b^2+{\mathrm{\Lambda }}_2^f}{{\mathrm{\Lambda }}_1^f}{T}_1+D_f$$

(29-a)

$${\theta }_b(T_1)=\frac{{\mathrm{\Lambda }}_3{C}_b^2+2{\mathrm{\Lambda }}_3{C}_f^2+{\mathrm{\Lambda }}_2^b}{{\mathrm{\Lambda }}_1^b}{T}_1+D_b$$

(29-b)

in which the constants $D_f$ and $D_b$ can be evaluated from the initial conditions, so the amplitude and phase of the response are determined and the final form of response is obtained as

$$U_1(\xi ,T_0,T_1)=\left[f_f(T_1)e^{i({\theta }_f(T_1)+{\omega }_f{T}_0)}+f_b(T_1)e^{i({\theta }_b(T_1)+{\omega }_b{T}_0)}\right]\sin (n\pi \xi )$$

(30-a)

$${\gamma }_1(\xi ,T_0,T_1)=\mathrm{\Gamma }\left[f_f(T_1)e^{i({\theta }_f(T_1)+{\omega }_f{T}_0)}+f_b(T_1)e^{i({\theta }_b(T_1)+{\omega }_b{T}_0)}\right]\cos (n\pi \xi )$$

(30-b)

By equating the real and imaginary parts in the complex variable definitions (13), the lateral deflections and shear components can be obtained. Because $T_1=\epsilon T_0$ , the nonlinear natural frequencies in forward and backward whirling modes are defined as

$${\omega }_{\mathrm{Nf}}=\frac{{\mathrm{\Lambda }}_3{C}_f^2+2{\mathrm{\Lambda }}_3{C}_b^2+{\mathrm{\Lambda }}_2^f}{{\mathrm{\Lambda }}_1^f}\epsilon +{\omega }_f$$

(31-a)

$${\omega }_{\mathrm{Nb}}=\frac{{\mathrm{\Lambda }}_3{C}_b^2+2{\mathrm{\Lambda }}_3{C}_f^2+{\mathrm{\Lambda }}_2^b}{{\mathrm{\Lambda }}_1^b}\epsilon +{\omega }_b$$

(31-b)

It can be seen that the vibration frequency is a function of its amplitude which is caused by the nonlinearity. Furthermore, as it can be clearly seen, the linear natural frequencies in forward and backward whirling modes are simultaneously affected the nonlinear free vibration behavior of the system. It should be noted that if the dependent part of natural frequency to the vibration amplitude is neglected, the linear natural frequency of the system is obtained as is performed in Mirtalaie and Hajabasi. ³⁴

Numerical results

In the numerical analysis, the required parameters are set as $E=200\mathrm{GPa}$ , $\nu =0.3$ , $\rho =7860\mathrm{kg}/{\mathrm{m}}^2$ , $R=0.015\mathrm{m}$ , and $L=1\mathrm{m}$ . The shear modulus is determined as $G=E/2(1+\nu )$ in which $\nu$ is Poisson’s ratio. Also, the shear correction factor used in the Timoshenko beam theory is obtained for circular cross section as $\kappa =(6(1+\nu ))/(7+6\nu )$ . The forward and backward nonlinear natural frequencies of the rotor obtained from the present model are compared with the results of the nonlinear Euler–Bernoulli beam model examined in Hosseini and Zamanian ³² in which the effect of shear deformations has been neglected. The results are tabulated in Table 1. It can be seen that the effect of shear deformations decreases the natural frequencies of the rotor. Furthermore, the results are compared in Table 1 with the results obtained from the nonlinear dynamic analysis of spinning Rayleigh beams investigated in Frank Pai et al. ⁹ where the results show good agreement. The effect of rotating speed on the nonlinear natural frequencies has been demonstrated in Table 2 in which the results are compared with their linear counterparts presented in the earlier study of the authors. ³⁴ It can be seen that the natural frequencies are affected by the nonlinearity in the decimal digits. This effect is more evident in greater rotating speeds as well as upper natural frequencies. By inspecting the presented results, it can be concluded that the values of the rotor’s natural frequency are corresponding to the rotor speed. Also, it can be seen that the natural frequencies vary very slowly by the variation in the rotational speed of the shaft. The rotating speed affects the nonlinear higher natural frequencies more. Moreover, for a same rotating speed, the nonlinear natural frequencies in forward and backward whirling modes have small differences in their absolute values, but they differ at the higher modes from each other. Besides, it can be clearly seen that the nonlinearity leads to decrement of the absolute values of natural frequencies comparing with the linear model such that the nonlinear natural frequencies are lower than their linear counterparts. This difference is more intensified in the higher modes of vibrations and greater spin speeds.

Table 1.

Comparison of the nonlinear natural frequencies of the rotor with references ( $\Omega =1000\mathrm{rad}/\mathrm{s}$ ).

Model	${\omega }_1$ (rad/s)		${\omega }_2$ (rad/s)		${\omega }_3$ (rad/s)		${\omega }_4$ (rad/s)
	Backward	Forward	Backward	Forward	Backward	Forward	Backward	Forward
Present	−9.8310	9.8701	−39.0964	39.2523	−87.1273	87.4763	−152.8495	153.4657
Hosseini and Zamanian 32	−9.8358	9.8749	−39.1736	39.3301	−87.5154	87.8676	−154.0659	154.6920
Frank Pai et al. 9	−9.8256	9.8646	−39.3032	39.4591	−88.2406	88.5897	−156.3540	156.9704

Table 2.

Effect of rotating speed on the nonlinear natural frequencies of the rotor and comparison with the linear natural frequencies ³⁴ ( $R/L=0.02$ ).

$\Omega$ (rad/s)	Model	${\omega }_1$ (rad/s)		${\omega }_2$ (rad/s)		${\omega }_3$ (rad/s)		${\omega }_4$ (rad/s)
		Backward	Forward	Backward	Forward	Backward	Forward	Backward	Forward
0	Linear	−9.8553	9.8553	−39.2518	39.2518	−87.6914	87.6914	−154.3786	154.3786
	Nonlinear	−9.8505	9.8505	−39.1743	39.1743	−87.3016	87.3016	−153.1573	153.1573
1000	Linear	−9.8358	9.8749	−39.1736	39.3301	−87.5154	87.8676	−154.0659	154.6920
	Nonlinear	−9.8310	9.8701	−39.0964	39.2523	−87.1273	87.4763	−152.8495	153.4657
2000	Linear	−9.8163	9.8946	−39.0956	39.4086	−87.3399	88.0442	−153.7538	155.0060
	Nonlinear	−9.8115	9.8897	−39.0187	39.3305	−86.9533	87.6514	−152.5424	153.7748
3000	Linear	−9.7968	9.9142	−39.0177	39.4873	−87.1647	88.2212	−153.4423	155.3206
	Nonlinear	−9.7920	9.9093	−38.9411	39.4089	−86.7796	87.8268	−152.2358	154.0844
4000	Linear	−9.7774	9.9339	−38.9400	39.5661	−86.9898	88.3985	−153.1315	155.6359
	Nonlinear	−9.7726	9.9290	−38.8637	39.4874	−86.6063	88.0025	−151.9299	154.3947
5000	Linear	−9.7580	9.9537	−38.8624	39.6451	−86.8153	88.5762	−152.8213	155.9518
	Nonlinear	−9.7532	9.9487	−38.7865	39.5660	−86.4334	88.1786	−151.6246	154.7056
6000	Linear	−9.7387	9.9734	−38.7850	39.7242	−86.6412	88.7543	−152.5117	156.2684
	Nonlinear	−9.7339	9.9685	−38.7094	39.6448	−86.2608	88.3551	−151.3199	155.0171
7000	Linear	−9.7193	9.9933	−38.7078	39.8035	−86.4674	88.9327	−152.2028	156.5855
	Nonlinear	−9.7146	9.9883	−38.6324	39.7238	−86.0885	88.5319	−151.0158	155.3292

The results are also presented for the stationary beam, that is, $\Omega =0\mathrm{rad}/\mathrm{s}$ where the values represent the nonlinear natural frequencies of the nonlinear Timoshenko beam in which the nonlinearities are originated from the inertial terms in addition to the stretching of the beam centerline. Furthermore, the effect of radius-to-length ratio on the nonlinear natural frequencies has been demonstrated in Table 3 and the results are compared with the linear natural frequencies. ³⁴ These results demonstrate the effect of slenderness ratio of the shaft on its linear and nonlinear natural frequencies. It is clearly seen that the nonlinear natural frequencies decreased by the increase in the value of the radius-to-length ratio of the shaft. Furthermore, the effect is more observed in higher modes. It can be clearly seen that the predicted linear and nonlinear natural frequencies are nearly similar in the fundamental mode, but they deviate more in higher frequencies. For higher modes, the present results are less than the frequencies obtained from the linear model. This shows that the effect of stretching nonlinearity of the rotor decreases the stiffness of the system in the higher frequency vibrations. Also, this shows that the effect of stretching nonlinearity is more evident at higher frequencies.

Table 3.

Effect of radius-to-length ratio on the nonlinear natural frequencies of the rotor and comparison with the linear natural frequencies ³⁴ ( $\Omega =1500\mathrm{rad}/\mathrm{s}$ ).

$\frac{R}{L}$	Model	${\omega }_1$ (rad/s)		${\omega }_2$ (rad/s)		${\omega }_3$ (rad/s)		${\omega }_4$ (rad/s)
		Backward	Forward	Backward	Forward	Backward	Forward	Backward	Forward
0.01	Linear	−9.8514	9.8807	−39.3627	39.4801	−88.4066	88.6707	−156.7726	157.2422
	Nonlinear	−9.8502	9.8795	−39.3433	39.4606	−88.3086	88.5721	−156.4635	156.9313
0.02	Linear	−9.8260	9.8847	−39.1346	39.3694	−87.4276	87.9559	−153.9098	154.8489
	Nonlinear	−9.8212	9.8798	−39.0576	39.2914	−87.0402	87.5638	−152.6959	153.6202
0.03	Linear	−9.7937	9.8817	−38.7982	39.1504	−85.9366	86.7290	−149.5739	150.9826
	Nonlinear	−9.7829	9.8707	−38.6266	38.9756	−85.0791	85.8556	−146.9055	148.2638
0.04	Linear	−9.7544	9.8718	−38.3606	38.8302	−84.0040	85.0605	−144.1107	145.9890
	Nonlinear	−9.7353	9.8522	−38.0590	38.5211	−82.5077	83.5264	−139.4783	141.2350
0.05	Linear	−9.7084	9.8552	−37.8305	38.4175	−81.7145	83.0352	−137.8927	140.2406
	Nonlinear	−9.6787	9.8245	−37.3653	37.9377	−79.4203	80.6662	−130.8054	132.9098

The variations in the nonlinear natural frequencies in forward and backward whirling modes with respect to the spin speed of the rotor are demonstrated in Figures 4–6 for various radius of the rotor to its length ratio in the first three modes of vibration. It can be concluded that the effect of nonlinearity is intensified at the higher modes of vibration. The nonlinear natural frequencies are increased as the spin speed and radius-to-length ratio of the rotor increased. The variation in the nonlinear natural frequencies with respect to the radius-to-length ratio of the shaft has a linear form. Also, it can be seen that for some values of the rotating speed lower than a specific value, the forward nonlinear natural frequencies are increased as the radius-to-length ratio of the rotor decreased, and for upper values of the rotating speeds than the specific value, the changes in the parameter affect the forward nonlinear natural frequencies inversely. Furthermore, the above-mentioned specific value is increased for the higher vibrational modes. This result is not observed for backward nonlinear natural frequencies. The effect of this parameter on the forward nonlinear natural frequencies is intensified as the higher modes are excited also for higher spin speeds of the rotor (greater than the specific value for forward nonlinear natural frequencies).

Figure 4.

The effect of radius-to-length ratio and spin speed on the first nonlinear natural frequency.

Figure 5.

The effect of radius-to-length ratio and spin speed on the second nonlinear natural frequency.

Figure 6.

The effect of radius-to-length ratio and spin speed on the third nonlinear natural frequency.

Free vibration response

Keeping in mind the boundary conditions (14), by assuming the mode shapes of the system, the solutions of equations of motion (11) is defined as

$$u(\xi ,\tau )=\sum _{n=1}^N{u}_n(\tau )\sin (2n\pi \xi )$$

(32-a)

$$U(\xi ,\tau )=\sum _{n=1}^N{U}_n(\tau )\sin (n\pi \xi )$$

(32-b)

$$\gamma (\xi ,\tau )=\sum _{n=1}^N{\mathrm{\Gamma }}_n(\tau )\cos (n\pi \xi )$$

(32-c)

Inserting these assumed solutions in the governing equations and keeping in mind that the mode shapes are orthogonal to each other lead to

$$\begin{array}{c}(-\frac{1}{2}{n}^3{\pi }^3+\epsilon n^5{\pi }^5)V_n^2+(-\frac{1}{2}{n}^3{\pi }^3+\epsilon n^5{\pi }^5)W_n^2-4n^2{\pi }^2{u}_n-\epsilon {\stackrel{··}{u}}_n\hfill \\ +[{\epsilon }^2{n}^3{\pi }^3{\stackrel{··}{V}}_n-\epsilon n^4{\pi }^4{\mathrm{\Gamma }}_{\mathrm{zn}}+2\Omega {\epsilon }^2{n}^3{\pi }^3{\stackrel{·}{W}}_n+{\epsilon }^2{n}^2{\pi }^2({\stackrel{··}{\mathrm{\Gamma }}}_{\mathrm{yn}}+\Omega {\stackrel{·}{\mathrm{\Gamma }}}_{\mathrm{zn}})\sin (\Omega \tau )+{\epsilon }^2{n}^2{\pi }^2(-{\stackrel{··}{\mathrm{\Gamma }}}_{\mathrm{zn}}+\Omega {\stackrel{·}{\mathrm{\Gamma }}}_{\mathrm{yn}})\cos (\Omega \tau )]V_n\hfill \\ +[{\epsilon }^2{n}^2{\pi }^2(-{\stackrel{··}{\mathrm{\Gamma }}}_{\mathrm{zn}}+\Omega {\stackrel{·}{\mathrm{\Gamma }}}_{\mathrm{yn}})\sin (\Omega \tau )-{\epsilon }^2{n}^2{\pi }^2({\stackrel{··}{\mathrm{\Gamma }}}_{\mathrm{yn}}+\Omega {\stackrel{·}{\mathrm{\Gamma }}}_{\mathrm{zn}})\cos (\Omega \tau )-2{\epsilon }^2{n}^3{\pi }^3\Omega {\stackrel{·}{V}}_n+{\epsilon }^2{n}^3{\pi }^3{\stackrel{··}{W}}_n-\epsilon n^4{\pi }^4{\mathrm{\Gamma }}_{\mathrm{yn}}]W_n=0\hfill \end{array}$$

(33-a)

$$\begin{array}{c}-\beta n^3{\pi }^3{\stackrel{·}{u}}_n{V}_n-\frac{1}{\epsilon }[-2\epsilon n^5{\pi }^5{u}_n-{\epsilon }^2{n}^3{\pi }^3{\stackrel{··}{u}}_n+\epsilon n^4{\pi }^4+n^3{\pi }^3{u}_n]W_n+n^3{\pi }^3{\mathrm{\Gamma }}_{\mathrm{yn}}+n^2{\pi }^2\beta {\stackrel{·}{V}}_n\hfill \\ -[\epsilon n\pi \Omega -\epsilon n^2{\pi }^2\Omega u_n]\sin (\Omega \tau ){\stackrel{·}{\mathrm{\Gamma }}}_{\mathrm{yn}}-[-\epsilon n\pi \Omega +\epsilon n^2{\pi }^2\Omega u_n]\cos (\Omega \tau ){\stackrel{·}{\mathrm{\Gamma }}}_{\mathrm{zn}}\hfill \\ -[-\epsilon n\pi +\epsilon n^2{\pi }^2{u}_n]\cos (\Omega \tau ){\stackrel{··}{\mathrm{\Gamma }}}_{\mathrm{yn}}-[-\epsilon n\pi +\epsilon n^2{\pi }^2{u}_n]\sin (\Omega \tau ){\stackrel{··}{\mathrm{\Gamma }}}_{\mathrm{zn}}\hfill \\ -[\epsilon n^2{\pi }^2-2\epsilon n^3{\pi }^3{u}_n+\epsilon ]{\stackrel{··}{W}}_n+2\epsilon n^3{\pi }^3{\stackrel{·}{u}}_n{\stackrel{·}{W}}_n-2\beta n^3{\pi }^3{\stackrel{·}{V}}_n-n^4{\pi }^4{\mathrm{\Gamma }}_{\mathrm{yn}}{u}_n=0\hfill \end{array}$$

(33-b)

$$\begin{array}{c}\beta n^3{\pi }^3{\stackrel{·}{u}}_n{W}_n+\frac{1}{\epsilon }[{\epsilon }^2{n}^3{\pi }^3{\stackrel{··}{u}}_n-n^3{\pi }^3{u}_n+2\epsilon n^5{\pi }^5{u}_n-\epsilon n^4{\pi }^4]V_n+2\beta n^3{\pi }^3{\stackrel{·}{W}}_n{u}_n\hfill \\ +[-\epsilon n\pi \Omega +\epsilon n^2{\pi }^2\Omega u_n]\cos (\Omega \tau ){\stackrel{·}{\mathrm{\Gamma }}}_{\mathrm{yn}}+[-\epsilon n\pi \Omega +\epsilon n^2{\pi }^2\Omega u_n]\sin (\Omega \tau ){\stackrel{·}{\mathrm{\Gamma }}}_{\mathrm{zn}}\hfill \\ +[-\epsilon n\pi +\epsilon n^2{\pi }^2{u}_n]\sin (\Omega \tau ){\stackrel{··}{\mathrm{\Gamma }}}_{\mathrm{yn}}+[\epsilon n\pi -\epsilon n^2{\pi }^2{u}_n]\cos (\Omega \tau ){\stackrel{··}{\mathrm{\Gamma }}}_{\mathrm{zn}}-\beta n^2{\pi }^2{\stackrel{·}{W}}_n\hfill \\ +[-\epsilon n^2{\pi }^2+2\epsilon n^3{\pi }^3{u}_n-1]{\stackrel{··}{V}}_n+[n^3{\pi }^3-n^4{\pi }^4{u}_n]{\mathrm{\Gamma }}_{\mathrm{zn}}+\beta n^3{\pi }^3{\stackrel{·}{u}}_n{W}_n+2\epsilon n^3{\pi }^3{\stackrel{·}{u}}_n{\stackrel{·}{V}}_n=0\hfill \end{array}$$

(33-c)

$$\begin{array}{c}[\epsilon n^2{\pi }^2{\stackrel{··}{u}}_n\sin (\Omega \tau )+\epsilon n^2{\pi }^2\Omega {\stackrel{·}{u}}_n\cos (\Omega \tau )]V_n-\epsilon {\stackrel{··}{\mathrm{\Gamma }}}_{\mathrm{yn}}-(\lambda +n^2{\pi }^2){\mathrm{\Gamma }}_{\mathrm{yn}}-\epsilon n\pi \Omega \sin (\Omega \tau ){\stackrel{·}{W}}_n\hfill \\ +[n^3{\pi }^3-n^4{\pi }^4{u}_n-\epsilon n^2{\pi }^2{\stackrel{··}{u}}_n\cos (\Omega \tau )+\epsilon n^2{\pi }^2\Omega {\stackrel{·}{u}}_n\sin (\Omega \tau )]W_n+2\epsilon n^2{\pi }^2\sin (\Omega \tau ){\stackrel{·}{u}}_n{\stackrel{·}{V}}_n\hfill \\ +[-\epsilon n\pi +\epsilon n^2{\pi }^2{u}_n]\sin (\Omega \tau ){\stackrel{··}{V}}_n+[\epsilon n\pi -\epsilon n^2{\pi }^2{u}_n]\cos (\Omega \tau ){\stackrel{··}{W}}_n-\epsilon n\pi \Omega \cos (\Omega \tau ){\stackrel{·}{V}}_n\hfill \\ +\epsilon n^2{\pi }^2\Omega \sin (\Omega \tau )u_n{\stackrel{·}{W}}_n+\epsilon n^2{\pi }^2\Omega \cos (\Omega \tau )u_n{\stackrel{·}{V}}_n-2\epsilon n^2{\pi }^2\cos (\Omega \tau ){\stackrel{·}{u}}_n{\stackrel{·}{W}}_n=0\hfill \end{array}$$

(33d)

$$\begin{array}{c}[n^3{\pi }^3-n^4{\pi }^4{u}_n-\epsilon n^2{\pi }^2{\stackrel{··}{u}}_n\cos (\Omega \tau )+\epsilon n^2{\pi }^2\Omega {\stackrel{·}{u}}_n\sin (\Omega \tau )]V_n-\epsilon {\stackrel{··}{\mathrm{\Gamma }}}_{\mathrm{zn}}-\epsilon n\pi \Omega \sin (\Omega \tau ){\stackrel{·}{V}}_n\hfill \\ -[\epsilon n^2{\pi }^2{\stackrel{··}{u}}_n\sin (\Omega \tau )+\epsilon n^2{\pi }^2\Omega {\stackrel{·}{u}}_n\cos (\Omega \tau )]W_n+[\epsilon n\pi -\epsilon n^2{\pi }^2{u}_n]\cos (\Omega \tau ){\stackrel{··}{V}}_n\hfill \\ +[\epsilon n\pi -\epsilon n^2{\pi }^2{u}_n]\sin (\Omega \tau ){\stackrel{··}{W}}_n+\epsilon n^2{\pi }^2\Omega \sin (\Omega \tau )u_n{\stackrel{·}{V}}_n-\epsilon n^2{\pi }^2\Omega \cos (\Omega \tau )u_n{\stackrel{·}{W}}_n\hfill \\ -2\epsilon n^2{\pi }^2\cos (\Omega \tau ){\stackrel{·}{u}}_n{\stackrel{·}{V}}_n-2\epsilon n^2{\pi }^2\sin (\Omega \tau ){\stackrel{·}{u}}_n{\stackrel{·}{W}}_n+\epsilon n\pi \Omega \cos (\Omega \tau ){\stackrel{·}{W}}_n-(\lambda +n^2{\pi }^2){\mathrm{\Gamma }}_{\mathrm{zn}}=0\hfill \end{array}$$

(33-e)

The series solutions are truncated and $5N$ coupled linear variable-coefficient ordinary differential equations are obtained that are solved numerically using Runge–Kutta method, and the free vibration response of the system to initial conditions is determined. The numerical results and simulation of the system are represented in this section. The time response of the system to initial conditions using three modes is shown in Figures 7–9, where it is compared with the results of the multiple scales method. It can be obviously seen that the results coincide with each other. A beating phenomenon is observed in the time histories of the system which is originated from the fact that the differences between the absolute values of the nonlinear natural frequencies of the forward and backward whirling modes are negligible. By comparing the results in various modes, it is clearly seen that the frequency of the beating is increased as the higher modes are excited.

Figure 7.

Time response of the first mode at $\xi =0.3$ .

Figure 8.

Time response of the second mode at $\xi =0.3$ .

Figure 9.

Time response of the third mode at $\xi =0.3$ .

Conclusion

The nonlinear variable-coefficient governing differential equations of motion and corresponding boundary conditions were derived using Hamilton’s principle based on the new expression of cross-sectional rotation angles in Timoshenko beam theory. The nonlinearities were originated from the stretching of the beam centerline, and the model includes the superposition of the shear deformations on the elastic deformation of the beam centerline which is represented by the Euler angles. The nonlinear free axial–lateral vibrations of the flexible shear deformable rotor with nonlinearities in stretching were investigated using the multiple scales method, and exact values of the nonlinear natural frequencies in the forward and backward modes were determined. The accuracy of the solution method was inspected by comparing the results with the results of two other beam theories, and also, the comparison of the response of the free vibrations of the system with the numerical integration of the governing equations was done. The effect of various parameters such as the rotational speed and radius-to-length ratio of the shaft on the free vibrations of the system was inspected. The conclusions are expressed as follows:

• The mutual interaction of the shear and Euler angles leads to the appearance of variable-coefficient terms in the nonlinear differential equations of the system and generally affects the global dynamic behavior of the Timoshenko rotor model.

• Due to nonlinear nature of the system, the nonlinear natural frequencies are functions of vibration amplitude.

• The nondependent part of the nonlinear natural frequency to the vibration amplitude represents the natural frequency of the linear model.

• The backward and forward whirling modes of linear vibrations are excited simultaneously and affect the nonlinear free vibration behavior of the system.

• The nonlinearity affects the natural frequencies in the decimal digits. This effect is more evident in greater rotating speeds as well as upper natural frequencies.

• For same rotating speeds, the absolute values of the nonlinear natural frequencies in the forward and backward modes of whirling have negligible differences, and therefore, the beat phenomenon is appeared in the time histories of the free vibrations. They differ at the higher modes from each other.

• The absolute values of nonlinear natural frequencies in forward modes are a little greater than their backward counterparts.

• The absolute values of nonlinear natural frequencies are increased in forward modes by the spin speed of the rotor, while they decreased for the backward nonlinear natural frequencies.

• The nonlinearity leads to decrement in the absolute values of natural frequencies comparing with the linear model such that the nonlinear natural frequencies are lower than their linear counterparts. This difference is more intensified in the higher modes of vibrations and greater spin speeds.

• The predicted linear and nonlinear natural frequencies are nearly similar in the fundamental mode but they deviate more in higher frequencies. For higher modes, the present results are less than the frequencies obtained from the linear model. This shows that the effect of stretching nonlinearity of the rotor decreases the stiffness of the system in the higher frequency vibrations. Also, this shows that the effect of stretching nonlinearity is more evident at higher frequencies.

• The nonlinear natural frequencies are increased as the radius-to-length ratio of the rotor increased. The variation in the nonlinear natural frequencies with respect to the radius-to-length ratio of the shaft has a linear form.

• For some values of the rotating speed lower than a specific value, the forward nonlinear natural frequencies are increased as the radius-to-length ratio of the rotor decreased, and for upper values of the rotating speeds than the specific value, the changes in the parameter affect the forward nonlinear natural frequencies inversely. This specific value is increased for the higher vibrational modes. This result is not observed for backward nonlinear natural frequencies. The effect of this parameter on the forward nonlinear natural frequencies is intensified as the higher modes are excited also for higher spin speeds of the rotor (greater than the specific value for forward nonlinear natural frequencies).