Mzaki Dakel

Sébastien Baguet

Régis Dufour
email: regis.dufour@insa-lyon.fr

Non-linear dynamics of a support-excited flexible rotor with hydrodynamic journal bearings

Keywords: Rotordynamics, asymmetric rotor, on-board rotor, finite element method, hydrodynamic journal bearing, support motion, parametric excitation, non-linear dynamics, transient motion, dynamic instability, bifurcation, chaos

The major purpose of this study is to predict the dynamic behavior of an on-board rotor mounted on hydrodynamic journal bearings in the presence of rigid support movements, the target application being turbochargers of vehicles or rotating machines subject to seismic excitation. The proposed onboard rotor model is based on Timoshenko beam finite elements. The dynamic modeling takes into account the geometric asymmetry of shaft and/or rigid disk as well as the six deterministic translations and rotations of the rotor rigid support. Depending on the type of analysis used for the bearing, the fluid film forces computed with the Reynolds equation are linear/non-linear. Thus the application of the Lagrange's equations yields the linear/non-linear equations of motion of the rotating rotor in bending with respect to the moving rigid support which represents a non-inertial frame of reference. These equations are solved using the implicit Newmark time-step integration scheme. Due to the geometric asymmetry of the rotor and to the rotational motions of the support, the equations of motion include time-varying parametric terms which can lead to lateral dynamic instability. The influence of sinusoidal rotational or translational motions of the support, the accuracy of the linear 8-coefficient bearing model and the interest of the non-linear model for a hydrodynamic journal bearing are examined and discussed by means of stability charts, orbits of the rotor, time history responses, fast Fourier transforms, bifurcation diagrams as well as Poincaré maps.

Introduction

Rotating machines are among the indispensable parts of the modern engineering industries. Due to the nature of the excitations applied to the machine rotors, most of them can be considered as on-board rotors affected mainly by both the mass unbalance and the support motions. Generator, pump, compressor and gas turbine rotors installed in power plants as well as vehicle turbochargers are examples of rotors on moving support. The unavoidable mass unbalance is due to the eccentricity of the center of mass along the rotor axis. The rotor balancing aims at reducing the mass unbalance but is not able to cancel it completely (see Xu et al. [START_REF] Xu | The optimization technique-based balancing of flexible rotors without test runs[END_REF] and Kang et al. [START_REF] Kang | Optimal balancing of flexible rotors by minimizing the condition number of influence coefficients[END_REF]). Rotors can also be subject to the frequent external movements of their support which can increase the flexural vibration of the rotors and create an unstable dynamic behavior.

In order to enhance the performances, the safety factor, the air-gap and the disk-stator gap have to be reduced. Therefore the prediction of dynamic behavior of the rotating machine components must be performed more and more carefully. In this context, the literature comprises numerous books studying a wide variety of phenomena related to the dynamics of symmetric/asymmetric rigid/flexible rotor systems mounted on linear/non-linear elastic bearings in the case of a fixed support. Just few of them can be cited here (Lalanne and Ferraris [3], Genta [START_REF] Genta | Dynamics of Rotating Systems[END_REF], Bachschmid et al. [START_REF] Bachschmid | Cracked Rotors[END_REF]). Some works observed the instability of parametrically excited systems. Dufour and Berlioz [START_REF] Dufour | Parametric instability of a beam due to axial excitations and to boundary conditions[END_REF] accomplished the study of Berlioz et al. [START_REF] Berlioz | Dynamic behavior of a drill-string: Experimental investigation of lateral instabilities[END_REF] on the dynamics of a drill-string. They analyzed the time-varying parametric equations of motion of the system by employing the Rayleigh-Ritz method, the Floquet theory and the Friedmann approach. The design using a computer opened doors for a technique called "finite element method" to be applied to the rotor modeling and numerous studies contributed to this topic. The real breakthrough in the finite element modeling of a rotor was presented by Nelson and McVaugh [8] who studied a model called "Rayleigh beam" and taking into account the rotary inertia of the shaft and the gyroscopic effects dependent on the speed of rotation of the rotor. Kang et al. [START_REF] Kang | Investigation on the steady-state responses of asymmetric rotors[END_REF] modeled rotor-bearing systems using Timoshenko beam finite elements and by taking into account the asymmetry of disk, shaft and/or bearing. They showed that the resonant speeds could change due to various angles between major axes of disk and shaft, the shaft asymmetry as well as the bearing characteristics. Nandi and Neogy [START_REF] Nandi | An efficient scheme for stability analysis of finite element asymmetric rotor models in a rotating frame[END_REF] proposed an efficient analysis of stability for finite element models of asymmetric rotors and investigated whether an unstable rotor could be stabilized using an isotropic viscous damper. Theoretical and experimental investigations for the isolated and internal resonances of non-linear forced and parametric oscillations of an asymmetric rotor with non-linear spring characteristics were carried out by Ishida et al. [START_REF] Ishida | Vibrations of an asymmetrical shaft with gravity and nonlinear spring sharacteristics (isolated resonances and internal resonances)[END_REF].

Hydrodynamic bearings play a key role in the design of rotor systems and an accurate prediction of the vibration characteristics of the rotor-bearing systems must be made. Therefore the hydrodynamic bearing performance is predicted and reported in the extant literature by mathematical formulations firmly established. The fluid film forces are strongly non-linear functions of the displacement and velocity of the rotor. Although they act locally on the system, rotordynamics is significantly affected and wholly non-linear. Khonsari and Chang [START_REF] Khonsari | Stability boundary of non-linear orbits within clearance circle of journal bearings[END_REF] showed that there was a stable region in the short bearing clearance circle outside of which, any initial conditions for the non-linear transient dynamic problem would yield an unstable orbit even if the hydrodynamic bearing should be stable according to the linearized stability analysis. Zhao et al. [START_REF] Zhao | Stability and response analysis of symmetrical single-disk flexible rotor-bearing system[END_REF] predicted the linearized stability and presented a comparison between linear and non-linear mass unbalance responses of a flexible rotor-hydrodynamic bearing system modeled by the lump mass method to reduce the degrees of freedom. Li and Xu [START_REF] Li | A method to determine the periodic solution of the non-linear dynamics system[END_REF] studied a Jeffcott rotor supported on oil film infinite-length bearings to obtain periodic orbits, their periods and their stability using the generalized shooting method. Harmonic, sub-harmonic, quasiperiodic and chaotic responses and their Poincaré maps were investigated for a rigid rotor on hydrodynamic short bearings by Brown et al. [START_REF] Brown | Chaos in the unbalance response of journal bearings[END_REF], for a finite element rotor model on hydrodynamic elliptical bearings by Zheng and Hasebe [START_REF] Zheng | An efficient analysis of high-order dynamical system with local nonlinearity[END_REF] and for a rigid rotor on hydrodynamic elliptical bearings by Shen et al. [START_REF] Shen | Nonlinear behavior analysis of a rotor supported on fluid-film bearings[END_REF]. Kishor and Gupta [START_REF] Kishor | On the dynamic analysis of a rigid rotor-gear pair-hydrodynamic bearing system[END_REF] used non-linear analytical characteristics for hydrodynamic long bearings and introduced them in a rigid rotor-bearing model with a simplified spur gear model. Baguet and Velex [START_REF] Baguet | Influence of the nonlinear dynamic behavior of journal bearings on gear-bearing assemblies[END_REF] combined a finite element shaft model based on the Timoshenko beam theory with both gear and bearing non-linearities which were represented by spur gear mesh stiffness and hydrodynamic short bearing forces. Baguet and Jacquenot [START_REF] Baguet | Nonlinear couplings in a gear-shaft-bearing system[END_REF] extended the study presented in [START_REF] Baguet | Influence of the nonlinear dynamic behavior of journal bearings on gear-bearing assemblies[END_REF] to include helical gears as well as hydrodynamic finite-length bearings and analyzed parametrically the shaft orbits, the dynamic tooth loading factor as well as the hydrodynamic forces.

Some works focused on the dynamic behavior of a rotor under seismic, random or shock excitations of its support. Srinivasan and Soni [START_REF] Srinivasan | Seismic analysis of a rotor-bearing system[END_REF] studied the effect of spin, support rotation as well as axial force and axial torque on the seismic response of a rotor-bearing system. Samali et al. [START_REF] Samali | Random vibration of rotating machines under earthquake excitations[END_REF] used the Monte Carlo simulation to simulate the non-stationary earthquake ground motions and to determine the statistics of rotating machinery response. Hori and Kato [START_REF] Hori | Earthquake-induced instability of a rotor supported by oil film bearings[END_REF] examined a seismic response of a Jeffcott rotor mounted on oil film bearings to a real seismic wave and investigated its stability. Suarez et al. [START_REF] Suarez | Seismic response of rotating machines[END_REF] observed that even for strong rotational inputs, the parametric terms in the equations of motion of a finite element rotor model subject to six components of support excitations could be neglected without affecting its seismic response. Subbiah et al. [START_REF] Subbiah | Response of rotors subjected to random support excitations[END_REF] obtained the amplitude power spectral density due to random excitations of the support for studying the rotor response using modal analysis methods. Lee et al. [START_REF] Lee | A finite element transient response analysis method of a rotor-bearing system to base shock excitations using the state-space Newmark scheme and comparisons with experiments[END_REF] proposed a finite element rotor model based on the Timoshenko theory by considering a shock excitation of the support and the statespace Newmark method and focused on the experimental behavior of a rotor under this excitation.

The effect of the flexible foundation on rotordynamics was investigated from theoretical and experimental points of view, see, for example, Bonello and Brennan [START_REF] Bonello | Modelling the dynamic behaviour of a supercritical rotor on a flexible foundation using the mechanical impedance technique[END_REF]. Da Silva Tuckmantel et al. [START_REF] Da Silva Tuckmantel | An analysis on the supporting structure representative model in rotating systems[END_REF] represented the supporting structure (foundation) of a rotating system by coupled as well as uncoupled modes and tested the methods of mechanical impedance and mixed coordinates to calculate the system responses. The experimental tests of Feng and Hahn [START_REF] Feng | Rotor-model-based identification of foundations in rotating machinery using modal parameters[END_REF] showed that even with input data truncated to two significant digits, satisfactory identification were possible for a flexibly supported undamped rigid block foundation in rotating machinery.

Few works are referenced in the literature on the investigation of dynamic behavior of a rotor in the case of a harmonic motion of its support. Duchemin et al. [START_REF] Duchemin | Dynamic behavior and stability of a rotor under base excitation[END_REF] observed the motion stability of a rotor under a sinusoidal rotation of the support by employing the Rayleigh-Ritz and multiple scales methods. They presented also experimental results to validate the analytical study. The work of Driot et al. [START_REF] Driot | Theoretical and experimental analysis of a base-excited rotor[END_REF] was based on the model presented in [START_REF] Duchemin | Dynamic behavior and stability of a rotor under base excitation[END_REF], described the numerical orbits of the rotor and compared them with experimental ones. El-Saeidy and Sticher [START_REF] El-Saeidy | Dynamics of a rigid rotor linear/nonlinear bearings system subject to rotating unbalance and base excitations[END_REF] derived the equations of motion of a rigid rotor linear/non-linear bearing system subject to rotating mass unbalance plus harmonic excitations of the support along or around lateral directions. They presented analytical frequency responses in the case of linear bearings, while they discussed numerical results with regards to the time domain, the fast Fourier transform as well as the Poincaré map in the case of a bearing cubic non-linearity. Das et al. [START_REF] Das | Active vibration control of unbalanced flexible rotor-shaft systems parametrically excited due to base motion[END_REF] investigated the active vibration control of a flexible rotor system modeled by Rayleigh beam finite elements and excited by a mass unbalance as well as a periodic rotational motion of the support. They performed the control with electromagnetic control force provided by an actuator placed on the stator in a suitable plane around the rotor and it was successful for avoiding the lateral parametric instability due to the sinusoidal rotation of the support.

A significant number of previously stated papers treat support-excited rotor systems. Nevertheless, few studies focus on the harmonically excited on-board rotors. Moreover, these papers focus on constructing simplified models. The adoption of simplified models facilitates the understanding of the behavior of such dynamic systems but it is limited in the practical applications and reduces the required accuracy between the theoretical rotor models and the realistic ones under the operating conditions. The most common simplifications and limitations are linked either with the geometry of the rotor (for example, Jeffcott rotor [START_REF] Hori | Earthquake-induced instability of a rotor supported by oil film bearings[END_REF]), with the utilization of the Rayleigh-Ritz method instead of the finite element method for the spatial discretization of the rotor model (for example, [START_REF] Duchemin | Dynamic behavior and stability of a rotor under base excitation[END_REF][START_REF] Driot | Theoretical and experimental analysis of a base-excited rotor[END_REF]), with the rotors supported by elastic bearings with constant damping and stiffness coefficients (for example, [START_REF] El-Saeidy | Dynamics of a rigid rotor linear/nonlinear bearings system subject to rotating unbalance and base excitations[END_REF][START_REF] Das | Active vibration control of unbalanced flexible rotor-shaft systems parametrically excited due to base motion[END_REF]) or with the kind of rotor support motion (simple rotation or translation).

This paper overcomes the previous limitations and contributes to the construction of an improved on-board rotor model. Namely, a rotor with geometric asymmetry due to the shaft and/or the rigid disk is discretized using the finite element method based on the Timoshenko beam theory, mounted on hydrodynamic journal bearings of external forces calculated using the Reynolds equation [START_REF] Frêne | Hydrodynamic Lubrication Bearings and Thrust Bearings[END_REF], and excited by combined deterministic motions (three rotations and three translations) of the support. The vibration analysis focuses on the case of bending of the flexible rotor. Eulerian angles are thus employed to describe the orientation of disk and shaft with respect to the moving support. The kinetic and strain energies in addition to the virtual work of the flexible rotor components are computed. Using the short bearing theory in order to obtain the hydrodynamic fluid forces from the Reynolds equation, two models are considered for the hydrodynamic journal bearing: a simplified linearized model with eight constant damping and spring coefficients and a non-linear model where the bearing is treated as external non-linear forces acting on the shaft. The Lagrange's equations are applied to establish the linear/non-linear differential equations of vibratory motion of the rotor in bending with respect to the rigid support representing a non-inertial reference frame. The governing equations of motion display periodic parametric coefficients due to the asymmetry of the rotor and time-varying parametric coefficients due to the support rotations. These parametric coefficients are considered as generators of internal excitation and can create lateral dynamic instability. The equations of motion are solved by the implicit Newmark time-step integration algorithm which must be combined with the Newton-Raphson iterative procedure in the case of a non-linear bearing model. The influence of sinusoidal motions of the support on the dynamics of the onboard rotor-bearing system is analyzed by means of stability charts, orbits of the rotor, time history responses, fast Fourier transforms (FFTs), bifurcation diagrams and Poincaré maps. In addition, a comparison between the linear and non-linear systems is performed to assess the validity range of the linearized 8-coefficient bearing model.

Problem statement and basic considerations

Let Fig. 1(a) present the basic components of the rotor: the disk (1), the shaft (2), the mass unbalance (3), the bearing (4) and the support [START_REF] Bachschmid | Cracked Rotors[END_REF]. The mass unbalance is distributed continuously as well as in a generic way and cannot completely be canceled. The excitations taken into account are due to the inevitable mass unbalance and to the imposed movements of the support. The transverse deflections of the elastic line of the rotor shaft are studied.

The following assumptions are used: the disk is rigid, the shaft is deformable and modeled by Timoshenko beam elements of constant cross-section solicited in bending along two orthogonal directions, the shaft and/or rigid disk can be asymmetric, the mass unbalance is modeled by concentrated masses, the bearings supporting the shaft are flexible of hydrodynamic short type, the support is rigid and mobile, the rotor rotates at a constant speed Ω .

The consideration of the support motion may modify the form of the equations of motion of a rotor in bending compared to those obtained in the case of a fixed support. In order to develop a modeling as simple as possible, the approach presented by Duchemin et al. [START_REF] Duchemin | Dynamic behavior and stability of a rotor under base excitation[END_REF] is adopted. Three frames of reference are utilized to take into account the movement of the rotor with respect to the support and that of the support with respect to the ground.

The derivation of the equations of motion is inspired from the theory described by Lalanne and Ferraris [START_REF] Lalanne | Rotordynamics Prediction in Engineering[END_REF]. First, the coordinate systems needed to describe the motion of an on-board rotor are defined and the vectors expressing the rotations between them are calculated. Then, the kinetic energies for a disk, a shaft element and a concentrated mass unbalance, the strain energy for a shaft element as well as the virtual work for a bearing are calculated. Moreover, the rotating flexible rotor is discretized using the Timoshenko beam finite elements. Finally, the differential equations of motion of an on-board rotor in bending are deduced from the Lagrange's equations being applied with respect to the generalized coordinates i q by

d d i q i i i T T U F t q q q   ∂ ∂ ∂ - + =   ∂ ∂ ∂   ɺ (1) with ( ) dof dof 1 n i n ≤ ≤
the number of degrees of freedom, T and U the kinetic and strain energies respectively, i q F the external forces corresponding to i q . The symbol " i " refers to the differentiation with respect to time t .

(a) (b)

Kinematics of the on-board rotor

Three principal Cartesian frames of reference shown in Fig. 1(a) are introduced to take into consideration the movement of the rotor support:

( )

g g g g
, , R x y z is the Galilean frame of reference fixed to the ground, ( ) In the classical rotordynamics theory, the transformation between the frames R and l R shown in Fig. 1 R with respect to its support R . The angular velocity vector of the rotor l R with respect to the ground g R measured in the frame l R is described by

, , R
g g l l l l l l T , , R R R x y z R R R R ω ω ω   = + =   ω ω ω (2)
where the superscript T is the matrix transposition symbol. The components ( l x ω , l y ω , l z ω ) are formulated as a function of (ψ , θ , φ ) and their time derivative as well as ( x ω , y ω , z ω ). Further details about the kinematics of the on-board rotor can be found in [START_REF] Dakel | Steady-state dynamic behavior of an on-board rotor under combined base motions[END_REF]. Let us consider an arbitrary point init C along the elastic line. Namely, it represents the geometric center of the non-deformed shaft. Its coordinates in the frame R are ( ) 0, ,0 y . Let ( ) , u y t and ( ) , w y t be its time-dependent dynamic lateral displacements due to bending along the Ox and Oz axes, see Fig. 1(b). As its coordinate along the axial axis is assumed to be a constant with respect to the rotor support, the axial displacement along the Oy axis is not taken into account.

Energy and virtual work calculations of the on-board rotor-bearing system components

This section includes formulations for the disk, the shaft, the rotating mass unbalance and the hydrodynamic journal bearing. The kinetic and strain energies are measured by an observer situated at the ground and their terms are written with respect to the frame linked with the rotor rigid support R . The rotary inertia, the gyroscopic inertia, the shear deformation of shaft as well as the geometric asymmetry of shaft and/or rigid disk are taken into account.

Disk

According to the assumptions presented in Section 2, the disk is assumed to be rigid and only its kinetic energy d T is calculated. This energy consists of the scalar sum of the absolute translational kinetic energy and the absolute rotational kinetic energy about the disk mass center located at the generic abscissa d y along the Oy axis, i.e., ( ) ( )

g g g g l l l l d d d d d d d d d d d d T T mo di mo di d d 1 with diag diag 2 2 R R R R x y z y m m m m m m m m m m O O R R m T I I I I I I I I     = + = = + -     v v ω I ω I (3)
where d m is the mass of the disk,

T T T d d d d d d , , , , , ,
O O O O O O R R R
x y z u y w u x y y w z

= + = + = + + + g l g l d d O O O O OO (4)
Using the transport theorem relating the time derivative of a vector measured from a fixed frame to that of the same vector but measured from a moving frame, the derivation of the position vector ω for the disk, the expression of its kinetic energy becomes

g l d O O yields g g g l l l l d d d d T d d , , d d R R R R R O O O O R u v w t t   = + ∧ = =   g l g l g l d d d O O O O ω O O v ɺ ɺ ɺ
( ) ( ) ( ) ( ) l l l l l l l l d d d d d d 2 2 2 mo 2 2 2 di 2 2 d d d d d d d 1 2 2 x z y y x z m m m O O O m T u v w I I I ω ω ω ω ω = + + + + + + - ɺ ɺ ɺ (6)
The rotor runs at a constant speed of rotation Ω . So the spinning angle φ and its time derivative

θ ψ θ θ θ ψ ψ ψ - - ≃ ≃ ≃ ≃ (7)

Shaft

The shaft is assumed to be flexible and is modeled by beam elements of constant cross-section and constant moments of inertia. The kinetic energy of the shaft element is obtained by taking a shaft elementary volume similar to a disk of very small thickness dy . ω , the expression for the kinetic energy sh T of the shaft element, whose end-points have the abscissas 1 y and 2 y along the Oy axis, is written in the following form

Considering the components

( ) ( ) ( ) ( ) l l l l l 2 2 2 2 l l l sh sh sh 1 1 1 1 2 2 2 mo 2 2 mo 2 di 2 2 sh sh sh sh sh sh 1 d d 2 d d 2 2 y y y y x z y x z S S S O O O y y y y S T u v w y I y I y I y ρ ρ ω ω ρ ω ρ ω ω = + + + + + + - ∫ ∫ ∫ ∫ ɺ ɺ ɺ (8)
where sh ρ , sh S and sh l ( sh 2 1 l y y = -) are respectively the density, the cross-sectional area and the length of the shaft element.

( ) The rigid support motion relative to the ground has no influence on the strain energy of the shaft because the latter depends only on the stresses and therefore on the transverse deflection of the shaft with respect to the rotor support R . In addition to the bending deformation, the shear effects introduced by Timoshenko and the second-order (non-linear) terms of the strain field are accounted. The non-linear strains are linked with the centrifugal stiffening being proportional to the square of the support angular velocity around the transverse axes and corresponding to the geometric potential energy. Thus the strain energy sh U of the shaft element is expressed by

sh
κ ψ θ ψ θ ψ θ κ ψ θ                    ∂ ∂ ∂ ∂ ∂ ∂        = + + + + - - -                    ∂ ∂ ∂ ∂ ∂ ∂                        ∂ ∂ - + - -     ∂ ∂     ∫ ∫ ∫ ( ) ( ) ( ) 2 2 sh 1 1 2 2 1 1 di sh 2 2 2 di
ψ θ Ω ρ κ ψ θ Ω ω ω     ∂ ∂    -    ∂ ∂               ∂ ∂ ∂ ∂   - + - + + - + +           ∂ ∂ ∂ ∂           ∫ ∫ ∫ ∫ (9)
where sh E , ( )

sh sh sh 2 1 G E ν = +
and sh ν are respectively the Young's modulus, the shear modulus and the Poisson's ratio of the isotropic shaft material.

( )

mo sh sh sh 2 x z κ κ κ = + and ( ) di sh sh sh 2 x z κ κ κ = -
are respectively the mean shear correction factor of the crosssectional area and that relative to the section asymmetry of the shaft.

Mass unbalance

The mass unbalance is modeled by a concentrated mass mu m placed at a point mu P with a distance mu r from the geometric center of the cross-sectional area of the shaft. Its initial angle with the Oz axis of the frame R at rest is mu η . The mass unbalance remains in a plane perpendicular to the Oy axis with a constant abscissa mu y along the Oy axis relative to the rotor support R . Its kinetic energy mu T is defined as follows ( )

g g mu mu T mu mu 2 R R P P m T = v v (10)
The translational velocity vector

O O O R u x r t y y w z r t Ω η Ω η = + + = + + + + + + +     g g l l mu mu mu mu O P O O OO O P (11)
Applying the transport theorem, the time derivative of the position vector g mu O P is of the following form

g g g mu d d d d R R R R R P t t = + ∧ = g g g mu mu mu O P O P ω O P v (12)
The external force vectors obtained by the application of the Lagrange's equations to the kinetic energy of the mass unbalance contain the so-called "centripetal and centrifugal force vectors", sinusoidal tangential force vectors and sinusoidal force vectors fixed in space with respect to the support R . Thus in the presence of sinusoidal rotations of the support, the excitation due to the mass unbalance combined with the support rotations is sinusoidal and has not only a frequency equal to the speed of rotation Ω of the rotor but also a combination of frequencies of speed of rotation of the rotor and support rotation. O are respectively the bearing center and the shaft geometric center. The radius, length and clearance of the bearing are respectively be r , be l and be be sh OC . In the dynamic regime, the fluid film pressure distribution in the hydrodynamic bearing is governed by the Reynolds equation and its integration over the bearing yields the external hydrodynamic forces T be be be ,

Hydrodynamic journal bearing

u w R F F   =   F
acting on the shaft and projected in the Cartesian frame connected to the rotor support R . Denoting the dynamic variation of shaft center position in the bearing (i.e., the shaft center describes an orbit situated within a circle whose radius is equal to the radial clearance of the bearing) due to the dynamic variation of forces applied to the rotor (mass unbalance, fluid film forces and support excitations), the dynamic transverse displacement vector of the shaft center is

[ ] T be be be , R u w = δ
. The dynamic polar coordinates of the shaft center in the bearing expressed in the frame ( ) be , R ta ra are the relative dynamic eccentricity be ε (relative dynamic radial displacement) and the dynamic attitude angle be ϕ , see Fig. 2(a). These coordinates and their time derivatives are expressed as 2 2 be be be be be be be be be be be be be be be be 2 2 2 be be be be be be be

; tan ; ; u w e u u u w w u w w u c c w c c ε ϕ ε ϕ ε ε + + - = = = = = ɺ ɺ ɺ ɺ ɺ ɺ (13)
In the present paper, the short bearing theory is considered (i.e., be be 1 8

l d ≤
where be be 2 d r =

). Assuming an isothermal and laminar flow of an isoviscous incompressible fluid, the dynamic polar (tangential and radial) components of the hydrodynamic force vector can be analytically deduced from the integration of the hydrodynamic fluid pressure field based on the Reynolds equation using the Gümbel boundary conditions which discard the negative pressure from the force integral. These dynamic components in the polar frame linked with the hydrodynamic bearing be R are defined by the following formulations [START_REF] Frêne | Hydrodynamic Lubrication Bearings and Thrust Bearings[END_REF] ( )

( ) ( ) ( ) ( ) 2 3
3 be be 2 2 be be be be be be be be be be be be be be 2 2 2 2 2 2 2 be be be be be

1 2 4 1 2 ; 2 2 2 1 2 1 2 1 ta ra r l r l F F c c πε ε µ ε µ π ε ε Ω ϕ ε Ω ϕ ε ε ε   +     = + - - = - + -       - - -   ɺ ɺ ɺ ɺ (14)
where be µ is the fluid film dynamic viscosity. The equation calculating the static solution can be reached by setting the time derivatives of be ε and be ϕ in Eq. ( 14) to zero and by setting the magnitude of the vector sum for the force vector components in Eq. ( 14) to be equal to the static load r W . Thus the resulting non-linear equation is solved by an iterative Newton-Raphson method and provides the relative static eccentricity 0 be ε and then the static radial displacement 0 be e of the shaft center in the bearing. The components of the static equilibrium position vector 0 be δ and those of the dynamic fluid force vector be F are obtained by a classical change of basis (see Fig. 2(a)) given as ( ) ( ) be be be be be be 0 0 T T be be 0 0 0 be be be be be 0 0 0 be be be be be be

0 cos sin cos sin ; ; ; sin cos sin cos ta R R R R ra R R R R R R F e F ϕ ϕ ϕ ϕ ϕ ϕ ϕ ϕ -     -       = - = - = =                   δ R F R R R (15)
where be 0 R R R and be R R R are the rotation matrices of the polar frame be R with respect to R in the static and dynamic cases respectively. According to Eqs. ( 13)-( 15), since the hydrodynamic force vector components depend non-linearly on the dynamic polar coordinates and on their time derivatives, they are non-linear functions of the dynamic transverse displacement and velocity vectors ( ) ) ( ) 0 be be be be be be be be be

, , = - - F δ δ F δ 0 c ∆δ k ∆δ ɺ ɺ (17)
with ( ) ( ) 0 0 be be 0 be be be be be be be be be be be be be be be be be be be , ,

;

; ;

xx xz xx xz zx zz zx zz c c k k c c k k     ∂ ∂ = - = = - = = - =     ∂ ∂     δ 0 δ 0 F F c k ∆δ δ δ ∆δ δ δ δ ɺ ɺ ɺ (18)
The bearing force vector Their analytical expressions can be found in [START_REF] Frêne | Hydrodynamic Lubrication Bearings and Thrust Bearings[END_REF] and deduced from Eq. ( 14). In addition, these expressions are presented in Appendix A. Lastly, it is necessary to execute the transformation between the frames be R and R in the static condition to produce the matrices be c and be k in R as follows ( ) ( ) be be be be be be T T 0 0 0 0 be be be be

; R R R R R R R R R R = = c R c R k R k R (19)

Overall dynamic behavior of the on-board rotor-hydrodynamic bearing system

Finite element formulation

The finite element method is chosen to discretize the on-board rotor in bending. Each node of nodal interpolation includes four degrees of freedom (two transverse translations and two rotations). The disk, the mass unbalance and the hydrodynamic journal bearing are modeled by one-node finite elements situated at the geometric center of the cross-sectional area of the shaft. On the other hand, the shaft is modeled by esh n two-node beam finite elements based on the Timoshenko beam theory and thus has esh 1 n + nodes.

Let us consider a non-deformed shaft finite element i esh which consists of two nodes i n and 1 i n + of coordinates sh where ( )

esh 1 i i n ≤ ≤
is the shaft finite element number. Its nodal displacement vector n i esh δ projected in the coordinate system attached to the rotor rigid support R is defined by

1 1 1 1 T n , , , , , , ,
i i i i i i i i i i i i i i i i i n n n n n n n n esh esh esh esh esh esh esh esh esh R u w u w θ ψ θ ψ + + + +   =   δ (20)
The above vector is separated into two vectors n, xOy and yOz , then

1 1 1 1 T T n, n, , , , ; , , ,
i i i i i i i i i i i i i i i i i i n n n n n n n n u w esh esh esh esh esh esh esh esh esh esh R R u u w w ψ ψ θ θ + + + +     = =     δ δ (21)
The dynamic transverse displacements ) along the Oy axis of the frame R and include the shear effects (see Chen and Ku [START_REF] Chen | Dynamic stability of a cantilever shaft-disk system[END_REF]).

ψ θ ψ θ = = = = F δ F δ F δ F δ (22) where ( )

Equations of motion

The Lagrange's equations are applied to the energies for the disk, the shaft finite element and the mass unbalance (see Eqs. ( 6), ( 8), ( 9) and ( 10)) as well as to the virtual work of the hydrodynamic bearing (see Eq. ( 16)) and the resulting elementary matrices and vectors are appropriately assembled. Due to the presence of hydrodynamic bearings, the boundary conditions do not exist and the corresponding row and column entries of the assembled matrices and vectors are not eliminated. Thus the obtained global matrices and vectors are respectively of dimension g g n n ×

and g 1 n × where ( )

g esh 4 +1 n n =
is the total number of degrees of freedom. The obtained non-linear second-order differential equations describing the dynamic behavior of the on-board rotor-bearing system in bending are written with respect to the non-inertial frame connected to the rigid support R as follows

( ) ( ) ( ) ( ) ( ) r r r r r r r be r r , t t t t + + = + M δ C δ K δ F F δ δ ɺɺ ɺ ɺ (23)
where ( )

r t M , ( ) r t C and ( ) r t K
are the global mass, damping and stiffness matrices with periodic and time-varying parametric terms due to the geometric asymmetry of the rotating rotor and to the rotations of its rigid support respectively. These matrices are not affected by the translational motions of the support. r δ ɺɺ , r δ ɺ and r δ are the global acceleration, velocity and displacement vectors conforming to the connectivity of the finite elements. ( )

r t F
is the global external force vector containing the excitations due to the influence of the mass unbalance as well as to that of the rotational and translational motions of the support. Lastly, the external nonlinear hydrodynamic force vector

Ω Ω = + + M M M M (24)
( )

+ + + + + K K K K K K K K K K K K K ɺ ɺ ɺ ɺ ( ) 2
+ + + + + + + + + K K K K K K K K K K ɺ ɺ ( ) 2 re, ,s sh,su sin 2 x z x z t ω ω ω ω Ω (26)
The subscripts " d ", " sh " and " su " refer to the disk, the shaft as well as the support respectively and express the contribution to the phenomena represented by the corresponding matrix. The superscripts " 2 c " and " 2 s " denote the geometric asymmetry of the rotor expressed in terms of the time-varying trigonometric functions cos 2 t Ω and sin 2 t Ω . The superscript " g " stands for the rotor gyroscopic effect, " e " for the shaft elasticity corresponding to the bending and shear deformations, " re " for the rotational effects due to the support rotations (these effects come from the kinetic energies of the disk and the shaft) and " gse " for the geometric stiffening effects corresponding to the centrifugal stress due to the support rotations (these effects come from the strain energy of the shaft).

The force vector ( ) r t F indicated in Eq. ( 23) is expressed as follows

( ) ( ) ( ) ( ) ( ) ( ) r 2
Ω Ω Ω Ω = - = + F V F V V (28)
( ) ( ) (

2 2 2 1 1 1 1 1 1 2 2 2 1 1 1 1 1 ,c ,c ,c ,c ,c ,
+ + + + + F V V V V V V V V V V V V ɺ ɺ ɺ ɺ ) 1 ,s sin z x z t ω ω Ω (29) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) (
θ ψ ω ω ω ω ω ω ω ω ω ω ω ω ω ω ω ω ω ω ω ω ω ω ω ω ω ω ω ω ω ω = - + - + + - - - + - + - + + - - - + - + + - - + - - F V V V V V V ɺ ɺ ɺɺ ɺ ɺ ɺ ɺ ɺ ɺ ɺɺ ɺ ɺ ɺ ɺ ( ) ( ) ( ) d,sh,su d,sh,su x y y x x y y z y z ψ θ ω ω Ωω ω ω Ωω ω ω - + + + V V (30) ( ) ( ) ( ) 2 2 2 c ,c ,c d,sh,su d,sh,su d,sh,su 2 2 x z y z z x x y t θ ψ ω Ωω ω ω ω Ωω ω ω = - - - + + + F V V ɺ ɺ (31) ( ) ( ) ( ) 2
= - - + + + F V V ɺ ɺ (32)
where the load vectors mu V , d,sh,su

V

and mu,su

V

and their corresponding force vectors mu F , d,sh,su F and mu,su F are associated respectively with the mass unbalance, the inertia force due to support motions and that due to coupling between both phenomena. The superscripts " 1 c " and " 1 s " signify the components of the mass unbalance force expressed in terms of the time-varying trigonometric functions cos t Ω and sin t Ω . The superscripts u , w , ψ and θ denote the direction of the action force components associated with the rotor support motions.

For the sake of conciseness, the elementary matrices involved in ( )

r t M , ( ) r t C , ( ) r t K
and the elementary vectors involved in ( ) r t F as well as the validation of the finite element calculations against the Rayleigh-Ritz formulation are not detailed here and can be seen in [START_REF] Dakel | Steady-state dynamic behavior of an on-board rotor under combined base motions[END_REF].

As already stated above, the equations of motion are non-linear. This is due to the hydrodynamic bearings (local components) even if they do not concern all the degrees of freedom of the rotor system. In the two following sections (i.e., Sections 5.3 and 5.4), the treatment of these equations is detailed.

Non-linear dynamic analysis

Classical transient analysis

For large hydrodynamic bearing forces and thereby large displacements of the rotor, the linear dynamic analysis is no longer valid. In this case, the non-linear analysis of a bearing is considered. Solving the complete non-linear dynamic problem requires the simultaneous solution of the equations of motion (i.e., Eq. ( 23)) of the rotor-bearing system and Eqs. ( 13)-( 15) calculating the nonlinear bearing forces obtained by the Reynolds equation in the dynamic regime. As a consequence, the analysis uses a complex solving algorithm combining the Newton-Raphson incremental-iterative procedure with the Newmark time integration scheme. The semi-discrete equations of motion (Eq. ( 23)) applied at time

( ) ( ) ( ) ( ) ( ) ( ) 1 1 1 1 1 1 r r r 1 r r 1 r r 1 r r 1 be r r , i i i i i i i i i i t t t t + + + + + + + + + + = + + - - = G δ M δ C δ K δ F F δ δ 0 ɺɺ ɺ ɺ (33)
The velocity and displacement vectors of the rotor nodes between i t and 1 i i t t t ∆ + = + , t ∆ being the time step, are approximated by the following general formulations of the implicit Newmark scheme

( ) ( ) ( ) ( ) ( ) 2 1 1 1 1 r r 1 r 1 r r r r 2 r 2 r 1 ; 2 1 2 2 i i i i i i i i i t t t ∆ ∆ ξ ξ ∆ ξ ξ + + + + = + + - = + + + - δ δ δ δ δ δ δ δ δ ɺ ɺ ɺɺ ɺɺ ɺ ɺɺ ɺɺ (34)
where the parameters 1 1 2 ξ = and 2 1 4 ξ = define the constant average acceleration scheme and assume an unconditional numerical stability and second-order accuracy. Since Eq. ( 33) is non-linear, an iterative strategy is used to treat the non-linearity. Namely, the consistent linearization of Eq. ( 33) is applied by finding a first-order Taylor series of this equation and the Newton-Raphson iterative algorithm is taken into account as follows

( ) ( ) ( ) ( ) ( ) 1 r 1 1 1 NR 1 NR 1 r r r r r r r r r r r 1 r d d k i k i k i k i k i i + + + + + + + + = ⇒ = - δ G G δ G δ ∆δ 0 J δ ∆δ G δ δ ≃ (35) where NR 1 1 1 r r r k i k i + + + = - ∆δ δ
δ is the iterative increment of the displacement vector and k is the Newton-Raphson iteration. The Jacobian matrix r J is a function of the displacement vector 1 r i+ δ at iteration k and is given by

( ) ( ) ( ) ( ) ( ) ( ) ( ) 1 1 1 r r r 1 1 1 1 1 r r r r r 1 r r r 1 r 1 r 1 2 1 1 1 1 1 1 r r r r r
= = + + - = = δ δ δ F G δ δ δ δ J δ M C K I I δ δ δ δ δ δ ɺ ɺɺ ɺ ɺɺ ɺ ( 36
)
where I is the identity matrix. The derivatives of the acceleration and velocity vectors 1 with respect to the displacement vector 1 r i+ δ are deduced from the approximations of the Newmark scheme (i.e., Eq.( 34)). The derivative of the non-linear hydrodynamic forces in the Jacobian matrix of Eq. ( 36) depends on both the displacement and velocity vectors ( )

1 1 r r , k i k i + + δ δ ɺ
, is updated at each iteration k and is expressed in what follows

( ) ( ) ( ) 1 1 1 1 1 1 r r r r r r 1 be be be r 1 1 1 1 r r r r , , , d d d d k i k i k i k i k i k i i i i i i + + + + + + + + + + + ∂ ∂ = + ∂ ∂ δ δ δ δ δ δ F F F δ δ δ δ δ ɺ ɺ ɺ ɺ ɺ (37)
The algorithm of the Newton-Raphson incremental-iterative procedure combined with the Newmark time integration scheme is presented in Fig. 3. The initial dynamic displacement vector is set to the static equilibrium position 0 r δ (or to other user-defined vectors) and the initial velocity vector is set to zero. When the residual vector ( )

1 r r k i+ G δ
is calculated by Eq. ( 33), a new displacement increment NR r ∆δ is produced by Eq. ( 35) and thus the displacement, velocity and acceleration vectors are updated as shown in Fig. 3.

If the ratio ( ) ( )

1 1 r r r 1
k i i t + + + G δ F
is less than a convergence criterion NR nl ε , then a converged solution is reached and the algorithm shifts to the next time step of the Newmark scheme. Otherwise the non-linear bearing forces are re-evaluated and a new Newton-Raphson iteration begins.

In order to reach the steady-state dynamic response in the periodic regime, the Newmark integration scheme must be used over a large number of time periods in the case of a harmonic external excitation of the system. This results in a huge computational time when the indicated scheme is repeated for a large number of fundamental excitation frequencies, i.e., each fundamental frequency corresponds to a large number of time periods. Computational time can be saved by computing directly the periodic solution by means of the shooting algorithm which is discussed in the next section.

Direct computation of the periodic solution (shooting method)

The following two-point boundary-value problem defined by the periodicity condition is considered, then

( ) ( ) ( ) ( ) ( ) T r r r r r r r r r r 0 , 0 , 0 with , t t τ τ   = = = - = =   H X X X X 0 X δ δ ɺ (38)
where r τ is the minimal period of the rotor response. The dependence of the system on the fundamental period ex τ of the harmonic exciting forces is assumed, i.e., r ex j τ τ = (with 1, 2,3 j = … ). Integer multiples of the fundamental period are used to compute subharmonic responses occurring, for instance, after a period-doubling bifurcation. In practice, in order to find the minimal period r τ of the response, the shooting algorithm is first run with r ex τ τ = (harmonic response) and if there is no convergence, it is then run with … consecutively (sub-harmonic response). The state-space vector r X in Eq. ( 38) represents the solution of Eq. ( 23) at the end of one period r τ for an approximated initial solution ( ) r 0 X . This vector is obtained with the non-linear Newmark time integration scheme over one period r τ (see Section 5.3.1 and Fig. 3). An appropriate linearization of Eq. ( 38) is applied by building a first-order Taylor series expansion of this equation and a Newton-Raphson iterative correction is performed, then

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) r r 1 SH SH SH r r r r r r r r r r r r r r r r 0 , 0 , 0 , 0 , 0 , 0 k k k k k τ τ τ τ τ + ∂ + = ⇒ = - ∂ X H H X H X ∆X 0 J X ∆X H X X ≃ (39) where ( ) ( ) SH 1 r r r 0 0 k k + = - ∆X X X
is the iterative increment of the state-space vector and k is the Newton-Raphson iteration. The Jacobian matrix SH r J is a function of the state-space vector ( ) r 0 X at iteration k and is given by

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) r r r r SH r r r r r r r 0 , 0 , 0 , 0 0 k k k τ τ τ ∂ ∂ = = - ∂ ∂ X X H X J X I X X (40)
The shooting method algorithm solving simultaneously Eqs. ( 23) and ( 38) is presented in Fig. 4. The approximated initial solution is set to a user-defined state-space vector ( )

1 r 0 X . When the residual vector ( ) ( ) r r r 0 , k τ H X
is calculated by Eq. ( 38), a new increment SH r ∆X of the state-space vector is produced by Eq. ( 39) and thus the approximated initial solution ( )

1 r 0 k + X is updated as shown in Fig. 4. If the ratio ( ) ( ) ( ) ( ) r r r r r r 0 , 0 , k k t τ τ = H X X X
is less than a convergence criterion SH ε , then a converged solution is reached and the algorithm is terminated. Otherwise the state-space vector

( ) ( ) r r r 0 , k t τ = X X
is re-calculated and a new Newton-Raphson iteration begins.

Linear dynamic analysis

Transient analysis

The solution technique can be significantly simplified and the computational time can be reduced if the assumption of a linearized 8-coefficient hydrodynamic bearing model is considered to be valid, i.e., the dynamic displacement vector of the rotor remains small in the vicinity of the static position 0 r δ and a first-order Taylor expansion of the dynamic fluid forces

( ) be r r , F δ δ ɺ is built by ( ) ( ) 0 be r r be r be r be r , , = - - F δ δ F δ 0 C ∆δ K ∆δ ɺ ɺ (41) 1 1 1 1 1 1 1 1 r r r r r r 0 ; ; 0 0 , 0 0 - = = = + - - r r r r r r be δ δ δ 0 δ M F F δ δ C δ K δ ɺ ɺɺ ɺ ɺ Time increment 1 i i t t t ∆ + = +
Initialization of the correction iterations (Newmark approximations)

Predictor step:

2 1 1 1 1 1 1 r r r r r r r r r ; ; 2 i i i i i i i i i t t t ∆ ∆ ∆ + + + = = + = + + δ δ δ δ δ δ δ δ δ ɺɺ ɺɺ ɺ ɺ ɺɺ ɺ ɺɺ
Evaluation of the hydrodynamic bearing forces ( )

1 1 1 1 be r r , i i + + F δ δ ɺ
Evaluation of the residual vector ( )

( ) ( ) ( ) ( ) ( ) 1 1 1 1 1 1 1 1 r r r 1 r r 1 r r 1 r 1 1 1 1 r 1 be r r , i i i i i i i i i i t t t t + + + + + + + + + + = + + - - G δ M δ C δ K δ F F δ δ ɺɺ ɺ ɺ
Evaluation of the Jacobian matrix ( )

1 r r k i + J δ
Evaluation of the correction increment ( ) ( )

1 NR 1 r r r r r k i k i + + = - J δ ∆δ G δ
Correction of the dynamic vectors at the Newton-Raphson iteration k ( )

1 1 1 NR 1 1 1 NR 1 1 1 NR 1 r r r r r r r r r 2 2 2 1 ;
;

k i k i k i k i k i k i t t ξ ξ ∆ ξ ∆ + + + + + + + + + = + = + = + δ δ ∆δ δ δ ∆δ δ δ ∆δ ɺ ɺ ɺɺ ɺɺ
Evaluation of the hydrodynamic bearing forces ( )

1 1 1 1 be r r , k i k i + + + + F δ δ ɺ
Evaluation of the residual vector ( )

1 1 r r k i + + G δ Convergence ( ) ( ) 1 1 NR r r r 1 nl k i i t ε + + + ≤ G δ F
Fig. 4. Shooting method algorithm used for the non-linear dynamic analysis of the on-board rotor. Correction of the initial solution vector at the Newton-Raphson iteration k : ( ) ( )

1 SH r r r 0 0 k k + = + X X ∆X Convergence ( ) ( ) ( ) ( ) SH r r r r r r 0 , 0 , k k t τ τ ε = ≤ H X X X
Computation of the state-space vector

( ) ( ) r r r 0 , k t τ = X X
with the algorithm in Fig. 3 Newton-Raphson increment ( )

( ) ( ) ( ) SH SH r r r r r r r 0 , 0 , k k τ τ = - J X ∆X H X
Iterations SH l ≤

Numerical simulations and result discussion

Description of the investigated on-board rotor-bearing system

The symmetric on-board rotor-hydrodynamic short bearing system is presented in Fig. 5. The origin O of the reference frame R coincides with the shaft end-point positioned at the center of the bearing # 1. The symmetric shaft is modeled with eight identical two-node Timoshenko beam finite elements, i.e., the total number of degrees of freedom is g 36 n = . The symmetric disk is located at node 5 and the mass unbalance is placed on it. The hydrodynamic short bearings # 1 and # 2 are located at nodes 1 and 9 respectively, i.e., the corresponding localized non-linear degrees of freedom are 1-4 and 33-36. The material characteristics and geometric properties of the rotor and the bearings are listed in Tables 1 and2.

The investigated rotor is subject to rotating mass unbalance excitation combined with sinusoidal rotational or translational motions of the support: the rotation around the Ox axis is given by the angular velocity 18) and data from Table 2

    × × × × = =     × × × ×     c k (43)
The equations of motion of the on-board rotor (see Eqs. ( 23)-( 32)) show that the support rotations have a parametric influence on these equations, while the support translations have only an influence on the external force vector. Due to the complexity of the rotor motion, the overall dynamic analysis of the on-board rotor is performed by means of stability charts, rotor orbits, time history responses, fast Fourier transforms (FFTs), bifurcation diagrams and Poincaré maps.

Role of the hydrodynamic bearing models in the stability analysis 6.2.1. Linearized model

In this section, the linearized equations of motion are considered (see Eq. ( 42)) and the dynamic stability of the on-board rotor is examined through the Floquet theory. The instability comes from the parametric excitation due to the sinusoidal rotation of the rotor support, see Eqs. ( 24)- [START_REF] Lee | A finite element transient response analysis method of a rotor-bearing system to base shock excitations using the state-space Newmark scheme and comparisons with experiments[END_REF]. For a certain speed of rotation Ω of the rotor, two factors affect the parametric excitations generated in the symmetric rotor: the amplitude ,

Non-linear model

According to the stability analysis in the case of linearized 8-coefficient hydrodynamic bearings, any operating point below the transition curve between the stable and unstable zones is considered to be stable. Namely, the transient dynamic analysis of the linearized equations of motion of the rotor (see Eq. ( 42)) with no external forces except the rotor weight yields the static equilibrium position 0 be C of the shaft center in the bearings for any initial conditions. Unlike the previous case, due to the non-linear bearings, the chosen initial conditions for the transient analysis of the non-linear equations of motion (see Eq. ( 23)) with only the rotor weight as external forces can lead to an unstable response instead of the static position. This phenomenon is a form of instability occurring in hydrodynamic bearings and called "self-excited vibrations" or "non-loaded bearing vibrations", see Khonsari and Chang [START_REF] Khonsari | Stability boundary of non-linear orbits within clearance circle of journal bearings[END_REF]. Fig. 7 presents the relative dynamic radial displacement of the rotor at the bearing # 2 (i.e., node 9) in the presence of a non-linear model and a sinusoidal rotation of the support of frequency 500 Hz

x Ω =

. The dynamic displacement in the time domain corresponding to the point A (Fig. 6(b)) with the initial conditions be be 2c δ (see Fig. 7(a)) settles into a stable point corresponding to the static position 0 be C but becomes unstable with the initial conditions 0 be be 2c -δ (see Fig. 7(b)). In the second case of initial conditions, the non-loaded bearing vibrations appear. Moreover, a contact between the rotor and the bearing occurs at the instant 0.012 s t = and the transient dynamic computation is stopped.

Comparison of linear and non-linear bearing models 6.3.1. Influence of the mass unbalance

All the orbits presented in this paper are made dimensionless with respect to the bearing clearance be c and are thus plotted within a unit circle.

In this section, the dynamic simulation is performed to present a comparison between rotor orbits obtained with linearized 8coefficient and non-linear hydrodynamic bearing models due to three different values of mass unbalance. Since the bearings are identical and the disk is positioned at the middle of the shaft (i.e., node 5), the orbits at both bearings (i.e., nodes 1 and 9) are identical in the case of a fixed support. These orbits illustrated in Fig. 8 show the transient motion of the shaft geometric center within the bearings from the static position 0 be C to the steady-state harmonic solution which has a period r τ equal to the inverse of the mass unbalance frequency Ω . The periodic non-linear orbits plotted in Fig. 8 are also obtained with the shooting method. As expected, they perfectly coincide with the steady-state periodic motion obtained by the transient time integration and thus validate the shooting implementation. For a small mass unbalance ( mu mu 750 g mm m r =

), the linear and non-linear orbits of the rotor exhibit very similar evolutions (see ), the dynamic displacement of the rotor becomes much greater and no longer remains in the vicinity of the static position and a high difference between the linear and non-linear orbits shown in Fig. 8(c) is noted. As a consequence, the assumption of constant damping and stiffness coefficients of the bearings is no longer valid and the linearized 8-coefficient bearing model yields inaccurate results. ). Unlike the previous section, the sinusoidal rotation of the support makes the orbits at the locations of the bearings different due to its normal and tangential accelerations. Focus is thus put on the bearing # 2 which is more affected by the support rotation.

Influence of the support motion frequency

Since the excitations due to the mass unbalance and to the support motions are sinusoidal and commensurable, the external excitation of the rotor is harmonic of period ex τ equal to the inverse of the greatest common factor of the mass unbalance frequency . In addition, the linear and non-linear orbits of the rotor in Figs. 9 and 10 are harmonic of periods r,l τ and r,nl τ equal to ex τ . As expected when the amplitude of the support motions is small, the shape of the linear and non-linear orbits in Figs. 9 and 10 changes a lot with frequency and is made more complicated, while their magnitudes are slightly modified when compared to the orbits shown in Fig. 8(b) corresponding to the fixed support. Lastly, the linear and non-linear orbits are rather similar, i.e., the dynamic motion of the rotor remains close to the static equilibrium position 0 be C and the linear bearing model can still be considered to be acceptable.

The relative linear and non-linear time history responses in the Oz direction of the rotor at the bearing # 2 due to the mass unbalance effect and to the sinusoidal rotational or translational excitations of the support are shown in Figs. 11 and12. They confirm the results provided in Figs. 9(c ), while its amplitude varies ( ). As explained previously, the excitations due to the mass unbalance and to the support motions are sinusoidal and commensurable. Thus the external excitation of the rotor is harmonic of period ex τ equal to the inverse of the greatest common factor of the mass unbalance frequency Ω and the support motion frequency x Ω or z Ω . In the case where , It can be noted that this period-doubling motion is also obtained when the transient time integration algorithm is used as a checking procedure. As expected, the linear bearing model is unable to reproduce this period-doubling motion, i.e., the linear motion of the rotor remains ex 1τ -periodic.

. Bifurcation diagrams

The response of the linear rotor system to a harmonic external excitation of period ex τ is harmonic of same period, i.e., r ex τ τ = , while the corresponding response of the non-linear rotor system can be non-periodic or can have a different period due to bifurcations as seen in Section 6.3.3. Such changes in the non-linear dynamic regime of the system can be investigated by means of diagrams of bifurcation of periodic solution. 3 and4 present a summary of the different dynamic regimes observed on the bifurcation diagrams of the on-board rotor-bearing system excited by sinusoidal rotational or translational motions of its support as a function of the bifurcation parameters. Influence of the initial conditions on the orbits of the rotor in the presence of a mass unbalance (mmurmu=1500 g mm) and sinusoidal translations of the support: zO=ZOcosΩ z t with ZO=5.2×10 -5 m and Ω z =80 Hz.

Sensitivity to initial conditions

A very small change in the initial conditions for a non-linear rotor system can make an important change in its corresponding dynamic trajectory. This phenomenon is called "sensitivity to initial conditions" and is one of the main ways in which a non-linear system can be recognized. This is demonstrated in Figs.

Conclusions

In order to predict and analyze numerically the role of the sinusoidal rotational or translational motions of the support in the dynamic behavior of an on-board rotor-non-linear hydrodynamic journal bearing system, a new finite element model based on the Timoshenko beam theory is presented in this paper. This model includes the effects corresponding to the rotary inertia, the gyroscopic inertia, the shear deformation of shaft, the geometric asymmetry of shaft and/or rigid disk as well as six types of deterministic motions (rotations and translations) of the rotor rigid support. In addition, the rotational effects and the geometric stiffening effects associated with the centrifugal stressing caused by the support rotations are taken into consideration.

The derivation of the equations of motion displays that the rotational motions of the rotor support bring not only a parametric contribution in the damping and stiffness matrices of the rotor but also additional external force vectors, while the translational motions of the support only add external force vectors. A robust algorithm combining the Newton-Raphson incremental-iterative procedure with the Newmark time integration scheme is used to obtain the non-linear transient dynamic motion of the on-board rotor. Moreover, a robust shooting algorithm is used to obtain the periodic non-linear responses of the rotor and is validated against the transient time integration algorithm.

When the rotor is excited by a sinusoidal rotation of the support, it is shown that a dynamic instability of the rotor can occur because of time-varying parametric coefficients in the linearized equations of motion and that the linearized hydrodynamic bearings have an influence greater than the rigid bearings on the size of the instability regions of the rotor.

It is noted that the shape and the magnitude of the rotor orbits can be significantly affected by the frequency and amplitude of the support motions. As expected, the linearized bearing model gives accurate results only for small amplitudes of the support motions.

Also, this simplified bearing model cannot account for non-linear phenomena and does not permit to predict the change of dynamic regime such as period-doubling motion and transition to chaos. Conversely, these phenomena are observed by means of diagrams of bifurcation of dynamic responses when the non-linear bearing model is considered.

It is observed that depending on the amplitudes of the sinusoidal motions of the support, the non-linear rotor system can exhibit sub-harmonic, quasi-periodic and chaotic motions. Lastly, it is well demonstrated that different initial conditions for the non-linear transient dynamic motion can lead to different responses of the rotor in the presence of the support motions.

Appendix A

The analytical expressions of the components of the damping and stiffness matrices be be R c and be be R k in the polar frame linked with the bearing be R are given by be be be be be be r r be be be be be be be be

;

tata tara tata tara R R rata rara rata rara c c k k W W c c c c k k Ω     = =         c k (A.1) with ( ) ( ) ( ) ( ) ( ) ( )
( )
be be be be be be be be 0 2 0 2 0 2 be be be 0 0 0 0 be be be be be be be be be be be be 0 0 2 0 0 2 be be be be 0 0 2 2 be be be

c k c k c k c k k f k f k f k f f π ε ε π ε ε ε ε ε ε ε ε ε ε ε π = = - = - = - - + + - = = = = - - = + ( ) ( ) 1 0 2 2 be 1 ε - - (A.2)

Fig. 1 .

 1 Fig. 1. (a) Components and frames of reference of the on-board rotor model, (b) Euler angles for the transformation between the frame linked with the moving support R and the local frame R l .

ω

  x y z is the non-inertial frame of reference attached to the moving rigid support, z is the non-inertial local frame of reference (moving with the rotor during its operation) attached to the mass center of the rigid disk or the geometric center of a section of the deformable shaft. The origins of the frames of reference g R , R and l R are g O , O and l O respectively. The elastic line of the non-deformed shaft is along the Oy axis on which the mass center of the disk is located. Moreover xOy and yOz are both planes of symmetry for the disk and the non-deformed shaft. Thus Ox , Oy and Oz are principal axes of inertia.The translational motions of the rotor support are described by the coordinates ( ) frame attached to the support R . The rotational motions of the rotor support are described by the components of the rigid support R with respect to the ground g R projected in the frame R .

Iv

  are the principal moments of inertia of the disk mass about the Ox , Oy and Oz axes respectively. In addition, seperate the effects due to the mean moment of inertia of the disk mass and those due to the moment of inertia modeling the disk geometric asymmetry. Thus when the disk is asymmetric, the moments of inertia d is found through the absolute position vector g l d O O of the disk mass center positioned at the arbitrary abscissa d y with respect to the frame R . The latter is observed from the ground and expressed in the frame fixed to the support R by

  the components of the translational velocity vector of the disk center. Considering the components

  φ ɺ are replaced by t Ω and Ω respectively. The translational displacements () their time derivatives are very small. Thus the trigonometric functions are replaced by their Taylor series expansion and the resulting expressions are limited at the order 2

  respectively the mean moment of inertia of the cross-sectional area and the moment of inertia characterizing the asymmetry of the shaft.

  the absolute position vector g mu O P of the mass unbalance. The latter is measured from the ground and written relative to the frame linked with the support R as

Fig. 2 (

 2 Fig. 2(a) shows a simple representation of a hydrodynamic bearing which consists of a fixed journal (1) including a rotating shaft (2) and a lubricant (3) separating these two components. The points O and lO are respectively the bearing center and the shaft geometric center. The radius, length and clearance of the bearing are respectively be r , be l and be be sh

  speed of rotation Ω of the rotor and for a constant static load r W created by the rotor weight, the shaft geometric center l O in the bearing occupies a static equilibrium position center in the journal and the static attitude angle 0 be ϕ between the r W load direction and the line of centers 0 be

  . As the bearing contribution appears as forces be F external to the rotor, their virtual work be W δ must be established to apply the Lagrange's equations() the dynamic displacements be δ of the shaft elastic line are assumed to be small in the vicinity of the static position 0 be δ , the linear analysis of a bearing is applied by constructing a first-order Taylor expansion of the dynamic fluid forces ( )

(

F

  δ 0 is related to the static displacement vector 0 be δ and must equilibrate the constant load r W imposed by the rotor weight. The damping and stiffness matrices be c , be k of the linearized 8-coefficient hydrodynamic bearing (i.e., four damping and four stiffness coefficients) (see Fig. 2(b)) depend only on the static position 0 be δ . Actually, they are easier to obtain in the polar frame be R , i.e.,

Fig. 2 .

 2 Fig. 2. (a) Schematic diagram of a hydrodynamic bearing, (b) Linearized damping and stiffness coefficients of a bearing model.

the

  Oy axis. This finite element is of density

  motions of the nodes i n and 1 i n + in the planes

  point init C due to bending respectively along the Ox and Oz axes of the frame R as well as the corresponding rotations i esh θ and i esh ψ respectively around the Ox and Oz directions are formulated as a function of the nodal displacement vectors n,

F

  are the shape function vectors of a Timoshenko beam finite element. They are thirddegree and second-degree polynomial equations, functions of the dimensionless variable y (

F

  δ δ ɺ of the bearing is opposite to the rotor motion.

Fig. 5 .

 5 Fig. 5. On-board rotor-hydrodynamic bearing system configuration.

Table 1 .Table 2 .

 12 Main characteristics of the on-board rotor. Density of both disk and shaft material Radius, thickness and location of the disk Radius and length of the shaft Young's modulus and Poisson's ratio of the shaft Shear correction factor of the shaft Mass unbalance phase ρd=ρsh=7800 kg/m 3 rd=0.15 m, ed=0.03 m, yd=0.2 m rsh=0.04 m, lsh=0.4 m Esh=2×10 11 N/m 2 , νsh=0.3 κsh x = κsh z =0.8864 [38] ηmu=0° Main characteristics of the hydrodynamic short bearings. Radius and length of the bearings Locations of the bearings Radial clearance of the bearings Oil film dynamic viscosity rbe=0.04 m, lbe=0.01 m ybe.# 1=0 m, ybe.# 2=0.4 m cbe=2×10 -4 m µbe=288×10 -4 Pa s

.

  unbalance frequency) and its first four natural frequencies for the considered speed of rotation and a fixed support in the case of linearized bearings are 33.6 Hz, 52.7 Hz, 524.8 Hz and 1573 Hz. The relative coordinates of the static equilibrium position 0 be C of the shaft geometric center in the fluid film bearings are given by The damping and stiffness matrices of the linearized 8-coefficient bearing model obtained with Eq. (

ωFig. 6 (.,Fig. 6 .Fig. 7 .

 667 Fig. 6(a) and (b) compare the stability charts investigated in the presence of the rigid short bearings (for which the first four natural frequencies of the rotor at 1200 rpm Ω =

Fig. 8 (

 8 Fig. 8(a)). According to Fig. 8(b), the chosen mass unbalance ( mu mu 1500 g mm m r =

Figs. 9 Fig. 9 ,5 6 2

 996 Figs. 9 and 10 display a comparison between linear orbits of the rotor and non-linear ones (obtained by the transient Newmark integration and the shooting algorithm) at the bearing # 2 (node 9) in the presence of the mass unbalance excitation ( mu mu 1500 g mm m r =

Fig. 8 .Fig. 9 .Fig. 10 .

 8910 Fig. 8. Comparison of linear and non-linear orbits of the rotor at both bearings for three different values of mass unbalance: mmurmu= (a) 750 g mm, (b) 1500 g mm, (c) 7500 g mm. Linear bearing model Non-linear bearing model

Ω

  and the support motion frequency x Ω or z Ω . For example, one period of the first external excitation of the rotor corresponds to one period of the mass unbalance excitation ( 20 Hz Ω =) and to four periods of the support motion for x

  ) and 10(c), i.e., these responses exhibit similar evolutions and their periods r,l τ and r,nl τ are equal to the period ex τ of the external excitation of the rotor. Moreover, the FFTs of the linear response (see Figs. 11(a) and 12(a)) show two frequency components due to the mass unbalance excitation (20 Hz) and the support motions (200 Hz), while the FFTs of the non-linear response (see Figs. 11(b) and 12(b)) show additional frequency components which correspond to multiples of 20 Hz (super-harmonics) relative to the harmonic motion of the rotor. The peaks at these additional frequencies are relatively small compared to those of 20 Hz and 200 Hz. This validates the assumption of small displacements of the rotor in the vicinity of the static position 0 be C . 6.3.3. Influence of the support motion amplitude

Figs. 13 Fig. 11 .Fig. 12 .

 131112 Figs.[START_REF] Zhao | Stability and response analysis of symmetrical single-disk flexible rotor-bearing system[END_REF] and 14 introduce a comparison between linear orbits of the rotor and non-linear ones (obtained by the transient Newmark scheme and the shooting algorithm) at the bearing # 2 located at node 9 for the mass unbalance excitation ( mu mu 1500 g mm m r =

Figs. 13

 13 Figs.[START_REF] Zhao | Stability and response analysis of symmetrical single-disk flexible rotor-bearing system[END_REF](a) and 14(a), the linear and non-linear orbits of the rotor are harmonic of periods r,l τ and r,nl τ equal to the period ex τ but they are not similar. This justifies that the assumption of small displacements of the rotor in the vicinity of the static equilibrium position 0 be C is not verified anymore and the validity range of the linearized bearing model is more affected by the amplitude than by the frequency of the support motion. For larger amplitudes of the support motions ,

Fig. 13 .Fig. 14 .

 1314 Fig. 13. Comparison of linear and non-linear orbits of the rotor at the bearing # 2 for two different amplitudes of sinusoidal rotations of the support: ω x =ω x,a cosΩ x t with Ω x =80 Hz and ω x,a = (a) 8×10 -2 rad/s, (b) 12×10 -2 rad/s.

Figs. 15 )Fig. 15 .Fig. 16 .

 151516 Figs. 15 and 16 illustrate the linear and non-linear time history responses along the Oz axis of the shaft geometric center within the bearing # 2 for the mass unbalance effect and the sinusoidal rotational or translational excitations of the support. They verify the results pointed out in Figs. 13(b) and 14(b), i.e., the non-linear response is ex 2τ -periodic because its period r,nl τ is two times higher than the period r,l τ of the linear response. In addition, the FFTs of the linear response (see Figs. 15(a) and 16(a)) display two frequency components due to the mass unbalance (20 Hz) and the support motions (80 Hz), while the FFTs of the non-linear response (see Figs. 15(b) and 16(b)) highlight much more frequency components and contain harmonics, sub-harmonics and superharmonics relative to the period-doubling motion of the rotor. Fig. 17 gives the linear and non-linear orbits of the rotor within the bearing # 1 (node 1) due to the mass unbalance excitation ( mu mu 1500 g mm m r =

Fig. 17 .

 17 Fig. 17. (a) Linear, and (b) non-linear orbits of the rotor at the bearing # 1 in the presence of a sinusoidal rotation of the support: ω x =ω x,a cosΩ x t with ω x,a =12×10 -2 rad/s and Ω x =80 Hz.

6 . 4 .

 64 Route to chaos in the rotor-non-linear hydrodynamic bearing system 6.4.1

Figs. 18

 18 and 19 present the bifurcation diagrams of the rotor at the bearing # 2 (node 9) for the mass unbalance effect ( mu mu 1500 g mm m r = ) and the sinusoidal rotational or translational motions of the support of frequency x Z of the support motions are used as bifurcation parameters to provide essential information about the Poincaré sections. The steps used for the bifurcation parameters , x a ω and O Z are 0.05×10 -2 rad/s and 0.05×10 -5 m respectively. The period of the harmonic external excitation of the rotor is ex 1 20 s 0.05 s τ = = which corresponds to the inverse of the greatest common factor of the mass unbalance frequency numerical simulations are performed with 1000 rotor revolutions (i.e., ex 1000τ ) and 512 time steps per period ex τ . The data of the first 500 revolutions of the rotor are not used in the bifurcation diagrams in order to exclude the transient regime. The Poincaré sections are obtained by taking successive intersections of the x coordinates of the dynamic displacements of the rotor with the instants corresponding to the multiples of the period ex τ , i.e.,

2 Fig. 18 .Fig. 19 .

 21819 Fig.[START_REF] Kishor | On the dynamic analysis of a rigid rotor-gear pair-hydrodynamic bearing system[END_REF]. Bifurcation diagram of the rotor at the bearing # 2 for a sinusoidal rotation of the support: ω x =ω x,a cosΩ x t with Ω x =80 Hz.

Fig. 21 .Fig. 22 .

 2122 Fig. 21. Orbits and Poincaré maps of the rotor at the bearing # 2 due to sinusoidal translations of the support: zO=ZOcosΩ z t with Ω z =80 Hz and ZO= (a) 7.3×10 -5 m, (b) 8.6×10 -5 m, (c) 9.5×10 -5 m.

  Fig. 23. Influence of the initial conditions on the orbits of the rotor in the presence of a mass unbalance (mmurmu=1500 g mm) and sinusoidal translations of the support: zO=ZOcosΩ z t with ZO=5.2×10 -5 m and Ω z =80 Hz.

. 23 , 5 5

 235 [START_REF] Samali | Random vibration of rotating machines under earthquake excitations[END_REF] and 23 which provide different orbits of rotor and their projected Poincaré maps at the bearing # 2 (node 9) produced due to different initial conditions for the transient dynamic motion.In Fig.22, the rotor support is excited by a sinusoidal rotation of amplitude , orbit corresponds to a sub-harmonic motion ( ex 2τ -periodic), see Fig.22(a), when the initial conditions are set to the coordinates of the static equilibrium position corresponds to a harmonic motion ( ex 1τ -periodic), see Fig. 22(b), when the initial conditions are slightly modified and set to the coordinates [ ] the rotor support is excited by a sinusoidal translation of amplitude to a quasi-periodic motion, see Fig. 23(a), when the initial conditions are set to corresponds to a sub-harmonic motion ( ex 4τ -periodic), see Fig. 23(b), when the initial conditions are set to [ ]

  )

	d,sh,su	t	d,sh,su u		O x		2	z	O	y	2	O y	z		z	O	y	x z	O y	z	x y		O x	y	2	z	2
			d,sh,su w	z	O	2	O y	x	2	O x	y	y	O	x	y z	O x	y	x z	z	O	x	2	y	2
			yw		x		y z		yu		z			x y		x		y z			z
			d,sh,su								d,sh,su					d,sh,su			d,sh,su

Table 3 .

 3 Types of motion identified in the bifurcation diagram due to the sinusoidal rotation of the support. Type of the rotor motion Support rotation amplitude ω x,a (10 -2 rad/s)

	1τex-periodic	[1 -11.7]
	2τex-periodic	[11.75 -12.05], [12.4 -12.8]
	4τex-periodic	[12.1], [12.85]
	Quasi-periodic

[12.15 -12.35

],

[12.9 -13]

). The global damping and stiffness matrices be C , be K deduced from the bearing linearization (see Section 4.4) are generated by the assembly of the elementary matrices. These matrices are kept constant all along the dynamic simulation. Substituting Eq. (41) into Eq. ( 23), the equations of motion of the on-board rotor-bearing system become linear as follows

The transient dynamic motion of the rotor nodes is obtained by solving Eq. ( 42) by means of the implicit Newmark time-step integration algorithm based on the average acceleration. The displacement vector 0 r δ of the static equilibrium position 0 be C is used to initialize the transient dynamic problem. The final integration time is selected such that the transient effects have disappeared and the steady-state regime has been reached. Lastly, it is necessary to establish precisely the validity range of the linear bearing model by comparing its results with those of the non-linear bearing model.

Stability analysis

The motion stability of the linear rotor-bearing system is investigated by applying the Floquet theory for one period pe τ of the periodic parametric excitation, see Dufour and Berlioz [START_REF] Dufour | Parametric instability of a beam due to axial excitations and to boundary conditions[END_REF] and Dakel et al. [START_REF] Dakel | Steady-state dynamic behavior of an on-board rotor under combined base motions[END_REF]. In other words, by calculating the so-called "Floquet transition matrix" (or monodromy matrix) during one period pe τ of the parametric excitation and by finding its eigenvalues, the zones of instability are identified, i.e., if the modulus of at least one of these eigenvalues is greater than 1, the system is unstable. As reported in [START_REF] Peletan | A comparison of stability computational methods for periodic solution of nonlinear problems with application to rotordynamics[END_REF], the Floquet transition matrix can be very efficiently computed using the linear version of the Newmark algorithm.

Sub-harmonic, quasi-periodic and chaotic bifurcation responses

In this section, some orbits of the rotor and their projected Poincaré maps at the bearing # 2 (node 9) corresponding to values of interest of the bifurcation parameters

When the amplitudes ,

x a ω of the support rotation are successively equal to 12.1×10 -2 rad/s (Fig. 20(a)), 12.15×10 -2 rad/s (Fig. 20(b)) and 13×10 -2 rad/s (Fig. 20(c)), the rotor orbits highlight a sub-harmonic motion ( ex 4τ -periodic) and quasi-periodic motions respectively. When the amplitudes O Z of the support translation are successively equal to 7.3×10 -5 m (Fig. 21(a)), 8.6×10 -5 m (Fig. 21(b)) and 9.5×10 -5 m (Fig. 21(c)), the rotor orbits display a chaotic motion, a sub-harmonic motion ( ex 4τ -periodic) and a quasiperiodic motion respectively.