1 Introduction

In nonlinear dynamic systems, the response can change qualitatively with small changes in system parameters, so being able to systematically locate these so-called bifurcations is vital in nonlinear mechanical systems. This complication is accentuated in rotating machinery because of gyroscopic coupling which makes the response modes complex and whirl frequencies dependent on the rotor speed. These systems exhibit rich phenomena including period doubling, aperiodicity, and even chaos in the system’s response [1, 5, 37].

Nonlinearity in rotordynamics can originate from sources such as rotor–stator interactions, faults such as cracks, nonlinear bearing forces, and geometric nonlinearity when shafts are flexible [13, 16, 24]. The rotor–stator contacting phenomenon in turbomachinery can result from a loss of blade, supports being misaligned, fatigue cracks shafts [2, 17], and depending on the bearing type can lead to destructive oil whips [4]. By contrast with turbomachinery, long drilling shafts with low stiffness have been documented to demonstrate nonlinearity due to rotor–stator contact [3, 31, 36]. Ehrich [10] described the following phenomena occurring in the nonlinear response of rotating systems: asynchronous response, changing amplitude for constant rotor speed, discontinuous amplitude response as rotor speed changes. Synchronous partial contact response have been observed [5, 8].

The response of an isotropic rotor system is typically in the form of a constant radius circular motion known as whirling—in linear and weakly nonlinear cases this is generally synchronous, meaning it occurs at the speed of the shaft rotation, excited by unbalanced mass. While the presence of friction can lead to a dangerous case known as dry backward whirl [13], more intriguing are responses where the amplitude of shaft displacement from its rotation centre continually oscillates, also known as bouncing motions or intermittent contact cycles [34]. These can be synchronous, in that the period of the response is in a real ratio with the period of shaft rotation. However, the more challenging case is that of asynchronous responses, where the response frequencies bear no obvious relation to the shaft speed [8, 39]. It has been shown that the intermittent contact can be asynchronous, yet it can be estimated by an internal resonance condition, relating integer multiples of whirl frequencies [34]. It is crucial to note that because the underlying linear frequencies depend on the rotor speed, different underlying modes may take part in internal resonance combinations at different drive speeds.

Many approaches have been used to analyse such response. Some research considered fully rigid models for the stator contact. Cole and Keogh [8] analysed the asynchronous intermittent contact phenomenon. The fundamental frequency of the harmonic balance method (HBM) they used was set to the contacting frequency which was found from the contact energy conservation. They concluded that in the high frequency range the number of possible contact modes decreased with increasing damping. The study gave upper and lower limits of unbalance to achieve the periodic contacting motions. Mora et al. [23] employed a rigid impact model for the contacting interactions, and explained the transition into bouncing motions as a grazing bifurcation.

Further studies use smooth or piecewise linear functions to represent nonlinearity due to contact. Karpenko et al. [18] analytically solved the two-degree-of-freedom (2-dof) Jeffcott system. Depending on the accuracy of Taylor expansion during contact, time simulations were approximated very closely. Gyroscopic effects, friction, and snubber mass were not considered. Karpenko et al. [19] observed that higher damping and lower stiffness eliminated multi-periodic and chaotic solutions on a similar model. A preloaded snubber ring in [27] shifted the bifurcation plots, added new bifurcations, and decreased the response amplitudes of the retainer ring. The position of the massless snubber ring was found from the energy spent on the ring during contact, with the help of a Lagrange multiplier. The experimental validation of this mathematical model was published by Karpenko et al. [20].

Zilli et al. [39] investigated a snubbed overhung rotor with 2-dof, and time-simulated a blade loss scenario by introducing an unbalance. The analytical model had a piecewise linear stiffness to define contact. A synchronisation condition was derived using phasors of system frequencies. The frequencies where the intermittent contacts occurred in the bifurcation diagrams were closely approximated with this synchronisation condition. In a rapid deceleration test, higher transient energy was found more likely to yield sustained bouncing motion.

Shaw et al. [34] expanded the internal resonance idea to higher-degree-of-freedom systems by setting the condition to involve any two whirl modes. In the study, soft contact with no damping and stiff contact with damping were investigated on a finite element (FE) model that is further reduced with some of the linear eigenmodes. As observed in [39], the authors reported the stiffening effect that causes the predicted rotor speeds to be lower than the numerical ones. Chipato et al. introduced friction [6] and gravity [5] to the overhung rotor system in [39]; they noted that the synchronisation condition holds in both cases. With gravity, multi-periodic and even chaotic responses appeared. However, increased rotor stiffness nullified the effect of gravity. Varney and Green [37] presented a rotor speed threshold after which the gravity can be ignored.

Analytical and semi-analytical methods to solve the nonlinear equation of motions include the use of normal forms, harmonic balance method (HBM), and alternating frequency/time domain methods (AFT), [35]. In the study by Von Groll and Ewins [38], after the equations of motion were converted to algebraic equations by HBM, arclength continuation was performed to efficiently find solutions, with the help of AFT to update contact forces. The method was demonstrated on a Jeffcott rotor in penalty–contact interaction with a stator of nonzero mass. Subharmonics were allowed in HBM formulations which was stated to be required by the intermittent contact motion.

Quasi-periodic solutions were investigated by Peletan et al. [29] who used HBM with pseudo-arclength continuation, both of which were modified to automate the selection of harmonics and continuation of quasi-periodic solution branches. A previous study [28] showed that for the stability analysis, the time domain methods were more accurate than the frequency domain methods that are based on the modified Hill’s method. A study [30] assessed the limitations of HBM in the rotor–stator contact modelling.

Salles et al. [32] used HBM and numerical continuation to obtain the bifurcation diagram of an entire aeroengine with three shaft and four snubber elements. Component mode synthesis was used to reduce the FE model, which was further reduced to the nonlinear degrees of freedom. A special treatment of the method allowed them to find isolated solutions in bifurcation diagram. However, the analysis of the asynchronous response was not their objective.

The fact that the intermittent contact response can be seen as the synchronisation of resonant frequencies of the rotating system, make the use of the normal forms technique a useful procedure to solve the nonlinear equations [35]. Shaw et al. [35] employed the normal forms method on a reduced order model together with an AFT method to update the nonlinear forces. The proposed method was demonstrated on a 2-dof overhung snubbed rotor system and a multi-disc system involving rotor–stator contact. The importance of the use of rotating coordinates was emphasised to get a periodic response for the bouncing cycles, as also mentioned by Cole and Keogh [8].

Geometric nonlinearity in rotors has been discussed in the context of rotordynamics by Nagasaka et al. [25]. The mathematical model of a continuous slender shaft with nonlinear cross section was constrained longitudinally through bearings at each end. The gyroscopic effects were small due to the lack of a disc. Thus, the 1:1 internal resonance in the stationary frame was investigated for both major and secondary critical speeds. The cause of the nonlinear response such as Hopf bifurcations and quasi-periodic motions was linked to the internal resonance. Findings were experimentally confirmed. Previous works accompanied that study such as [14, 15, 26]. In the book by Ishida and Yamamoto [16], further discussions of geometric nonlinearity and their solutions were given. The book [16], represented the weak nonlinearity in stiffness as one for which the potential energy can be expressed as a power series, whose derivatives with respect to the dimensions of the system give restoring forces. The harmonic and subharmonic responses, as well as combination resonances, were investigated in the stationary coordinate frame.

Aside from experimental investigations, the primary tools for understanding rotordynamics have thus far been time simulation [5,6,7, 12, 34, 37, 39] and analytical and semi-analytical approaches based on the harmonic balance method [29, 30, 32, 35, 38]. The time simulation approach is slow; solutions must be continued for some arbitrary amount of time to converge to a steady state which make it excessively time-consuming for low damping cases [33]. Also, many simulations must be run with different initial conditions if multiple solutions are to be found [34]. Even when this has been done, there is no guarantee that all responses have been found, making time simulation quite an inefficient method to investigate these systems. Analytical and semi-analytical approaches require sometimes complex formulations to transform the equations of motion to algebraic form, and still leave doubt that all solutions have been found, depending on the numerical solver’s start point. However, the use of arclength continuation does offer significant benefits in exploring the space around a given solution in the sense that it allows a complete solution family to be found without the interruption of jumping phenomenon, a common phenomenon for time simulations. It can pass the turning points, the so-called folds, possibly progressing to the unstable region [30], which make the method feasible for high-resolution bifurcation diagrams [33]. However, an initial time to set up the continuation parameters such as those that define accuracy is required in continuation [9]. Given the rich complexity of solutions in rotor systems, it is perhaps surprising the method of numerical continuation, where continuation is applied directly to the equations of motion, has not previously been applied more extensively.

Considering the literature above, the present work seeks the analysis of a 2-dof mathematical model of an overhung rotor with an isotropic cubic stiffness nonlinearity. The main focus is asynchronous periodic orbits that are caused by internal resonance in the case of geometric nonlinearity, and we see similarities to the behaviours seen with contacting nonlinearity. The results are compared to studies with discontinuous stiffness contact definitions. The continuation method is applied directly to the systems of ODE, in contrast to Refs. [29, 32, 38]. The autonomous rotating coordinate frame equations of motion are used in the numerical continuation procedures; the observed periodic solution families correspond to quasi-periodic solutions in the stationary frame.

In the rest of this paper, Sect. 2 presents the theoretical background: the equations of motion are derived using Lagrange’s method, nondimensionalised and transformed to the rotating frame. Section 2 ends with a plot of the natural frequency map featuring signed natural frequencies [35]. Section 3 contains time simulations set-up using the MATLAB function ode45. The numerical continuation procedures are explained in Sect. 4 with a summary of the continuation theory. Then the results of the continuation runs and time simulations are presented in Sect. 5 with an extra validation of the continuation results using time simulations. Conclusions from the results, along with the future work ideas, are given in Sect. 6.

2 Theoretical background

2.1 Derivation of equations of motion

The system subject to analysis here is a two-degree-of-freedom (2-dof) overhung rotor system, as shown in Fig. 1. Here, a and b are distances from point O, up to the disc and the stiffness element, respectively; \(k_r\) is the isotropic linear stiffness; \(k_3\) is the isotropic cubic stiffness; \(c_r\) is the linear damping coefficient; \(\theta ,\psi ,\phi \) are the angles of rotation about x, y, z axes, respectively; \(\varOmega \) is the rotor speed \((\phi =\varOmega t)\); \(I_p\) and \(I_d\) are polar and diametral moments of inertia, respectively. This system is in some ways similar to the contacting system in [39], although the piecewise linear contact is replaced by a smooth cubic stiffness \(k_3\). As an isotropic stiffening effect, this cubic stiffness can represent effects such as a bend–stretch coupling due to fixed positions of the bearings, for example, in a pinned–pinned shaft [16]. Note that the derivation of the equations of motion in the next section assumes constant rotational speed, \({\dot{\varOmega }}=0\), and small transverse angles \(\theta \) and \(\psi \).

Fig. 1
figure 1

A schematic representation of the overhung rotor with cubic stiffness nonlinearity

Full size image

2.2 Equations of motion

In this section, the derivation of equations of motion is given which takes the transverse angles \(\theta \) and \(\psi \) as the generalised coordinates. The derivation uses Lagrange’s energy method [13].

2.2.1 Kinetic energy

Kinetic energy resulting from the translational motion of the disc is,

$$\begin{aligned} T_{\mathrm{tra}} = \frac{1}{2}m_d({\dot{x}}_g^2+{\dot{y}}_g^2) \end{aligned}$$
(1)

where \(m_d\) is the disc mass, \(x_g=x_c+e\cos {\phi }\) and \(y_g=y_c+e\sin {\phi }\) are the transverse translations of the disc centre of gravity,

$$\begin{aligned} {\dot{x}}_g=\dot{x_c}-e\varOmega \sin {\phi }, \quad {\dot{y}}_g=\dot{y_c}+e\varOmega \cos {\phi } \end{aligned}$$
(2)

and \(x_c\) and \(y_c\) point to the disc centre of rotation. For small motions, the angles \(\theta ,\psi \) can be defined in terms of \(x_c,y_c\) and linearised to give,

$$\begin{aligned} x_c=a\psi , \quad y_c=-a\theta . \end{aligned}$$
(3)

It is convenient to express the rotational kinetic energy with respect to the body-fixed frame of reference, i.e. a coordinate system that is fixed onto the disc. The method requires that the Euler angles be obtained from position with zero angles of rotation by applying rotations. Let the order of rotations be:

  • \(\psi \) about y-axis,

  • \(\theta \) about resulting x-axis,

  • \(\phi \) about resulting z-axis.

The transformation matrices needed for these transformations are:

$$\begin{aligned} R_{\theta }= \begin{bmatrix} 1 &{} 0 &{} 0 \\ 0 &{} \cos {\theta } &{} \sin {\theta } \\ 0 &{} \sin {\theta } &{} -\cos {\theta } \end{bmatrix}, \quad R_{\phi }= \begin{bmatrix} \cos {\phi } &{} \sin {\phi } &{} 0 \\ -\sin {\phi } &{} \cos {\phi } &{} 0 \\ 0 &{} 0 &{} 1 \end{bmatrix}.\nonumber \\ \end{aligned}$$
(4)

These matrices are be used to find the angular speeds in the body-fixed frame of reference as,

$$\begin{aligned} \begin{aligned} \begin{pmatrix} \omega _{{\bar{x}}} \\ \omega _{{\bar{y}}} \\ \omega _{{\bar{z}}} \end{pmatrix} =&\begin{pmatrix} 0 \\ 0 \\ \varOmega \end{pmatrix} + R_\phi \begin{pmatrix} {\dot{\theta }} \\ 0 \\ 0 \end{pmatrix} + R_\phi R_\theta \begin{pmatrix} 0 \\ {\dot{\psi }} \\ 0 \end{pmatrix} \\=&\begin{pmatrix} {\dot{\theta }}\cos {\phi }+{\dot{\psi }}\sin {\phi }\cos {\theta } \\ -{\dot{\theta }}\sin {\phi }+{\dot{\psi }}\cos {\phi }\cos {\theta } \\ {\varOmega - {\dot{\psi }} \sin {\theta }} \end{pmatrix} \end{aligned} \end{aligned}$$
(5)

where \(\omega _{{\bar{x}}},\omega _{{\bar{y}}},\omega _{{\bar{z}}}\) are the angular speed with respect to the body-fixed frame. As a result, the rotational kinetic energy can be expressed as,

$$\begin{aligned} \begin{aligned} T_{\mathrm{rot}}&= \frac{1}{2} I_d \left( \omega _{{\bar{x}}}^2 + \omega _{{\bar{y}}}^2 \right) {+\frac{1}{2} I_p \omega _{{\bar{z}}}^2 } \\&= \frac{1}{2} I_d ( {\dot{\theta }}^2 + {\dot{\psi }}^2 \cos ^2{\theta } ) \\&\quad + \frac{1}{2} I_p ( \varOmega ^2 -2\varOmega {\dot{\psi }}\sin {\theta } +{\dot{\psi }}^2\sin ^2{\theta } ). \end{aligned} \end{aligned}$$
(6)

Therefore, the total kinetic energy of the disc becomes \(T_{d}=T_{\mathrm{tra}}+T_{\mathrm{rot}}\), with small angles \(\theta \) and \(\psi \):

$$\begin{aligned} \begin{aligned} T_d&= \frac{1}{2} I_d^o \left( {{\dot{\theta }}}^2 + {{\dot{\psi }}}^2 \right) + \frac{1}{2} m_d e^2 \varOmega ^2 \\&\quad - m_d a e \varOmega \left( {\dot{\psi }}\sin {\phi }+{\dot{\theta }}\cos {\phi } \right) \\&\quad + \frac{1}{2} I_p \varOmega ^2 - I_p\varOmega {\dot{\psi }}\theta \end{aligned} \end{aligned}$$
(7)

where \(I_d^o=I_d+m_da^2\) is the diametral moment of inertia of the disc with respect to centre of the fixed-frame of reference.

2.2.2 Potential energy

All potential energy is due to the nonlinear spring, acting in the radial direction. The conservative force can be expressed as the superposition of the linear and cubic stiffness contributions,

$$\begin{aligned} f_r = - k_r r - k_3 r^3 \end{aligned}$$
(8)

where \(r=b \sqrt{\theta ^2+\psi ^2}\). Therefore, the potential energy is given by:

$$\begin{aligned} V = -\int _{0}^{r}{f_r \mathrm {d}r} = \frac{1}{2}k_r r^2+\frac{1}{4}k_3 r^4. \end{aligned}$$
(9)

2.2.3 Lagrangian equations of motion

The Lagrangian function \(L=T-V\) is found from Eqs. (7) and (9) as,

$$\begin{aligned} \begin{aligned} L&= \frac{1}{2}I_d^o\left( {{\dot{\theta }}}^2+{{\dot{\psi }}}^2\right) + \frac{1}{2}m_de^2\varOmega ^2 \\&\quad -m_dae\varOmega \left( {\dot{\psi }}\sin {\phi } + {\dot{\theta }}\cos {\phi }\right) \\&\quad +\frac{1}{2}I_p\varOmega ^2 - I_p\varOmega {\dot{\psi }}\theta {- \frac{1}{2}k_r r^2 - \frac{1}{4}k_3 r^4}. \end{aligned}\end{aligned}$$
(10)

The following nonconservative forces are introduced to include the effects of linear damping in the restoring force:

$$\begin{aligned} \begin{aligned} Q_1&=-c_rb^2{\dot{\theta }} \\ Q_2&=-c_rb^2{\dot{\psi }}. \end{aligned} \end{aligned}$$
(11)

The 2nd type Lagrange equations,

$$\begin{aligned} \frac{\mathrm {d}}{\mathrm {d}t}\left( \frac{\partial L}{\partial {{\dot{q}}}_k}\right) -\frac{\partial L}{\partial q_k}=Q_k, \quad k=1,2 \end{aligned}$$
(12)

yield the following equations of motion in the angular coordinates, \((q_1,q_2)=(\theta ,\psi )\),

$$\begin{aligned} \begin{aligned}&I_d^o\ddot{\theta }+I_p\varOmega {\dot{\psi }}+c_rb^2{\dot{\theta }}\\&\quad +k_rb^2\theta +k_3b^4\left( \theta ^2+\psi ^2\right) \theta =-m_dea\varOmega ^2\sin {\phi } \\&I_d^o\ddot{\psi }-I_p\varOmega {\dot{\theta }}+c_rb^2{\dot{\psi }}\\&\quad +k_rb^2\psi +k_3b^4\left( \theta ^2+\psi ^2\right) \psi =m_dea\varOmega ^2\cos {\phi }. \end{aligned} \end{aligned}$$
(13)

2.3 Nondimensionalised equations of motion

A nondimensionalisation of Eq. (13) is performed by dividing by \(k_t=k_rb^2\), which represents the angular equivalent of the radial stiffness \(k_r\), and it is related to the natural frequency of the underlying linear system through \(\omega _n^2=k_t/I_d^o\). A new term \(\gamma \) is introduced to express the ratio of linear and cubic returning forces:

$$\begin{aligned} \gamma =\frac{k_3b^2}{k_r}. \end{aligned}$$
(14)

If one keeps b constant, the ratio \(\gamma \) is a direct measure of cubic stiffness, which is the geometric nonlinearity used here. Thus, the nondimensional equations of motion are,

$$\begin{aligned} \begin{aligned}&{\hat{\theta }}^{\prime \prime } + \widehat{I_p}\hat{\varOmega }{\hat{\psi }}^{\prime \prime } + 2\zeta {\hat{\theta }}^\prime \\&\quad + \hat{\theta } + \gamma \left( {\hat{\theta }}^2 + {\hat{\psi }}^2\right) \hat{\theta } = -{\hat{m}}{\hat{e}}{\hat{\varOmega }}^2\sin {\phi } \\&{\hat{\psi }}^{\prime \prime } - \widehat{I_p}\hat{\varOmega }{\hat{\theta }}^{\prime \prime } + 2\zeta {\hat{\psi }}^\prime \\&\quad + \hat{\psi } + \gamma \left( {\hat{\theta }}^2 + {\hat{\psi }}^2\right) \hat{\psi } = {\hat{m}}{\hat{e}}{\hat{\varOmega }}^2\cos {\phi } \end{aligned} \end{aligned}$$
(15)

where the following identities are used,

$$\begin{aligned} \begin{aligned}&\frac{I_p}{I_d^o}=\widehat{I_p},\ \ {{\hat{m}}=\frac{m_da^2}{I_d^o}},\ \ {\hat{e}}=\frac{e}{a},\ \ \zeta =\frac{c_rb^2}{2I_d^o\omega _n}\\&\tau =t{\omega _n},\ \ (\cdot )^\prime =\frac{(\cdot )^\cdot }{\omega _n},\ \ \hat{\varOmega }=\frac{\varOmega }{\omega _n},\ \ \hat{\theta }=\theta ,\ \ \hat{\psi }=\psi . \end{aligned} \end{aligned}$$
(16)

Hats on \(\theta \) and \(\psi \) were placed for consistency within the equations. Here, \((\cdot )^\prime \) is the derivative with respect to the nondimensionalised time, \(\tau \).

The parameters were chosen to provide a system that can be reasonably compared to the system presented in [39]. Clearly, there is no direct equivalence between the coefficients describing the smooth nonlinearity in this work with the piecewise nonlinear stiffness in [39], but the values chosen give oscillating limit cycles at similar amplitudes found by Zilli. A comparison of the nonlinear functions used in this work and [39] is shown in Fig. 2.

Fig. 2
figure 2

a Restoring force and b potential energy comparison between piecewise contact stiffness model that is used in [39] and cubic stiffness model used here

Full size image

2.4 Stationary frame equations in displacement coordinates

The displacement of the disc in the stationary frame can be derived from the equation of motion in Eq. (15) by simply applying the transformation,

$$\begin{aligned} \hat{\theta }=-{\hat{y}},\quad \hat{\psi }={\hat{x}} \end{aligned}$$
(17)

where \({\hat{x}}=\frac{x_c}{a},\ {\hat{y}}=\frac{y_c}{a}\) from Eqs. (3) and (16).

The nondimensionalised stationary frame displacement equations of motion are then obtained as:

$$\begin{aligned}&{{\hat{x}}}^{\prime \prime }+\widehat{I_p}\hat{\varOmega }{{\hat{y}}}^\prime +2\zeta {{\hat{x}}}^\prime +{\hat{x}}+\gamma \left( {{\hat{x}}}^2+{{\hat{y}}}^2\right) {\hat{x}}={\hat{m}}{\hat{e}}{\hat{\varOmega }}^2\cos {\phi }\nonumber \\&{{\hat{y}}}^{\prime \prime }-\widehat{I_p}\hat{\varOmega }{{\hat{x}}}^\prime +2\zeta {{\hat{y}}}^\prime +{\hat{y}}+\gamma \left( {{\hat{x}}}^2+{{\hat{y}}}^2\right) {\hat{y}}={\hat{m}}{\hat{e}}{\hat{\varOmega }}^2\sin {\phi }. \end{aligned}$$
(18)

In matrix form with \(\hat{{\mathbf {X}}}=[{\hat{x}},{\hat{y}}]^T\) this is,

$$\begin{aligned} \hat{{\mathbf {X}}}^{\prime \prime } - \hat{\varOmega }\widehat{I_p}{\mathbf {J}}{\hat{{\mathbf {X}}}}^\prime + 2\zeta {\hat{{\mathbf {X}}}}^\prime + \hat{{\mathbf {X}}} + \gamma {{\hat{r}}}^2\hat{{\mathbf {X}}} = {\hat{m}}{\hat{e}}{\hat{\varOmega }}^2\begin{Bmatrix}\ \cos {\phi }\\ \sin {\phi }\\ \end{Bmatrix} \end{aligned}$$
(19)

where \({\hat{r}}=\sqrt{{{\hat{x}}}^2+{{\hat{y}}}^2}=r/b\) is the distance of the centre of the disc from the equilibrium position in the nondimensional form, and \({\mathbf {J}}\) is the skew symmetric matrix:

$$\begin{aligned} {\mathbf {J}}=\left[ \begin{matrix}0&{}-1\\ 1&{}0\\ \end{matrix}\right] . \end{aligned}$$
(20)

2.5 Rotating frame equations of motion

A transformation in \(\phi \), the rotor angular position, has to be applied in the following form,

$$\begin{aligned} \hat{{\mathbf {X}}}={\mathbf {T}}\hat{{\mathbf {U}}} \end{aligned}$$
(21)

where \(\hat{{\mathbf {U}}}=\left[ {\hat{u}},{\hat{v}}\right] ^T\) is the rotating coordinate vector, and the transformation \({\mathbf {T}}\) is,

$$\begin{aligned} {\mathbf {T}}=\begin{bmatrix}\cos {\phi }&{}-\sin {\phi }\\ \sin {\phi }&{}\cos {\phi }\\ \end{bmatrix}, \quad {\mathbf {T}}^{-1}={\mathbf {T}}^T. \end{aligned}$$
(22)

During the transformations, the first and second normalised derivatives of \({\mathbf {X}}\) are needed,

$$\begin{aligned} \begin{aligned}&\hat{{\mathbf {X}}}={\mathbf {T}}\hat{{\mathbf {U}}},\quad \hat{{\mathbf {X}}}^{\prime } = {\mathbf {T}}^{\prime }\hat{{\mathbf {U}}}+T\hat{{\mathbf {U}}}^{\prime },\\&{\mathbf {X}}^{\prime \prime }={\mathbf {T}}^{\prime \prime }\hat{{\mathbf {U}}} + 2{\mathbf {T}}^{\prime }\hat{{\mathbf {U}}}^\prime + {\mathbf {T}}\hat{{\mathbf {U}}}^{\prime \prime } \end{aligned} \end{aligned}$$
(23)

where the first and second normalised derivative of \({\mathbf {T}}\) are,

$$\begin{aligned} {\mathbf {T}}^\prime =\hat{\varOmega }{{\mathbf {T}}}{{\mathbf {J}}}, \quad {\mathbf {T}}^{\prime \prime }=-{\hat{\varOmega }}^2{\mathbf {T}}. \end{aligned}$$
(24)

After substituting Eqs. (23) and (24) to Eq. (19), and multiplying on the left by \({\mathbf {T}}^T\), the rotating frame equations of motion are obtained as,

$$\begin{aligned} \begin{aligned}&\hat{{\mathbf {U}}}^{\prime \prime } + \big (-\hat{\varOmega }{\mathbf {J}}(\widehat{I_p}-2)+2\zeta \big )\hat{{\mathbf {U}}}^\prime \\&\quad + \big ({\hat{\varOmega }}^2(\widehat{I_p}-1)+1+2\zeta \hat{\varOmega }{\mathbf {J}}\big )\hat{{\mathbf {U}}} \\&\quad + \gamma {{\hat{r}}}^2\hat{{\mathbf {U}}} = {\hat{m}}{\hat{e}}{\hat{\varOmega }}^2\begin{Bmatrix}1\\0\\\end{Bmatrix}. \end{aligned} \end{aligned}$$
(25)

The state-space representation is:

$$\begin{aligned} {{\dot{q}}}_1&=q_3\nonumber \\ {{\dot{q}}}_2&=q_4\nonumber \\ {{\dot{q}}}_3&= - \hat{\varOmega }(\widehat{I_p} - 2)q_4 - 2\zeta q_3 + \big ({\hat{\varOmega }}^2(1 - \widehat{I_p}) - 1\big )q_1 \nonumber \\&\quad + 2\zeta \hat{\varOmega }q_2 + {\hat{m}}{\hat{e}}{\hat{\varOmega }}^2 - \gamma {\hat{r}}^2q_1\nonumber \\ {{\dot{q}}}_4&=\hat{\varOmega }(\widehat{I_p} - 2)q_3 - 2\zeta q_4 + \big ({\hat{\varOmega }}^2(1 - \widehat{I_p}) - 1\big )q_2 \nonumber \\&\quad - 2\zeta \hat{\varOmega }q_1 - \gamma {\hat{r}}^2q_2 \end{aligned}$$
(26)

where \({\mathbf {q}}_{1-4}\) are \(\left[ u,v,{\dot{u}},{\dot{v}}\right] ^T\).

2.6 Natural frequency map

The natural frequency of the underlying linear system changes with rotor speed due to the gyroscopic couplings, or additionally through the dependence of the system parameters such as stiffness and damping on the rotor speed [13]. A natural frequency map, plotting the natural frequency with respect to the rotor speed, can help identify the critical speeds, where harmonic, subharmonic, or internal resonance occur [16]. The natural frequencies are shown in Fig. 3, from an eigenvalue analysis of the linear conservative parts of the equations of motion. Figure 3a shows a traditional representation of the Campbell diagram, where unsigned stationary frame frequencies are plotted. Figure 3b shows the form suggested in [35] where a signed frequency in the rotating frame is presented. Anticlockwise whirling is indicated by a positive value, whereas clockwise whirling is indicated by a negative value. This form is invaluable for highlighting drive speeds where combination resonances can occur.

Fig. 3
figure 3

The natural frequency map in a stationary frame, and b rotating frame. 1st engine order (EO) is also shown in stationary frame frequency map. Whirl frequencies are signed positive and negative for forward and backward whirl, respectively. The 3:1 and 2:1 internal resonance points are marked in the rotating frame

Full size image

3 Time simulation set-up

The explicit time marching was performed in [22] using the built in ode45 function. In this study, three cases of simulations were designed: constant initial condition simulations, sweep tests, and starting data generation for AUTO.

For constant initial condition simulations, a randomly generated stationary frame initial condition was used in the entire drive speed range to initiate the simulation. The initial conditions of the rotating frame simulations were derived from this initial condition using the transformations at \(\tau =0\) (nondimensionalised time) given by,

$$\begin{aligned}&{\mathbf {T}}_{{\hat{x}}{\hat{y}}\rightarrow {\hat{u}}{\hat{v}}} = \begin{bmatrix} 1&{}0&{}0&{}0\\ 0&{}1&{}0&{}0\\ 0&{}-\hat{\varOmega }&{}1&{}0\\ \hat{\varOmega }&{}0&{}0&{}1\\ \end{bmatrix}, \ \begin{pmatrix} {\hat{x}}\\ {\hat{y}}\\ \dot{{\hat{x}}}\\ \dot{{\hat{y}}} \end{pmatrix} = {\mathbf {T}}_{ {\hat{x}}{\hat{y}} \rightarrow {\hat{u}}{\hat{v}} } \begin{pmatrix} {\hat{u}}\\ {\hat{v}}\\ \dot{{\hat{u}}}\\ \dot{{\hat{v}}} \end{pmatrix} \end{aligned}$$
(27)
$$\begin{aligned}&{\mathbf {T}}_{\hat{\theta }\hat{\psi }\rightarrow {\hat{x}}{\hat{y}}} = \left[ \begin{matrix} 0&{}-1&{}0&{}0\\ 1&{}0&{}0&{}0\\ 0&{}0&{}0&{}-1\\ 0&{}0&{}1&{}0\\ \end{matrix}\right] , \ \begin{pmatrix} \hat{\theta }\\ \hat{\psi }\\ \dot{\hat{\theta }}\\ \dot{\hat{\psi }} \end{pmatrix} = {\mathbf {T}}_{\hat{\theta }\hat{\psi }\rightarrow {\hat{x}}{\hat{y}}} \begin{pmatrix} {\hat{x}}\\ {\hat{y}}\\ \dot{{\hat{x}}}\\ \dot{{\hat{y}}} \end{pmatrix} \end{aligned}$$
(28)
$$\begin{aligned}&{\mathbf {q}}_{{\hat{u}}{\hat{v}}} = {\mathbf {T}}_{ {\hat{x}}{\hat{y}} \rightarrow {\hat{u}}{\hat{v}} }^{-1} {\mathbf {T}}_{ \hat{\theta }\hat{\psi } \rightarrow {\hat{x}}{\hat{y}} }^{-1} {\mathbf {q}}_{\hat{\theta }\hat{\psi }} \end{aligned}$$
(29)

where the subscripts \(\hat{\theta }\hat{\psi }\), \({\hat{x}}{\hat{y}}\), and \({\hat{u}}{\hat{v}}\) represent nondimensional coordinates in their respective frames of reference. Simulations were run with damping ratio \(\zeta =0.01\), duration of nondimensional time, \(\tau _{\mathrm{end}}=3000\), and the relative tolerance, \(\hbox {tol} = 1e-7\). This relative tolerance value was adequate as lower tolerances did not show qualitative changes to the converged solution while increasing the computational time. Then, initial conditions were varied randomly to capture the responses of the system with different response types. The parfor looping provided by MATLAB was used to utilise 12 CPU threads in parallel.

For sweep tests, the drive speed, \(\hat{\varOmega }\), was swept in the range of 0.01–7.01 in the forward and backward directions. Each time a simulation of a particular drive speed finished, the last datum of the trajectory was selected as initial condition for the next \(\hat{\varOmega }\). To do this in the rotating coordinates, the last datum was transformed to the stationary frame using Eqs. (27) and (28), and then back to the rotating frame at the next drive speed using Eq. (29). Because of the need of one simulation to finish for the next to start, parallel computing was not available for the sweep tests. Different initial conditions were used until the full synchronous solution families were obtained without jumps to periodic orbits. Also, the asynchronous responses found using the constant initial condition simulations were swept in the same manner. For sweep tests, the same damping ratio \(\zeta \), simulation duration \(\tau _{\mathrm{end}}\), and relative tolerance values were used as mentioned above for the constant initial condition simulations. Lastly, the synchronous and asynchronous sweep test results were plotted to form a single bifurcation diagram.

Starting data generation for AUTO was created in the rotating frame for different \(\zeta \); the duration of simulation was increased as the damping ratio, \(\zeta \), was decreased, to ensure that the transient decayed within at most half the simulation duration. For stationary solutions in the rotating frame, the last data points in the time series generated by the time simulations were used. The time and state variable data that belong to one single period of the converged periodic orbits were stored in a .dat file to be used as starting periodic orbits in AUTO runs. More details of .dat file generation can be found in “Appendices”.

4 Numerical continuation set-up

4.1 Background

Compared to the time simulations of ode45, numerical continuation theory provides a different way to extracting the response of the numerical system from a previously known solution. The method [33] is based on the implicit function theorem which states that a regular solution point \(x_0=\left( u_0,\lambda _0\right) \) satisfying the system of equations of form \(G(u,\lambda )=0,\quad G:R^n\times R\rightarrow R^n\) has a continuous family of solutions provided that the Jacobian \(G_u(u_0,\lambda _0)\) is nonsingular and the function is continuously differentiable at this point. So, the derivative of the system gives:

$$\begin{aligned} G_u\frac{\mathrm {d}u}{\mathrm {d}\lambda }+G_\lambda =0 \ \rightarrow \ \frac{\mathrm {d}u}{\mathrm {d}\lambda }=-G_u^{-1}G_\lambda \end{aligned}$$
(30)

where \(G_\lambda \) is the derivative of the equations \(G(u.\lambda )=0\) with respect to \(\lambda \). However, due to the problem that arises when the Jacobian is singular at the so-called folds of a solution family, its arclength, s, is introduced as a free parameter, with an accompanying parametrisation:

$$\begin{aligned} \begin{aligned} G(u_1,\lambda _1)&= 0 \\ (u_1-u_0)\cdot {{\dot{u}}}_0 + (\lambda _1-\lambda _0){{\dot{\lambda }}}_0 - \varDelta s&=0. \end{aligned} \end{aligned}$$
(31)

In this way, the extended Jacobian, \(G_x=[G_u,G_\lambda ]\), is kept full rank. An initial guess for the next point of a solution branch can be written as [33],

$$\begin{aligned} (u_1^{(0)},\lambda _1^{(0)})=(u_0,\lambda _0)+\sigma (\varDelta u_0,\varDelta \lambda _0) \end{aligned}$$
(32)

where \(\sigma \) is a step length. A correction with a Newton method is then,

$$\begin{aligned} \begin{aligned}&\left[ \begin{matrix} G_u(u_1^{(\nu )},\lambda _1^{(\nu )}) &{} G_\lambda (u_1^{(\nu )},\lambda _1^{(\nu )}) \\ {{\dot{u}}}_0^T &{} {{\dot{\lambda }}}_0 \end{matrix}\right] \left\{ \begin{matrix} \varDelta u_1^{(\nu )} \\ \varDelta \lambda _1^{(\nu )} \\ \end{matrix}\right\} \\&\quad = - \begin{Bmatrix} G(u_1^{(\nu )},\lambda _1^{(\nu )}) \\ (u_1^{(\nu )}-u_0)\cdot {{\dot{u}}}_0+(\lambda _1^{(\nu )}-\lambda _0)\cdot {{\dot{\lambda }}}_0-\varDelta s \\ \end{Bmatrix}. \end{aligned} \end{aligned}$$
(33)

Once the convergence is achieved, the next direction vector is found from the derivative of Eq. (31),

$$\begin{aligned} \begin{bmatrix} G_u^1 &{} G_\lambda ^1 \\ {{\dot{u}}}_0^T &{} {{\dot{\lambda }}}_0 \\ \end{bmatrix} \begin{Bmatrix} {\dot{u}}_1 \\ {\dot{\lambda }}_1 \end{Bmatrix} = \begin{Bmatrix} 0 \\ 1 \end{Bmatrix}. \end{aligned}$$
(34)

The new direction vector is rescaled by normalising its norm \({||}{\dot{x}}_1{||}^2={||}{\dot{u}}_1{||}^2+\lambda _1^2=1\).

A boundary value problem’s continuation procedure can be used [33] to follow periodic solution branches of an n-dimensional ODE, \(y^\prime =Tf(y,\lambda )\), with the derivative being with respect to normalised time, t with \(0\le ~t\le ~1\). The ODE problem then needs to be augmented with a periodicity condition, \(y(0)=y(1)\), and a phase condition, e.g. \(p(y(0),\lambda )={{\dot{y}}}_k(t=0)=0\). Accompanying these are the statements of the period and continuation parameter, \(\lambda \) not changing for a certain periodic solution: \(T^\prime =0\) and \(\lambda ^\prime =0\). The continuation of a branch of periodic solutions then take the form,

$$\begin{aligned} \begin{pmatrix} y \\ T \\ \lambda \end{pmatrix}^\prime = \begin{pmatrix} Tf(y,\lambda )\\ 0 \\ 0 \end{pmatrix} , \quad {\tilde{r}} = \begin{pmatrix} y(0)-y(1) \\ P(y(0),\lambda ) \\ p(y(0),\lambda ) \end{pmatrix}=0\nonumber \\ \end{aligned}$$
(35)

where P is a parametrisation such as arclength parametrisation of Eq. (31), and p is a phase condition. The discretisation of the ODE part will ultimately result in an algebraic equation of the form \(F(y,T,\lambda )=0\) [33].

Discussion of stability is covered well in the text books, for example [33]. AUTO, the numerical continuation software employed in this study, uses a method proposed by Fairgrieve and Jepson [11], where the Floquet multipliers are calculated without explicit calculation of the monodromy matrix.

4.2 AUTO files

The AUTO-07p open source numerical continuation software [9] was used to apply pseudo-arclength continuation. This software is operated through a Unix Command Language User Interface (CLUI), and can be written in a Python CLUI which supports Python 3 syntax. AUTO set-up basically includes filling three dedicated files: equations file, constants file, Python CLUI files.

4.2.1 Equations file

The equations file requires information on the equation of motion, as well as the starting conditions of parameters and a regular solution on a branch of solutions that is to be continued.

The rotating frame equation of motion (Eq. 26) was used since in the rotating frame the forcing becomes constant, which facilitates the analysis by making the system of equations autonomous. Also, asynchronous solutions appear periodic (although quasi-periodic in the stationary frame, as the Poincare map has closed loops). Therefore, the synchronous solutions are sought with the stationary continuation solver \((\hbox {IPS}=1)\); whereas asynchronous solutions are found by employing periodic continuation solver of AUTO \((\hbox {IPS}=2)\).

Amplitude was defined to be a solution measure in the PVLS subroutine of the equations file. The amplitudes of stationary and periodic solutions were calculated as,

$$\begin{aligned} \begin{aligned} \hbox {Amplitude for IPS:1}&= \sqrt{q_1^2+q_2^2} = \sqrt{{\hat{u}}^2+{\hat{v}}^2} \\ \hbox {Amplitude for IPS:2}&= \max {\sqrt{q_1^2(\tau )+q_2^2(\tau )}}\\&= \max {\sqrt{{\hat{u}}^2(\tau )+{\hat{v}}^2(\tau )}} \end{aligned} \end{aligned}$$
(36)

where the nondimensional time, \(\tau \), is defined between \(0\le \tau \le 1\).

The constants file is where the user mainly specifies the continuation software’s solution control parameters. The main adjustments of these constants concern those that define the accuracy. More detail can be found in “Appendices”.

4.2.2 Continuation solution process

The Python CLUI file contains the chronological flow of runs. The complete content of this file can be written in the AUTO CLUI shell. However, the usage of Python CLUI eases the use of the software by a complete view of the runs. The flow of runs is shown in Table 1. The references to the single- and double-loop families in Table 1 indicate the shape of the orbits corresponding to the 3:1 and 2:1 internal resonances, respectively, as seen in Fig. 4.

Table 1 Flow of continuation runs in Python CLUI file of AUTO set-up
Full size table
Fig. 4
figure 4

One sample from each of the periodic solution families, referred to as a single- and b double-loop solutions corresponding to 2:1 and 3:1 internal resonance conditions in the rotating frame

Full size image

Damping ratio, \(\zeta \), was selected as the main continuation parameter because the level of damping might explain the rotor speed range where the asynchronous solutions could be expected, e.g. earlier in the range and more likely to appear particularly with low damping (see in Sect. 5). In Table 1, A, B and C are different ways the continuation was carried out:

A::

A first continuation run with \(\zeta =0.01\) is performed through homotopy of the \({\hat{m}},\gamma \) and \(\hat{\varOmega }\) from zero, initialised with the known analytical values for synchronous solutions and simulation results for periodic solutions. Then, continuation in \(\zeta \) is performed from a selected point of the solution family. After, from the data sampled at various \(\zeta \) during this \(\zeta \)-continuation, new solution branches are found with continuation in \(\hat{\varOmega }\).

B::

Only time simulation starting data is used in \(\hat{\varOmega }\)-continuation for different \(\zeta \) values

C::

Folds in \(\hat{\varOmega }\) are continued in \(\zeta \) to trace the effect of the change of \(\zeta \).

and the cases 1, 2, and 3 are,

1::

Synchronous solution branch continuation

2::

Single-loop branch continuation

3::

Double-loop branch continuation

5 Results and discussion

Results of the numerical continuation are presented first; then, in Sect. 5.2, the time simulation data used to validate these results and provide starting solutions is discussed.

5.1 Numerical continuation

Figure 5 shows the bifurcation diagram in nondimensionalised amplitude versus rotor speed, for all of the solutions from the different \(\zeta \) values investigated. Labels SS, DL, and SL in the figure refer to the solution families of synchronous, double-loop asynchronous, and single-loop asynchronous solutions, respectively. Synchronous and asynchronous solutions are also referred to as stationary and periodic as only the rotating frame equations of motion were used in AUTO.

For the lower \(\zeta \), the SS branch exceeds the limit of the plot. A close-up view of the SS branch is shown in Fig. 5d. For the increased damping values, the synchronous peaks merge and form a stiffened nonlinear response peak. As the damping is further increased, the unstable solutions, which are shown with dashed lines, are lost. At higher drive speeds the amplitude of SS branch decreases very slightly with increased damping.

Fig. 5
figure 5

The bifurcation diagram generated from AUTO using the autonomous rotating frame equations of motion. Both stationary (synchronous) and periodic (asynchronous) solutions are included. The zoom views at b and c show the change in double-loop and single-loop solution families with \(\zeta \), respectively. The zoom view in d shows the synchronous region. The zoom view in e shows the cyan-coloured single-loop periodic solutions with \(\zeta =0.082\). In the figure, ’SS’, ’DL’, and ’SL’ stand for synchronous, double-loop, and single-loop solution families, respectively. The legends show \(\zeta \) values of periodic solutions in b and c, and of synchronous solutions in d. (Color figure online)

Full size image
Fig. 6
figure 6

Zeta sampling with continuation in zeta of a single- and b double-loop periodic solutions in the rotating frame. Each “U” symbol indicates a different sample that was continued to form a periodic solution family

Full size image

Figure 5 also shows that there is a region of nondimensionalised rotor speed between \(\hat{\varOmega }=1.5-2.3\) where jumps would only happen between the synchronous solutions. Chipato et al. [6] and Zilli et al. [39] also documented this phenomenon in their overhung models with frictionless contact definitions.

Fig. 7
figure 7

Zoom view at the bifurcation diagram where stationary and single-loop periodic solutions meet in the lower-\(\zeta \) extreme. Turquoise and purple lines show \(\zeta =1e-4\) and \(1e-5\), respectively, and “SS” stands for synchronous solutions. The stacked lines below the 0.413 amplitude level show the synchronous branches where higher amplitudes correspond to lower \(\zeta \). (Color figure online)

Full size image
Fig. 8
figure 8

Fold continuation in order to decrease the severity of eccentricity parameter on double-loop periodic solution family with \(\zeta =1e-5\)

Full size image

As damping increases, the fold at the lower end of both the DL and SL solution branches moves higher in terms of both amplitude and drive speed, with the SL branch moving beyond the range of drive speeds considered at higher damping values. The various \(\zeta \) value sampling of periodic solution families was done through the continuation of the \(\zeta =0.01\) orbits in AUTO, for Methods A2 and A3 (see Table 1). The results are given in Fig. 6. The higher damping ratio limit to the existence of SL and DL are \(\zeta =0.0828\) and \(\zeta =0.0117\), respectively. This shows that the increase in damping in the system would first eliminate DL solutions. SL solution families show that the increased \(\zeta \) shrinks and eliminates this type of periodic solution family out of the investigated drive speed region (see Fig. 5e). This can be observed in Fig. 5e where the brown, dark turquoise, and cyan coloured branches correspond to \(\zeta \) of 0.050, 0.075, and 0.082, respectively. Note that in this figure the dark blue dashed line shows the fold continuation (Method C2) of the SL solutions. Increased \(\zeta \) for DL branches, however, kept these solution families as isolated amplitude circles within the drive speed region of interest until the damping ratio is above the 0.0117 limit.

Fig. 9
figure 9

a Bifurcation diagram for the angular equations of motion with initial condition [0.0; 0.6; \(-1.246\); 0.0]. b and c show the double-loop response from the perspective of stationary and rotating frame of references, respectively. d and e show single-loop response in stationary and rotating frame, respectively. The orange dots in b and d show the stroboscopic Poincare sampling at the drive speed after the convergence. (Color figure online)

Full size image
Fig. 10
figure 10

Bifurcation diagrams from a constant initial condition simulation results and b stepped-sine sweep tests on the synchronous and asynchronous responses

Full size image

Figure 5b and c shows that the decrease in damping extends the left-hand side of the periodic solution families downwards, closing into the synchronous solutions at \(\hat{\varOmega }=2.31\) and \(\hat{\varOmega }=3.46\) for DL and SL solution families, respectively. Figure 5b also shows the locus of folds as zeta is freed in addition to rotor speed. At the place where the periodic and stationary solutions meet, the periodic solution branches become quite sharp, making it difficult to perform continuation. As a result, the SL solution family for \(\zeta =1e-5\) converged only with increased number of Newton–Chord correction iterations. The results are seen in Fig. 7.

The Campbell diagram of the linear, undamped system (Fig. 3) shows possible locations of the internal resonances at \(\hat{\varOmega }=1.81\) for 4:1, \(\hat{\varOmega }=2.18\) for 3:1, \(\hat{\varOmega }=3.33\) for 2:1, \(\hat{\varOmega }=5.84\) for 3:2. Therefore, DL solution families could be related to 3:1 internal resonances, and SL periodic solutions could arise from 2:1 internal resonances. To investigate the relationship of the 3:1 internal resonance to DL families, the fold with \(\zeta =1e-5\) was continued in eccentricity as an extra continuation parameter to the drive speed. The results are presented in Fig. 8. A decrease in eccentricity, \({\hat{e}}\) from 0.353 to 0.092 made DL solution family’s left tip frequency to go down to 2.19, much closer to the 3:1 frequency of 2.18 from the Campbell diagram (Fig. 3). As the forcing is decreased, the amplitude reduces. By reducing the amplitude, the stiffening effect of amplitude is also decreased, so that the underlying frequencies that are in internal resonance move closer to their linear values. And as a result, the solution moves closer to the critical drive speed predicted using the Campbell plot (Fig. 3). In the vicinity of these points, any synchronous response would take very little disturbance to move onto and beyond the unstable part of the periodic solutions, possibly converging onto the stable limit cycle. Nonetheless, the brute force bifurcation diagrams do not show the periodic solution families at the place they are expected, because the time simulations only show \(\zeta =0.01\) data, which appear with higher amplitudes, and hence with higher nonlinearity. In the limit of zero out of balance force and zero damping, the fold reached the intersection point on the linear Campbell diagram. In fact, as forcing and damping increase, the solutions become more isolated as the fold point increases in amplitude. Even in that asymptotic condition, the stiffening with amplitude would cause these points to diverge from what is expected from the Campbell diagram.

Fig. 11
figure 11

Validation of single-loop orbits with \(\zeta =0.01\) from AUTO and MATLAB. The bottom row shows the trajectory history (from left to right) during the transition from UPO to SPO at \(\hat{\varOmega }=3.47\). The green squares are stroboscopic Poincare sampling of the data at the rotor speed. Floquet multipliers indicate the stability. (Color figure online)

Full size image
Fig. 12
figure 12

Validation of single-loop orbits with \(\zeta =0.01\) from AUTO and MATLAB. The bottom row shows the trajectory history (from left to right) during the transition from UPO to SPO at \(\hat{\varOmega }=3.47\). The green squares are stroboscopic Poincare sampling of the data at the rotor speed. Floquet multipliers indicate the stability. (Color figure online)

Full size image

These results compare to those of contact nonlinearity studies. For example, in Ref. [35] on a similar 2-dof model with piecewise-linear stiffness contact model, the periodic solutions were observed to start in a similar range after \(\hat{\varOmega }=3.5\), with similar amplitudes for a 2:1 internal resonance. Zilli et al. [39] also observed the intermittent contact motions after \(\hat{\varOmega }=3.4\).

No asynchronous response under contacting level amplitude can be seen in the studies that model a clearance as the underlying system is linear. However, the cubic stiffness model here removes this limit and bridges the gap between the synchronous and asynchronous responses. This possibly reveals the transition from synchronous to asynchronous response, as discussed above. As a result, the stiffening effect that is seen in the contact-featuring studies is available at all amplitudes in the cubic stiffness continuous stiffness model.

The present isotropic cubic stiffness model did not produce period doubling in the continuation runs, which is an identified route to chaos in systems with stiffness asymmetry [18, 37]. Asymmetry due to gravity also resulted in the chaos in [6]. Li and Païdoussis [21] reported the quasi-periodic route to chaos in the presence of rubbing friction, which is not a part of the model in this study, although we observe quasi-periodicity. These indicate that the isotropic cubic stiffness model used here does not present a way to explain the chaotic behaviour. Thus, the numerical simulations with various initial conditions did not result in chaos. However, a homotopy parameter can be used to include the friction, gravity, contact (e.g. via hyperbolic tangent functions), making the cubic stiffness a base to explore these responses.

Contrary to the systems with impact definitions in the contact [8, 23], the smoothness of the trajectories of the cubic stiffness model is not broken. Also the effect of contact dissipation cannot be seen in the present model, on which Shaw et al. [34] noted that stiffer, dissipative contact definitions eliminated most of the sustained intermittent bouncing cycles.

5.2 Time simulations

5.2.1 Constant initial condition

The simulations with constant initial conditions in the entire rotor speed range yielded different bifurcation diagrams. The 12 thread CPU parallelisation decreased the simulation time by 5 times. Figure 9 shows the response of the system to the initial conditions [0.0; 0.6; \(-1.246\); 0.0] in angular coordinates (physical coordinates). Asynchronous solutions appear quasi-periodic in the stationary frame and periodic in the rotating frame. It can be noted here that the quasi-periodic solutions disappear at higher drive speeds, while easily achieved in the 3.5–5.0 region.

The double-loop orbits were only discovered as the initial condition was increased in amplitude (Fig. 9a and b), With the greater initial energy due to higher magnitude position and/or speed, the orbits settle more readily onto the asynchronous orbits, which were recorded in the transient study in Ref. [39].

In order to extract the most significant features of the complete bifurcation plot, a number of random initial conditions in the order of 100 was selectively merged and is illustrated in Fig. 10a. The stable synchronous high- and low-amplitude solution families, as well as two distinct asynchronous solution families can be identified. The difficulty of finding the double-loop orbits was noted in the process.

5.2.2 Sweep

Figure 10b shows the forward and backward sweep test results as a bifurcation diagram for both synchronous and asynchronous responses. Once the different solution families are identified with constant initial condition simulations (Fig. 10a), it is more likely for sweep tests to discover a given family of solutions (Fig. 10b). Also note that in backward sweep and periodic solution sweep tests, the unstable parts of the solution families are not seen, which are clearly visible in the AUTO numerical continuation results.

In the sweep tests, the next increment in the same family of solutions converges with less transient motion, thereby making the extraction of system response quicker than each time simulation in the case of constant initial condition simulations. However, the advantage of multi-core processors cannot be utilised because in sweep test the last point of the trajectory of the previous step is required in the next step. Nonetheless, in terms of computational speed, time simulations cannot compare to numerical continuations: for example, the computation of Fig. 5 took under 10 min, while time for the simulations easily reached an hour.

With these constant initial condition and sweep test results, it is clear that the numerical continuation is computationally more efficient and reveals more of the system response characteristics.

5.2.3 MATLAB: AUTO validation

The single-loop orbits were compared to those of AUTO for validation of both simulations and the continuation with each other. An unstable single-loop periodic orbit at \(\hat{\varOmega }=3.47\) with \(\zeta =0.01\), was given to MATLAB as initial condition for ode45 simulations. Results are illustrated in Fig. 11. The trajectory in the ode45 simulations stayed on the unstable periodic orbit (UPO) for the first few cycles, then diverged from this unstable periodic orbit to converge onto the stable periodic orbit (SPO). This diverging speed is much higher for the higher amplitude case (Fig. 12), where also the greater rotation of the orbit is visible. A close match was obtained between the two solution methods with respect to the stable and unstable orbits. As seen in Figs. 11 and 12, the absolute value of one Floquet multiplier is greater than unity for the unstable orbits.

Table 2 AUTO initial conditions for rot-EOM
Full size table

6 Conclusion

This study considered the responses of an overhung rotor with an isotropic cubic stiffness. Direct application of continuation method on the equations of motion was set up and validated with time simulations. The following conclusions were drawn:

  1. 1.

    It was shown in the present work that the previously reported bouncing cycle phenomenon for the models of hard impact and compliant snubbing rings can also be observed with a geometric nonlinearity of cubic stiffness replacing contact interaction definitions.

  2. 2.

    If a discontinuous stiffness was defined, then the bouncing cycles would not be present under the clearance level amplitude, due to the rotor not experiencing nonlinearity below this level. On the contrary, nonlinearity is present in all amplitude levels for the cubic nonlinearity; thus, the bouncing cycle amplitudes could get very close to the synchronous solutions.

  3. 3.

    Depending on the level of damping, the periodic orbit region can expand and shrink. Especially for higher \(\zeta \) damping ratio values, the periodic orbits’ region vanishes from the nondimensionalised rotor speed region of interest.

  4. 4.

    As damping tends to zero, the lower limit of drive speeds that see these responses tends towards the critical speed predicted with linear frequencies in the natural frequency map. Decreasing the eccentricity level (or out of balance mass) makes this match closer.

  5. 5.

    As damping tends towards zero, the periodic solution family approaches the synchronous orbit. It might be inferred that the gap closes asymptotically as \(\zeta \) and oscillation amplitude approaches zero.

  6. 6.

    This indicates that this region requires very little to no disturbance to jump to the periodic orbits because of the closeness of the unstable side of periodic solutions to the stable synchronous solutions.

  7. 7.

    As the results were compared to the previous studies featuring nonsmooth contact models, the synchronous solution jump, and the nondimensional rotor speed where periodic solutions first appeared were reproduced with the cubic stiffness model adopted here.

The application of numerical continuation theory directly to the system equations proved practical and accurate in this study when compared to the brute force bifurcation diagrams. The type of low order model studied, with a simple smooth stiffening nonlinearity, is also relevant to many rotating machines with flexible shafts, that may exhibit bend stretch coupling, or overhung rotors that may show nonlinear displacements in the case of large out of balance force.

In future work, we are looking to extend the work to consider smooth approximations to contact phenomena that are seen in wide range of industrial applications, particularly where bearings are present. A further aim is to explore systems with multiple degrees of freedom. If these aims can be achieved, numerical continuation will be proven as a highly efficient tool for exploring complex phenomena in industrial rotating machinery.