Experimental characterisation of asynchronous partially contacting motion in a multiple-degree-of-freedom rotor system Mechanical Systems and Signal Processing

Recent theory has predicted the onset of asynchronous bouncing motion at speeds beyond those of internal resonance in multi-degree-of-freedom rotating systems with intermittent contact. This paper provides the ﬁrst attempt to experimentally validate the theory. Vibrations incorporating rotor-stator contact are recorded from a vertically ounted rotordynamics test rig comprising two rigid shaft-disk assemblies that are axially joined by a bellows coupling. The upper rotor has an elastic bearing whereas the lower one is free but has a snubber ring located on its lower shaft, with clearance about 2.5% of the entire system’s length. Nonlinear vibrations are excited by a small measured eccentricity of about 20% of the clearance. Measurements are taken by a wireless accelerometer and processed to produce an experimental bifurcation diagram, for small steps in the rotor speed, allowing the transients to decay. Evidence is found of bistability between quiescent and violent motion over a wide range of rotational speeds, including those representing both fundamental and internal resonances. The results are shown to qualitatively match numerical simulations from a differential equation model that incorporates rigid impacts. bottom of the rotor corresponding to a disk and the sensor and at the top, the rotor is attached to a motor without any bellows couplings. The preliminary test was conducted at different rotational speeds of the motor. The mean acceleration in x and y in the rotating frame was obtained, to calculate the radial position of the sensor from the axis of rotation. The radial displacements calculated were r x ¼ (cid:2) 0 : 012 m and r y ¼ (cid:2) 0 : 001 m. The acceleration in x and y due to the radial position of the sensor are subtracted from the acceleration response of the sensor in each test conducted with the double rotor system. The acceleration in the x and y directions are calculated by multiplying r x and r y by the rotational speed of the motor squared in each test. The resulting acceleration in the rotating frame is transformed to the stationary frame by the transformation matrix T e A Butterworth high pass ﬁlter of order 2 with a cut-off frequency at 1 Hz is used to eliminate any drift when the acceleration signal is integrated. The use of a high-pass ﬁlter has the disadvantage of reducing the synchronous whirl solutions to zero because the synchronous whirl is stationary in the rotating frame. The ﬁltered acceleration signal in the stationary frame is then integrated twice to obtain the x 2 and y 2 displacements.


Introduction
Rotating machines are unlike other mechanical vibration systems because of the large rotational kinetic energy that can be transferred into lateral or torsional vibration modes whenever the rotor is out of balance. Reducing these vibrations is often the main goal of the design and control of rotating machines.
Extensive recent research has explored, both experimentally and through mathematical models, the nonlinear vibrational dynamics of rotating machines, see e.g. [1][2][3][4][5][6][7][8][9][10][11][12][13][14], including complex phenomena such as quasi-periodic motion and even chaotic motion observed by a number of other authors, see for example [15,16]. The typical mechanisms that allow energy transfer to vibration modes from the rotating mass is the presence of centrifugal forces due to out-of-balance or out-of-alignment eccentricity. These abnormalities may arise due to manufacturing and assembling deficiencies, or due to operational conditions and wear. Finite amplitude behaviour is typically then highly sensitive to resonance between vibrational modes and the rotor drivespeed, and is exacerbated by the presence of various kinds of nonlinearity. Perhaps the most common causes of nonlinear behaviour in any rotor are a wide variety of effects due to the dynamics of the bearings; see for example [1,17]. Here the specific source of bearing nonlinearity can vary from squeeze film effects [18][19][20][21][22], the so-called dry whip of journal bearings [23], to magnetic bearing singularity effects [24].
Another major class of nonlinearity in rotor dynamics is geometric nonlinearity of various kinds, see e.g. [25]. For example, Luczko [26] included the so-called von-Karman nonlinearity, curvature effects, and large displacement gyroscopic and shear effects. Moreover, Ferjaoui et al. [27] analysed the effect of a cracked shaft on the onset of a whirl instability.
A more violent source of geometric nonlinear behaviour is the presence of free masses that are able to rotate in order to form an automatic dynamic balancing mechanism, see [28,29] and references therein for comprehensive bifurcation analyses and comparison with experimental results. Similarly, van de Wouw et al. [30] developed an experimental test rig to understand the interaction between torsional ball bearing dry-friction effects used as a form of stabilizer. Basically all of these systems, which are able to balance a rotor beyond its fundamental resonance, also suffer from the possibility of mode-locking into large-amplitude sub-harmonic whirling limit cycles, similar to those described by van der Heijden [31,32] for a Jeffcott rotor. This latter work [32] was one of the first papers to consider the dynamics of rotor-stator contact from the point of view of the geometric theory of dynamical systems, including the concept of Arnol'd tongues.
More generally, and indeed the theme of this paper, probably the most extreme and most damaging form of nonlinear behaviour in rotordynamic systems is due to rotor-stator contact, see e.g. [33,16,34] for comprehensive studies. Such problems often occur due to interactions when overhung rotors, such as aircraft engines, come into contact with so-called snubber rings that are designed to limit vibration amplitude [35,36]. Such problems are also relevant for micro-electromechanical system (MEMS) rotors [37]. Furthermore, rotating flexible shafts in drillstrings [38,39] can present complex interactions with the borehole such as impacts, when combined with bit-bounce and stick-slip. Magnetic bearing systems also require good models for this behaviour when dealing with problems such as touch down [40,41].

Latin Letters a
Position of the stator measured from dis k 2 to the stator. b Length to the centre of disk 1. c Radial clearance of stator. C Damping matrix. c r1 ; c r2 Linear angular damping of bellows couplings 1 and 2. e 1 ; e 2 Eccentricities of disks 1 and 2. d 1 ; d 2 Diameters of disks 1 and 2. F i Excitation vector due to impacts between rotor 2 and the stator. F e Excitation vector due to the mass imbalance of each disk. G Gyroscopic matrix. g Acceleration due to gravity. h Equivalent radial clearance. K Stiffness matrix. k s Linear stiffness of the stator. k r1 ; k r2 Linear angular stiffness of bellows couplings 1 and 2.  Greek Letters X Rotational speed of the motor. w 1 ; w 2 Rotational coordinates of disks 1 and 2 about the y axis. h 1 ; h 2 Rotational coordinates of disks 1 and 2 about the x axis. f 1 ; f 2 Modal damping ratios due to bellows couplings 1 and 2.
The nonlinear behaviour associated with rotor-stator contact can largely be characterised as being due to one of two effects; namely sustained frictional contact or 'rub', and intermittent contact or 'bouncing.' Studies focused on rub include [42][43][44][45] which have shown that both forward and backward permanent-in-contact whirl motion is possible depending on whether the rotor slips or sticks around the stator. Under rub conditions, torsional effects become important [46] and can give rise to Hopf bifurcations and stick-slip behaviour [47].
In contrast, the dynamics associated with bouncing are potentially more violent than pure rub. Complex and hard-toexplain responses due to rotor stator contact have been reported over many years, and perhaps some of the earliest work highlighting the complexity of this problem was reported by Ehrich [48,49] and at a similar time by Black [50]. An excellent overview of similar works and phenomena can be found in [51], and a review with many more works concerning turbomachinery can be found in [52]. More recently authors began to consider the rotor stator-contact in relation to chaos and bifurcation theory [53,54]. Furthermore, non-smooth approaches are becoming widespread [55][56][57], with the latter reference highlighting a novel form of bifurcation that can occur with this approach. The role of internal resonance has also been considered, with Inoue and Ishida showing that a 1:1 resonance between forward and backward whirls has an important effect on shaping the response of contacting systems [58].
An important classification is made between synchronous and asynchronous bouncing cycles [50,56]. The former involve frequencies with integer ratios to the drive speed, and occur in rotor systems with anisotropic bearings, or other anisotropic effects such as significant gravity acting laterally to the shaft. On the other hand, asynchronous cycles are more complicated because response frequencies seem unrelated to the drive frequency, and they also can occur at drivespeeds removed from direct resonances of the system. (Examples of these phenomena can be found throughout the history of the literature, for example [23,59].) The work of Zilli et al. [10] was one of the first studies to propose a theory, backed by simulations for the onset of asynchronous bouncing orbits with numerical results. A synchronicity 1 condition proposed by Zilli et al. was later shown by Shaw et al. [60] to be equivalent to an internal resonance condition when viewed in the rotating frame. Bouncing motion that appears as a steady single-impact-per-period limit cycle within the rotating frame is observed to suddenly appear at drive speeds just beyond such internal resonances. This motion is observed to co-exist with quiescent non-contacting whirl motion. Nevertheless, in a later paper [61] it is shown how such highly nonlinear bouncing motion can be predicted semianalytically using a harmonic-balance and time-domain method provided the stator is not too stiff. Note that the analysis in [60] is quite general and can be applied to any multiple degree-of-freedom rotor system, and indeed provides the primary inspiration for the present study.
In view of these recent results, the aim of this study is to provide some experimental verification of the theory of asynchronous 'bouncing' orbits due to internal resonance. This is supported with an experimental bifurcation analysis to determine at which drive speeds the intermittent contacts occur. The experimental rig studied in this research comprised of two shafts with asymmetrically mounted non-identical disks positioned at a given length on each shaft, and both shafts were axially joined by a bellows coupling. The upper shaft was attached to an elastic bearing and the lower shaft is free at the bottom end but has a snubber ring, similar to the main example in [60]. The ultimate objective of this work is to gain understanding of these bouncing orbits with experimental results to help engineers when designing rotating machines. A proper design of a rotating machine may prevent intermittent contacts from occurring. The presence of these intermittent contacts in rotating systems may lead to a catastrophic failure of the system. The rest of this paper is outlined as follows. Section 2 describes the experimental test rig, and Section 3 goes on to describe its mathematical model and analyses its natural frequencies. Section 4 then presents experimental results, which are explained in the light of new simulation results. Finally, Section 5 draws conclusions and suggests avenues for future work.

Experimental test rig
The experimental test rig is first shown. Then a description is given on how the system parameters were estimated, before describing a linear frequency analysis.

Description of test rig
The experimental test rig is shown in Fig. 1(a); a vertically mounted rotor consisting of two disks and two rotor shafts connected by isotropic bellows couplings. Impacts are allowed to occur with a stator, taking the form of a concentrically mounted snubber-ring stator on the second (lower) shaft. The stator clearance is circular and symmetric with the same centre as both rotors. Note that, apart from its coupling to the first rotor, the second rotor is free and has no bearings, as can be seen in the zoom in Fig. 1(b). A flexible coupling was used to transfer the rotating motion of the motor to the system.
A brushless electric motor with integrated electronics BLDC58-50L is positioned at the top end of rotor 1, to provide a fixed rotational drive speed. This speed is controlled by a voltage input, which is supplied by a National Instruments (NI) cRio-9024 data acquisition system with output module NI-9263. The BLDC58-50L has an embedded system that monitors the speed of the motor; the signal is acquired using an input module NI-9215 on an FPGA architecture and the LabView software. To measure the vibrations of the stator during impact a PCB triaxial accelerometer 356A02 is attached to the top of the stator and the signal is acquired using an input module NI-9234. To measure the lateral motions at the end of rotor 2 a Lord Microstrain triaxial wireless accelerometer G-Link-200-8G is attached and the signal is acquired using Lord Microstrain WSDA-Base-101-LXRS gateway in conjunction with SensorConnect software for the wireless sensor configuration and data acquisition.

Preliminary experimental testing
Initial tests were conducted to determine the basic parameters of the rig and to characterise its underlying linear responses. First, by using a single free pendulum-like lateral motion and using a logarithmic decrement method, the stiffness and damping properties of the couplings were estimated to be k r1 ¼ k r2 ¼ 1:20 Nm and c 1 ¼ c 2 ¼ 0:258 kg/s respectively. From the damping coefficients given above, the first modal damping ratio f of the system was calculated to be f ¼ 0:05. The modal damping ratio was estimated by using the first natural frequency of the system obtained from experimental testing.
Next, the method of Costache et al. [62] was adopted to estimate the eccentricity of the disk plus triaxial accelerometer attached at the end of a vertical slender shaft. Individual tests were conducted with a single shaft and mass, attached to the motor as shown in Fig. 2, where the acceleration due to gravity acts in the z direction. The inertial reference frame xyz is positioned at the origin O, where the motor provides the rotational speed to the system. The top end of the rotor was attached to the motor which was running at a variety of different rotational speeds X in the range 70.0-540.0 rpm. The eccentricity of a disk plus triaxial accelerometer (which has dimensions of length) was estimated by The derivation of Eq. (1) is shown in A. In Eq. (1),His the angular displacement of the whirling observed at drivespeed X.
An equivalent disk diameter d and thickness v was estimated to be d ¼ 0:038 m and v ¼ 0:063 m. The equivalent disk diameter and thickness were obtained by taking an average of the disk and attached triaxial accelerometer dimensions. The total mass was measured to be m ¼ 0:385 kg. The length L for the single rotor used in this test was L ¼ 0:14 m. The stiffness k r and modal damping ratio f used during the test was k r ¼ 1:20 Nm and f ¼ cr 2mxe ¼ 0:05, respectively. The fundamental natural frequency of the system is represented as x e and is defined in equation (2).
Using this method, the estimated eccentricity of the disk including a triaxial accelerometer was approximately 0:002 m. This eccentricity was due to the geometrical imperfections of the rotor, assuming it takes the form of the disk being slightly off centre. The estimated eccentricity corresponds to disk 2 in the test rig. The eccentricity of disk 1 in the test rig, which consists only of a solid disk, was assumed to be 0 m.
The stator stiffness k s was defined to be the bending stiffness of the combined assembly of the snubber ring plus the four steel bars shown in Fig. 1(a) which are rigidly attached at their top end to a steel plate and at their bottom end to the snubber. This stiffness was obtained in two ways by measurement and by calculation. First, the assembly was allowed to vibrate freely after being hit with a hammer, and lateral vibrations were measured using the stator accelerometer. The measured fundamental natural frequency was x s ¼ 12:50 Hz. Via calculation; the steel bars each have a length L s ¼ 0:4125 m, diameter d p ¼ 0:003 m, modulus of elasticity of steel E s ¼ 210 Â 10 9 N/m 2 , and the total measured mass of the assembly was m s ¼ 0:093 kg. Using standard formulae, the stator bending stiffness was estimated to be k s ¼ 571:00 N/m. Using the measured mass m s , the calculated natural frequency was x s ¼ 12:47 Hz, giving an error of 0:23% when compared with the measured x s . Thus, the calculated bending stiffness was accepted.
Note that the operating frequency range of the motor speeds in the experimental testing was up to 10 Hz, which is below the natural frequency of the stator assembly. Thus, when deriving a mathematical model of the rig (see Section 3) it is acceptable to neglect the inertia of the stator.
In summary, the parameters of the experimental test rig are listed in Table 1.

Testing protocol
The experimental tests consisted of running the motor at a wide range of different rotor speeds from about 120 rpm to 600 rpm (2 Hz-10 Hz). Each test had a duration 60.0s which was found in practice to be sufficient time for transient behaviour to decay. At the end of each test, the rotor was stopped and a new test initiated. Two different versions of each test were run, with different initial conditions. Condition A consisted of allowing the rotor to be in its free-hanging equilibrium position, with no applied displacement or velocity. Condition B consisted of initially displacing rotor 1 and rotor 2 by hand and then letting go when the rotor reaches full speed, so that the shafts are both at an angle of approximately 0.05 rad to the vertical, with zero angular velocity. This particular angular displacement was chosen so that rotor 2 was in gentle contact with the stator.
To determine the displacements of x 2 and y 2 in the stationary frame from the accelerometer response measured in the rotating frame, a methodology was developed to transform the measured acceleration from the position of the sensor to the sensor axis of rotation. The preliminary testing configuration consisted of a single rotor with a discrete mass attached at the bottom of the rotor corresponding to a disk and the sensor and at the top, the rotor is attached to a motor without any bellows couplings. The preliminary test was conducted at different rotational speeds of the motor. The mean acceleration in x and y in the rotating frame was obtained, to calculate the radial position of the sensor from the axis of rotation. The radial displacements calculated were r x ¼ À0:012 m and r y ¼ À0:001 m. The acceleration in x and y due to the radial position of the sensor are subtracted from the acceleration response of the sensor in each test conducted with the double rotor system. The acceleration in the x and y directions are calculated by multiplying r x and r y by the rotational speed of the motor squared in each test. The resulting acceleration in the rotating frame is transformed to the stationary frame by the transformation matrix T e T e ¼ cos Xt À sin Xt sin Xt cos Xt

5: ð3Þ
A Butterworth high pass filter of order 2 with a cut-off frequency at 1 Hz is used to eliminate any drift when the acceleration signal is integrated. The use of a high-pass filter has the disadvantage of reducing the synchronous whirl solutions to zero because the synchronous whirl is stationary in the rotating frame. The filtered acceleration signal in the stationary frame is then integrated twice to obtain the x 2 and y 2 displacements. Fig. 3 gives a schematic description of a mathematical model of a double rotor system which consists of shafts of lengths L 1 and L 2 . The shafts are assumed to be joined by isotropic bellows couplings having linear angular stiffness and damping k r1 , k r2 ; c r1 , and c r2 respectively. The rotors have disks of masses m 1 and m 2 , positioned on each rotor, where b is the length measured from the excitation point where the motor is attached to the geometrical centre of disk 1 and disk 2 is position at the end of rotor 2 at the geometrical center of disk 2. Each disk has a diameter d 1 and d 2 with a thickness v 1 and v 2 . The system is subject to a radial clearance c at the stator in the form of a snubber ring fixed with a linear spring system represented by k s , positioned at length a measured from the geometrical center of disk 2 to the center of the snubber ring. The inertial reference frame xyz is positioned at the origin O, where the motor provides the rotational speed to the system and the subscript i ¼ 1; 2 represents the coordinates of disk 1 and 2.

Mathematical modelling
The acceleration due to gravity acts in the z direction as shown in Fig. 3. The angles w 1 ; h 1 ; w 2 , and h 2 represent the rotational displacements of disk 1 and disk 2. The angle U corresponds to the rotational displacement about the z axis, which can be related to the rotational speed X of the system as U ¼ Xt and _ U ¼ X. The equations of motion were derived using a Lagrangian formulation, as in [63]. For the purpose of the model, the shafts are assumed to be rigid and massless and the dynamics of the stator as an additional degree of freedom are neglected. Contact between rotor 2 and the stator is considered to be intermittent and the times at which such contacts are initiated or terminated have to be determined as part of the solution to the model.

Equations of motion
The entire mechanical system can be represented as where X is the rotor drive speed, and u t ð Þ is a column vector of displacements, defined as System parameters Acceleration due to gravity g ¼ 9:81m/s 2 Diameter of disk 1 Schematic of the mathematical model.
where an upper dot represents the derivative with respect to time t. Moreover, M; G; C, and K represent the mass, gyroscopic, damping, and stiffness matrices, respectively, which are defined in B; see Eqs. (B.1)-(B.10). The excitation vector is represented as F e which contains the excitation due to the mass imbalance of each disk, which is derived from the kinetic energy according to Friswell et al. [64] and F i is the nonlinear force due to the contact between rotor 2 and the stator, defined as : ð6Þ H r ð Þ is the Heaviside step function defined as The nonlinear force due to contact F i does not take into account the friction between the rotor and the stator. The backward whirl is a well known phenomenon that occurs as a result of frictional contact. This work concentrates in analyzing the asynchronous bouncing motion phenomena arising from intermittent impacts in rotating systems. These intermittent impacts have a short duration, thus the friction component can be neglected as shown by Zilli et al. [10] and Shaw et al. [60].
The radial displacement r of the second rotor at distance a is defined as It is convenient to transform Equation (4) from the stationary frame to one that rotates at the drive speed X. Following ; whereũ is the vector in the rotating frame. Differentiating Eq. (11) with respect to time t gives The following transformations can be derived : Applying these transformations to Eq. (4) gives the equation of motion in the rotating frame, expressed as The advantage of working in the stationary frame is to clearly identify the contribution of each parameter in the dynamics of a rotating nonlinear system. The term 2MJX represent the Coriolis acceleration and the term MX 2 represent the centripetal softening. The term T T F e in the rotating frame becomes a constant, and thus the system becomes autonomous. The term T T F i remains the same in the rotating frame because the system is isotropic. The coefficients of the stiffness matrix K, gyroscopic matrix G and mass matrix M were obtained from the parameters of the experimental rig, and are given in Appendix B (see specifically Eqs. B.1, B.5, B.6, and B.10). The matrix coefficients can be estimated by using the experimental parameters of the test rig defined in Table 1.
The equations of motion were solved numerically using the fourth-order Runge-Kutta algorithm ODE45 in Matlab, with event detection to locate changes between contacting and non-contacting motion. The event detection function halted the simulation every time the radial displacement of the contact node crossed the clearance threshold, and then the simulation would be restarted with optimal properties for the new state of contact or non-contact. The mathematical model was developed assuming the intermittent contacts between the rotor and stator are frictionless. Experimental results later show the contacting backward whirl occurred when the rotor is in permanent contact with the stator due to friction. The dry backward whirl is a well known phenomenon and the main objective in this paper was to investigate the novel phenomena of asynchronous bouncing motion.

Natural frequency analysis
A Campbell diagram of natural frequency versus rotor speed is shown in Fig. 4, which was computed using the linearised version of the above model versus the rotor speed. Here, red lines represent the forward whirl (FW), the dashed blue lines represent the backward whirl (BW) of the rotors, and the dashed green line represent the rotor speed. The linear natural frequencies can be estimated by obtaining the eigenvalues of the characteristic equation As is commonly understood, see e.g. [64], the crossings of a rotor's linear natural frequency in a Campbell diagram with the 45 line, corresponding to the rotor drive speed, represent the direct resonance frequencies or critical speeds. From Fig. 4 the critical speeds are obtained at motor speeds of 99.8 rpm (1st BW and FW), 245.80 rpm (2nd BW), and 258.00 rpm (2nd FW). Note that the first forward and backward whirl frequencies are practically identical due to the weak gyroscopic terms.
Alternatively, in Fig. 5 is shown a Campbell diagram in the rotating frame. The common convention used in this paper is that backwards (clockwise) whirling motion has a negative frequency. This sign convention can be used in the rotating frame because the frequencies are signed rotational velocities. In that case, the rotational velocities are obtained from the frequencies in the stationary frame by simply subtracting X. As expected, when X is zero, there is no difference in the FW frequencies between the stationary frame and rotating frame. As X increases, the frequencies in the rotating frame decrease. This is due to the sign convention adopted in this paper. When the FW lines pass through zero at 99.8 rpm and 245.80 rpm, this indicates the first and second critical drive speeds as demonstrated in Fig. 4, with the Campbell diagram in the stationary frame. The usefulness of such a rotating-frame Campbell diagram is that it is now straightforward to identify parameter values at which so-called internal resonances occur. For example, at X % 299:40 rpm and 441:60 rpm respectively the line that corresponds to twice the 1st FW mode crosses the 1st and 2nd BW lines, respectively. Also, at X % 424:00 rpm twice the 1st FW line crosses the 2nd FW line. Each of these crossings implies a potential 2:1 internal resonance between the modes concerned.
The investigations in [10,60] suggest that for rotor speeds just beyond these crossings, one should expect to see asynchronous partial contact limit cycles, and the drivespeeds are referred to as internal resonant drive speeds. Further harmonic lines, for example 3Â or 4Â the FW frequency, have been omitted from Fig. 5 for the sake of clarity. Nevertheless, is noted that these could lead to the identification of further internal resonances at which higher-order n : m bouncing asynchronous limit cycles could occur.

Results
In what follows, experimental results are presented for steady motion, after initial transients have decayed, using the protocol described in 2.3. Later on, two different kinds of bouncing motion solutions are analysed, before presenting equivalent results by numerical simulation.

Experimental bifurcation diagrams
Using the output for different rotor speeds, the motion is quantified, by taking the local maximum magnitude of the radial displacements of rotor 2 during the final time of 55 s in each run. An alternative, stroboscopic Poincaré-map analysis as used in [46] was found to be less helpful because of difficulties in choosing a sensible sampling frequency for the complete range of responses found.
The experimental bifurcation diagram shown in Fig. 6 reveals the bouncing orbits at supercritical drive speeds. The bifurcation diagram in Fig. 6 was obtained by developing a range of tests conducted with different rotor speeds from 121.50 rpm to 636.90 rpm. The black dots represent the local maximum radial displacements of rotor 2 between the time span of 40.0 s to 45.0 s and the dashed blue line represents the stator radial clearance to give an approximate indication of whether the amplitude of the disk is sufficient to cause contact. The results presented in Fig. 6 show a broad spread of points at nearly all drive speeds, indicating that this is a system with low effective damping and highly susceptible to noise. Fig. 6(a) represents the results for the more quiescent initial condition A. No evidence of contacting motion is found for any drive speed in this case, although a higher amplitude of response can be observed around the 2nd critical drive speed. Note that there is no corresponding observable high-amplitude radial displacement near the first critical speed, presumably due to the high pass filter applied to the output. Fig. 6(b) shows similar results for condition B, in which the rotor is initially in contact with the stator. Between 121.50 rpm to 208.90 rpm and from 333.72 rpm to 636.90 rpm the rotor is found to settle to synchronous continuouslyin-contact whirling motion. In contrast, the peaks in Fig. 6(b) from 218.79 rpm to 324.10 rpm show where partial contact solutions occur. Fig. 7 shows the fast Fourier transform (FFT) spectra recorded between 30.0 s and 59.0 s for each rotor speed between 121.5 and 636.9 rpm. Fig. 7(a) shows the results for condition A, which has a clear response at the drive frequency for all drive speeds, which indicates that the motion is largely synchronous. Nevertheless, at low drive speeds (<240 rpm) there appears to be significant sub-synchronous components, indicated by a further diagonal. Furthermore, at drive speeds above approximately 380 rpm, there is also a small but clear signal at the approximate natural frequencies of the rotor, denoted by the near horizontal lines. For example, at approximately 260 rpm, there is a peak that is consistent with the 2nd fundamental resonance speed given in Fig. 4. Finally, for speeds between about 280 rpm and 380 rpm there is more complex behaviour, with more asynchronous frequency components excited. Fig. 7(b) shows similar information for initial condition B, where the drive-speed x ¼ X is less prominent and the spectra are in general more complex, indicating significant asynchronous response. A description of two different regimes in this data is now presented in detail.

Internally resonant bouncing response
The high amplitudes in the data in Fig. 7(b) between rotor speeds of 256 rpm to 343.30 rpm are low frequencies which are close to the first backward and forward whirl frequencies of the system. High peaks at these low frequencies correspond to the occurrence of the partially contacting, bouncing responses.
The initial aim of this work was to obtain an experimental reproduction of asynchronous partial contact orbits for a multidisk system as predicted in [60]. Those bouncing orbits induced by internal resonance were predicted to have the following properties: 1. They are aperiodic in a stationary coordinate frame, but periodic in a coordinate frame that rotates with the shaft 2. They occur due to internal resonance between 2 modes, and hence occur at drive speeds just above the nonlinear critical speed for the form of internal resonance concerned, and cannot occur at drive speeds lower than this. 3. Their orbits in the rotating frame often include a characteristic 'double loop' in the rotating frame. Fig. 8 shows an orbit reconstruction of the response under initial condition B, for X ¼ 304:94 rpm. Fig. 8(a) shows the time series of the radial displacement. After approximately 35.0 s, the response appears to establish an approximately steady, albeit chaotic, response that makes repeated contacts with the stator. Fig. 9 shows the orbit of this limit cycle viewed in the stationary frame, and from this view alone it is hard to discern any order in the motion. However, with the aid of Fig. 8(b), which highlights and labels the final few maxima of the radial displacement, it may be seen that these maxima precess around the stator in a clockwise direction by following the sequence of the labels. Fig. 10(a) gives a clearer perspective by displaying the motion in the rotating frame, where a noisy form of the double loop orbits seen in the rotating frame for 2:1 resonances in [60] may be perceived, but in a quasiperiodic form. (Quasiperiodic responses can be found in internally resonant responses of static structures, for example in a beam with a 3:1 resonance in [65].).
A spectral analysis (see Fig. 5), shows that a drive speed of 304.9 rpm is just above a point of intersection between the second harmonic of the first forward whirl mode and the first harmonic or the first backward whirl mode. This confirms that the drive speed is just beyond that at which a 2:1 internal resonance occurs between these modes. The segment of data and the labels highlighted in Fig. 8(b) are shown in the rotating frame in Fig. 10(b), from which it may be seen that the maxima occur at roughly the same angular position in the rotating frame. This further confirms that the orbit is showing approximately cyclic behaviour in this coordinate frame. Fig. 11 shows the FFT spectra of this orbit in (a) the stationary and (b) the rotating frame, when the motor speed is 304.94 rpm for condition B. In (a), there is a peak at 5.083 Hz, equivalent to a motor speed of 304.94 rpm, and so this is the direct response. Apart from this, a group of 4 peaks is seen between 1.4 Hz and 2.0 Hz, which it should be noted are in the vicinity of the first forward and backward whirl frequencies (see Fig. 4). There is also some noise near zero Hz attributed to the process of integrating the signal from the accelerometer, and some other small peaks.
It is noteworthy that there is relatively little frequency content near the 2nd backward and 2nd forward whirl speeds despite the drive speed being much closer to these frequencies than the 1st backward and forward whirl speeds. As expected for asynchronous vibration, apart from the direct response no peaks are at the frequencies commensurate with the motor speed or any other peak frequency. When considering the rotating frame as shown in Fig. 11(b) it is expected that for an assumed isotropic and axisymmetric system the motor speed becomes effectively zero. However, there is a peak near the motor speed frequency which the authors attribute to the effects such as misalignment or bearing stiffness asymmetry that break these assumptions, noting that fixed features in the stationary frame will appear to rotate at the motor speed in the rotating frame.   Fig. 8(b) with corresponding maxima.
The transformation to the rotating frame can be regarded as simply a shift in frequencies (or rather rotational velocities), so therefore equivalence between peaks in each graph can be drawn by comparing their magnitudes and allows to identify which peaks are due to forwards whirling motions, and which are due to backward whirling. For example, the peak at  Fig. 8(b).  Fig. 11(a) clearly relates to the peak at 6.50 Hz in Fig. 11(b), because they both have similar magnitudes. This peak clearly corresponds to backward whirling motion, because it appears to whirl backwards more quickly in the rotating frame and it may be verified that 1.42 Hz + 5.08 Hz = 6.50 Hz.

Hz in
Similarly, a relatively slow forward whirl in the stationary frame will appear to be whirling backwards at the difference between the whirl speed and the motor speed, and can be verified using the peaks in each graph with magnitudes of approximately 3 Â 10 À3 m and that 5.08 Hz-1.62 Hz = 3.46 Hz. If all the 4 peaks found between 1.4 Hz and 2.0 Hz in (a) are considered in this way, it can be seen that the motion is dominated by two forward whirls with similar frequencies and two backward whirls with similar frequencies. The presence of two nearby frequencies is suggestive of a slowly modulated sine wave, and it may be noted that the forward frequencies and the backward frequencies are both separated by 0.42 Hz. This relates to a modulation period of 2.40 s, and inspection of Fig. 8(a) shows that this matches the approximate period between high amplitude peaks in the settled region. It is interesting to note that the frequencies of 3.46 Hz and 6.92 Hz that appear in Fig. 11(b) are in an exact 2:1 ratio, giving evidence that this motion is fundamentally related to this form of internal resonance in the rotating frame. Finally, it may also be commented that the small peak at 9.5 Hz in Fig. 11(b) is at the approximate 3rd harmonic frequency, and it is to be expected that as a response featuring significant nonlinearity there will be frequency content at all higher harmonics (see for example [13]), although since these do not coincide with resonance these are less pronounced.
The presence of distinct peaks in these FFT graphs allow us to conclude that this response is not fundamentally chaotic in nature, despite the haphazard appearance of Fig. 8. Furthermore, the fact that the irregular spacing of peaks seen in the stationary frame FFT becomes approximately regularly spaced in the rotating frame FFT (apart from the noise at the drive frequency) gives further confirmation that the motion is approximately periodic in the rotating frame. Furthermore, the fact that the amplitude of these groups of peaks becomes much smaller after the first two groups suggests that the first 2 dominate due to resonance, and this provides evidence that a 2:1 resonance is active in this motion.
FFT data for another case where this type of motion was observed is given in Fig. 12(a), although the peaks near the 1st whirling modes are less distinct and there is more motion near the 2nd whirling frequencies giving a spike at 4.375 Hz. It should also be noted that in Fig. 12(a) the peaks at 3.458 Hz and 6.958 Hz are very close to a 2:1 relationship.

Bouncing backward whirl motion
The highest amplitude responses in Fig. 7(b) occur between rotor speeds 400 rpm and 636.90 rpm and were observed to be a predominately backward whirling response. Confirmatory evidence is that the angular speed during roll can be esti- Again, there appears to be a small peak of synchronous response in the region of the 2 nd critical speed. This motion is shown in Figs. 13 and 14. It can be seen that the motion shows a clearer order than that in Fig. 9, and the labeled maxima show that the motion precesses far more rapidly around the stator. However, as can be seen in Fig. 15, the motion is cyclic in the rotating coordinate frame. The stationary frame frequency content shows greater order as well, with Fig. 16 showing that the response is dominated by a signal near the 1st backward whirl frequency. The backward whirling is confirmed by the corresponding peak moving to 7.041 Hz in the rotating frame as shown in Fig. 16(b), which is approximately the sum of the two frequencies shown in Fig. 16(a). This motion seems to be fundamentally different from the previously presented cases because the frequency content suggests that a single mode, the 1st backward whirl, is predominant. Fig. 17(a,b) shows the bifurcation diagram obtained with the simulation for rotor speeds from 5 rpm to 650 rpm. The condition A represented in Fig. 17(a) corresponds to the angular initial conditions of 0 rad. The condition B represented in Fig. 17 (b) corresponds to the initial angular displacement for rotor 1 and rotor 2 in the lateral in plane direction of 0.05 rad, in the lateral out of plane direction of 0 rad, and with 0 rad/s initial angular velocity. In Fig. 17(a,b) the markers represent local maxima of radial displacements of rotor 2 between the time span of 40.0 s to 45.0 s and the dashed blue line represents the stator radial clearance.

Numerical bifurcation diagram
Note from Fig. 17(a) how below 95 rpm and from 290 rpm to 650 rpm the simulations always settle to synchronous whirling motion, where this motion is linear forward (non-contacting) whirl or permanently-in-contact forward whirl. The high amplitude radial displacements shown in Fig. 17(a) from 100 rpm to 285 rpm occur near the fundamental resonance frequencies of the system. The cases shown in Fig. 17(b) below 85 rpm and from 280 rpm to 720 rpm always settle to synchronous non-contacting whirling motion. In Fig. 17(b) the high radial displacements between 90 rpm to 275 rpm are near the first and second backward and forward whirl resonance frequencies.
The separate markers shown in Fig. 17(b) between 90 rpm and 275 rpm for condition B show where the partial contact solution occurs. There is significant qualitative and quantitative agreement between this simulated onset of bouncing motion and the corresponding experimental results in Fig. 6(b).
There is a noticeable difference between the responses shown in Fig. 17(b) and Fig. 6(b) though at much higher drive speeds. In the experiments for initial condition B, (Fig. 6(b)) shows high-amplitude contacting response, whereas in the simulations ( Fig. 17(b)) shows only non-contacting solutions, as in condition A. Further investigation of this high-X experimental response showed that the motion is dominated by contacting backward whirl in which the rotor rolls around the stator, due to friction. Friction was not modelled in the simulation and hence this phenomenon was not observed there. Neverthe- less, the primary objective of this work is to investigate the novel phenomenon of bouncing motion rather than the relatively well-understood phenomenon of dry backward whirl.   Fig. 13(b).

Discussion
This work has presented experimental results concerning the bouncing motion, in a rotor-stator contact coupled double rotor system, designed to elucidate rotor-stator contact phenomena at low frequencies. This rig exhibits many coexisting responses including intermittent impacts, non-resonant, and contacting whirl modes. High amplitude responses were seen at both primary resonances or far away from the primary resonances of the rotating system. In real industrial applications this means that transient effects can push the rotor from non-contacting whirling motion into bouncing motions, and highlights the complexity of handling these responses in models of real machines.
The initial objective of the work was to provide experimental verification of asynchronous partial contact orbits due to internal resonance as predicted by Shaw et al. [60]. This objective has been largely achieved. Predicted features of bouncing responses such as modes excited far away from direct resonance with out of balance excitation, cyclic behaviour when viewed in a rotating coordinate frame and exact 2:1 relationships between some frequency components when viewed in the rotating coordinate frame were shown. However, the results were far from the exactly periodic results seen previously in simulation. Instead, the amplitude of responses modulated relatively slowly giving a quasiperiodic character to the response. Further work needs to identify specific factors that lead to this behavior.
Moreover, the results presented in this paper demonstrate at least excellent qualitative agreement between experimental and simulation results for the central phenomenon under investigation. First and foremost, this work has demonstrated the existence of extreme bi-stability in both the experimental and the simulation results. With benign, non-contacting initial conditions observed only on the forward whirl motion. Encouragingly, there is broad agreement between simulation results, experimental test results and the linear theory in [60] as to the parameter regions, where a 2 : 1 resonant bouncing motion is expected to be seen. Furthermore, the orbits from the experimental testing have similar shape as the one predicted from the mathematical model.
Clearly, the experimental results lack the pristine geometrical elegance of the results found in [60], which claims that asynchronous impacting orbits are quasiperiodic in the stationary frame, while being a periodic internal resonance in the rotating frame. However, the comparison to Fig. 5, the presence of 2:1 frequencies in (b) and the greater order visible in the rotating frame shown by a rotating coordinate system in Fig. 10(a) compared to the stationary frame orbits seen in Fig. 9 all suggest that many of the insights of [60] apply. The primary complication would seem to be that the rotating frame orbits are additionally quasiperiodic in the experimental results. Possible reasons for this could be the well-known torus bifurcation seen in many internally resonant systems, or the presence of symmetry-breaking or other unmodelled features in the experiment. A feature not included in the mathematical model is friction at the contact between the rotor and stator. As the impacts between the rotor and stator are instantaneous, friction was not considered when developing the mathemat- ical model due to the numerical results shown in [60]. Future work should include friction to see how the asynchronous bouncing motion occurs with intermittent contacts.
In addition to the orbits predicted by the 2:1 internally resonant orbits, a further experimental response was seen that is periodic in the rotating frame and involving repeated stator impacts. This appeared to involve participation by a single mode, as opposed to being due to internal resonance. These motions were primarily due to backward contacting whirl motion, which is known to occur as a result of frictional contact. The previous studies have shown how continuously in contact backward whirl motion can become unstable to a lift-off instability [47], which could provide another path to asynchronous bouncing motion besides internal resonance.
An interesting avenue for future experimental work would be to investigate how asynchronous bouncing motion evolves as the rotational speed of the motor is slowly varied. One should expect to see significant hysteresis between spin up and spin down. Other possibilities would be to systematically vary properties of the experimental test rig, such as greasing the snubber to reduce friction, or to adjust the stiffness of the bellows couplings.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

H:
This formula describes a direct linear relationship between the steady-state whirling amplitudeHof the system in Fig. 2 and the eccentricity of the disk. k 2 ¼ k r2 L 1 L 2 : ðB:12Þ k 3 ¼ k r2 þ m 2 gL 2 : ðB:13Þ The excitation vector due to the mass imbalance of each disk F e is defined in (B.14). : ðB:14Þ