A new model for analyzing the vibration behaviors of rotor-bearing system

Journal Pre-proof

A New Model for Analyzing the Vibration Behaviors of Rotor-bearing System Zi Wang , Caichao Zhu PII: DOI: Reference:

S1007-5704(19)30449-6 https://doi.org/10.1016/j.cnsns.2019.105130 CNSNS 105130

To appear in:

Communications in Nonlinear Science and Numerical Simulation

Received date: Revised date: Accepted date:

10 April 2019 2 November 2019 20 November 2019

Please cite this article as: Zi Wang , Caichao Zhu , A New Model for Analyzing the Vibration Behaviors of Rotor-bearing System, Communications in Nonlinear Science and Numerical Simulation (2019), doi: https://doi.org/10.1016/j.cnsns.2019.105130

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier B.V.

Highlights      

Adynamic model forrotor-bearing system with strong nonlinearity is proposed. Two kinds of bearing stiffness are classifiedbased on theirdifferent physical meanings. Thenumericalcomputation processbased onimproved Runge-Kutta methodis described in details. One type of cylindrical roller bearing in introduced for case study. Strong nonlinear phenomenon and behaviorsare observedand discussed.

Title: A New Model for Analyzing the Vibration Behaviors of Rotor-bearing System Authors:Zi Wang, Caichao Zhu Affiliation:State Key Laboratoryof MechanicalTransmissions, Chongqing University, P.R.China, 400044 Contact Information:Zi Wang: [email protected] Caichao Zhu : [email protected], Corresponding author. Abstract: Considering the strong non-linear properties of rolling element bearings (REBs), a new model for the dynamic analysis of rotor-bearing system, especially for the heavy-load and high-speed condition, is proposed in this work. Two types of bearing stiffness are classified and analyzed with physical meanings based on the load-deflection curves, which can full preserve the geometric and kinematic nature of REBs for the vibration model. The improved numerical integration process based on 4-th order Runge-Kutta method is described in details. One type of cylindrical bearing with internal clearance is brought in for case study. Complicated non-linear vibration behaviors caused by REBs instinct properties are observed. The effects from speed fluctuation, unbalance forces, external loads on the vibration responses of the system are investigated. This work may provide some references for the structure design, vibration prognosis and other engineering applications of the rotor-bearing system. Keywords: rolling element bearings; bearing stiffness; non-linear; rotor vibration;

1 Introduction Bearings are the most important components in rotary systems. As the requirements from modern machines for the higher-speed, heavy-load, high accuracy and longer life span are emerging, the vibration problems caused by inherent geometric and kinematic characteristics from bearings are gradually displayed in applications such as wind turbines, aircrafts, machine tools and etc.. Because the bearing failure is among the most primary causes for the rotary machine breakdown [1-3], massive researches on health assessment techniques were conducted recently for the detection, diagnosis and prognosis of rolling element bearings (REB). The existing maintenance techniques of rolling element bearings are mostly focused on exploring the signal processing methods based on vibration characteristics of rotary transmission system with defected REBs [2, 4]. However, the investigations on the nonlinear dynamics model for REB-rotor system, which are useful and practical for understanding the origins of nonlinearity occurring in the vibration behaviors of rotary systems, are progressing slowly.

The REBs excite varying compliance vibration even with defect-free geometrical structure, which is because the rolling elements carrying loads vary along the cage rotation. From the very beginning, Sunnersjö [5] studied the varying compliance vibrations of rolling bearings for the bearing with radial clearance under radial load by computer simulation. Subsequently, M.F.While [6] defined the bearing stiffness as the rate of change of load with deflection and calculated the stiffness fluctuation affected by rolling elements movements though the load zone. Fukata S., Gad E and their colleagues [7] found sub-harmonic, super-harmonic and chaos behaviors might appear in the vibration of ball bearing system which is excited by varying compliance frequency, which is also called as ball pass frequency. Lim T.C. and Singh R. [8, 9] studied the vibration transmission through REBs step by step; firstly, they proposed a mathematical model for the comprehensive bearing stiffness of REBs. Then a new model which combines the lumped parameter method and finite element approach for the vibration transmission of a shaft-bearing-plate-mount system through the comprehensive bearing stiffness matrix was proposed. And results were proved to be reasonable good when comparing with experimental data. Internal clearance is necessary for heat expansion of REBs. Based on the analytical model mentioned above, the researchers later developed the model which considers the internal clearances. Tiwari M. ,Gupta K. and Prakash O. [10] observed that the radial internal clearances contribute to the periodic, subharmonic and chaotic vibration behaviors of the rigid rotor supported by deep groove ball bearing. They also conducted the experiments [11] for investigating the effect from internal clearance on the rigid rotor vibration behaviors and results shows that bearings with large clearances have stronger super-harmonics and more sub-harmonics. Similar effect from clearance on the sub-harmonic response has been obtained by Upadhyay S.H., Harsha and Jain S.C. [12] and they categorized the nonlinear dynamic responses of the bearing with clearances as harmonic, sub-harmonic, quasi periodic and chaotic motions. The experiments described in [13] discovered two routes to chaos of vibration behaviors for the ball bearing with internal clearance. Kappaganthu K. and Nataraj C. [14] proposed a detailed bearing model for considering the clearances and found that the bifurcations and chaos behaviors are strongly influenced by clearances. From the conclusions of Zhang Z., Chen Y. and Cao Q. [15], even for a small clearance, the coupling nonlinear effects with the influences from Hertzian contact and varying compliance can cause quasi-period, chaotic and other complicated vibration resonances. A dynamic model was proposed [16] for the rotor-bearing system and the increase of clearances are found to cause more complicated dynamic behaviors. The unbalance mass effect of the rotor plays an important role in vibration motions of the bearing rotating systems [11, 12], which is one of the most common excitations for the high-speed rotating machines and cannot be avoided in practical engineering. The formulation for the rotor dynamic response subjected to the unbalance excitations was presented in detail[17] and periodic and

quasi-periodic resonances at the combinations of unbalance frequency and base excitations are observed. Zhang X, Han Q. and their cooperators [18] found that the width of instability regions can be changed by the unbalance force variation. Sharma A., Amarnath M. and Kankar P.K. [19] proposed an analytical model of a rigid-unbalanced-rotor –bearings system for exploring the nonlinear behaviors. Similarly, the resonant frequencies are observed to be associated with bass pass frequency and rotor rotation frequency. In addition, the high speed effect is also a key influence factor for bearing model properties [20, 21]. Most the techniques described in the above publications for calculating contact deformation between raceways with rolling elements are based on Hertz contact theory, with the variation of the number of rollers in load zone, which will result in the nonlinear relationships between the loads with deflections. Researchers [6, 8, 22] both observed this kind of nonlinearity instinct to the REBs both theoretically and experimentally. The determination of bearing stiffness keeps being an important issue for researchers for decades [8, 21-24]. While for large amount of publications [12, 14, 17, 19, 25], there is only restore force in the dynamic model for presenting the transmission capacity of bearing. For this kind of method, lots of mathematical derivations are needed for obtaining the bearing varying compliance properties. The stiffness matrix calculated by finite difference method is only proper for the low-load work condition or for the bearings with linear load-deflection curves. It is necessary to classify the definition for the general physical meanings of bearing stiffness, no matter the load condition or the types of bearings. Based on the descriptions above, the vibration behaviors caused by the geometric and kinematic characteristics of REBs with internal clearances are extra complicated. The varying compliance, internal clearance, mass unbalance, speed effect and nonlinearity between loads with deflections have both been proven to make dramatic effects on the vibration motions of rotor-bearing systems through the investigations mentioned above. However, those model or techniques are mostly based on an oversimplified Hertzian contact model even with overcomplicated mathematical formulations for achieving the varying compliance feature [8, 14-22]. Also the operation conditions are limited: some models are only prepared for low or no external loads [10, 12, 13, 15] and lots of works ignore the high speed effects and unbalance force [5-7, 9, 10, 14-16]. The effects from nonlinearity of relationships between loads with deflections are totally neglected [5, 7, 10]. For the heavily nonlinear system with complicated unpredicted excitations, for which the principal of superposition is not applicable any more, that makes the solutions for the equations of motion cannot be obtained through analytical approaches [10, 26-28]. This work proposes a REB-rotor dynamic analysis model, which is suitable for all types of bearings and also appropriate for high-load condition. The bearing forces are calculated by two kinds of bearings stiffnesses with different physical meanings, based on the nonlinear

load-displacements curves. The unbalance forces, not only from the static eccentricity but also from the translational shift progress in the dynamic progress, are also considered. The numerical integration process for this kind of nonlinear vibration system is described in detail. One type of cylindrical roller bearing is brought in for the case study. The effects from variation of speed, static eccentricity, external loads and unbalance force on the dynamic behaviors are studied.

2. Dynamic Analysis Model for Rolling-Element-Bearing-Rotor System 2.1 Equations of motion Y X

s O

Z

Fig.1 Coordination definition of rotor bearing system

The configuration of rotor supported by a rolling element bearing is commonly seen in the application of machine spindle, aero-engines or other transmission systems. The rotor is assumed to be rigid in this work, assembled to the inner raceway of bearing through interference fit. The outer raceway of bearing is fixed on a housing, which is attached to the coordinate O-XY showing in Fig.1. The coordinate attached to the mass center of rotor is O’-XY, described in Fig.2 The kinematic and geometric analysis of rolling element bearing is complex that involving four components, which are rolling elements, inner and outer raceways and cage. In addition, for the high-speed operation condition, the internal clearance in the design process is necessary for the thermal expansion. The consideration of strongly non-linear and time-dependent of load-deflection relationship can lead to the increase in complexity of the rotor-bearing system vibration problems. This work is focused on the effect from instinct properties of rolling element bearings on rotor vibration problems for severe working condition with high-speed and heavy-load. In order to isolate the vibrations behaviors caused by the bearing and reduce the complexity of computing process, the following assumptions are made below: (1) There is no manufacturing or assembling error for this system which is assumed to be defect-free. (2) The rolling elements are equi-spaced over the pitch circle. (3) There is no impact between cage with rolling element, raceways. The equations of motion for the rotor-bearing system can be written in Eq.(1).

Mu  (C  G)u  Fb (u)  W  Fe (u)

(1)

Where M=diag[m m m Ib Ib] denotes the mass matrix and m and Ib are mass and moment of inertia of the rotor, respectively; u is the displacements of the system and can be expressed as u=[ux uy uz θx θy]T which represents the five degree-of-freedom motions including plane motions on X and Y directions, axial direction Z, and angular distortion around axes of X and Y; C is the Rayleigh damping matrix for eliminating transient response, which is determined by: C   M   K where α and β are the proportional constants to mass and average stiffness matrix. K is the average effective bearing stiffness matrix. In which G and W are the gyroscopic and external loads matrices, which can be written by: W  Wx

0 0  G  0  0 0 

Wy

0 0 0 0 0

Wz

Mx

0 0 0 0 0 0 0 0 0  I p s

M y 

T

     I p s  0  0 0 0

Where Ip is the polar moment of inertia of rotor and ωs is the rotation speed of rotor; Fe is the unbalance force; Fb denotes the bearing support force.

3. Model Excitations Analysis 3.1 Unbalance Force In practical engineering, the mass center is always out of the geometric center of the rotor; in addition, the rotor plane motions caused by bearing elasticity will shift the center of mass of rotor when applied external loads. For the static and quasi-static state, the terms with time derivations are ignored and this kind of eccentricity will not draw extra effects of deflection or force on the system. However, for the dynamic analysis, the centrifugal force is obvious due to the mass eccentricity for high-speed condition which is also an excitation for the vibration response. The calculation process for the unbalance force during dynamic response process, caused by mass eccentricity e is showed as follows: Fe   Fex

Fey

0 0 0 

T

Fex  ms2 e sin 

(2)

Fey  m e cos  2 s

In which γ is the phase angle showed geometrically in Fig.2, given by:

 =arctan

u y    cos(s t ) ux    sin(s t )

(3)

Fig.2 Calculation of eccentricity

Fig. 2 describes the calculation scheme of eccentricity OO’, in which O’ represents the mass center of rotor. The dynamic eccentricity of the rotor at an arbitrary time step is denoted by e, which equals to the length of OO’. The length of O’D is denoted by term ε, which presents the static eccentricity due to the alignment or manufacture error. Thus the dynamic eccentricity e is gotten by:

e  ux2  u y2   2  2 ux sin(st )  u y cos(st ) 

(4)

3.2 Bearing Support Force The support force transferred from rolling element bearings is the uppermost important excitation for the system. The varying compliance and nonlinear load-deflection characteristics of a rolling element bearing with internal clearance on Y direction are showing in Fig.3, in which the highly nonlinear properties are found upon the load-deflection distributions. The similar phenomenon has been proven through theoretical derivation [21]. Based on the load-deflection curve for an arbitrary position over the ball pass cycle, the bearing support force at a transient dynamic equilibrium state can be obtained through static equilibrium point and the small motions about the equilibrium points, which is expressed in Eq.(4).

Fb (ue  Δu)  Fb  ue   Fb  ue , Δu 

(5)

Where ue denotes the displacement at the certain instantaneous equilibrium point; Δu is the small motion, which is also called small perturbation from equilibrium points and Fb is the disturbance force around the normal load zone. The curves of the load-deflection of bearing are supposed to be smooth enough and disturbance force can be expressed through Taylor series in Eq.(6). Fb (ue )  Δu

In which

Fb u

 O(Δu 2 )

(6)

u  ue

O(Δu2 ) is the truncation error and the superscript 2 presents this kind of physical model

in this work is second-order accurate.

External load W y / kN

10 8 6 4 2 u 0y  s1

0

30

35 40 45 Deflection u y / μm

50

Fig.3 Load-deflection curve on Y direction

Two different types of stiffnesses are classified for capturing the strong nonlinearity inside the load-deflection distribution of the REBs, which are defined as mean stiffness k and alternating stiffness k . 1) Mean stiffness definition The bearing mean stiffness is calculated from the definition that is the magnitude of the elastic deflection in the bearing under load, written as: k

W W  diag  x u - u0  uˆx

Wy uˆ y

Wz uˆz

Mx ˆ x

My    diag[k x ˆy 

ky

kz

k x

k y ]

(7)

In which W is the applied external load; u0 is the no-load deflection. Thus the bearing mean stiffness only can be activated while the deflection exceeds the no-load deflection on the corresponding directions. The no-load deflection is examined or calculated when the external load is tending toward to nil, which is a crucial parameter for judging whether the rolling elements are in contact with raceways. While under the static analysis that the centrifugal forces is neglected, for the starting point of ball pass cycle in the static state, the no-load deflection at the starting position over one ball pass cycle, which the position is showed and defined in the Fig.1, is calculated as:

 bc  u0  s1  0 0 0 0 2  

(8)

In which bc is defined as the internal radial clearance. The bearing mean stiffness presents the capacity for transferring load to the corresponding directions. Thus the bearing load Fb transferred from bearing to the rotor can be activated while the deflection vector over the no-load deflection and the following relationships can be gotten:

Fb  k  u - u0 

(9)

2) Alternating stiffness definition The bearing alternating stiffness is defined as the rate of change of load with deflection with respect to a normal deflection vector, is obtained by:

 Wx  u  x  Wy   u x  W W z k    u x u   M x  u x   M y  u  x

Wx u y

Wx u z

Wx  x

Wy

Wy

Wy

u y

u z

 x

Wz u y

Wz u z

Wz  x

M x u y

M x u z

M x  x

M y

M y

M y

u y

u z

 x

Wx   y   Wy   k xx    y   k yx Wz   k  y   zx   M x   k x x  y   k x   y M y   y 

k xy

k xz

k x x

k yy

k yz

k y x

k zy

k zz

k z x

k x y

k x z

k x x

k y y

k y z

k y x

k x y   k y y   k z y   k x y   k y y 

(10)

The bearing alternating stiffness shows the load fluctuation capacity around the normal equilibrium point. Discarding the remainder of Eq.(6), the disturbance bearing force is achieved by the alternating stiffness through:

Fb  Δu  k

(11)

4. Numerical Method Description Upon the description of section 3, the system is strong nonlinear with complicated unpredicted excitations. The techniques for solving the Eq.(1) are totally different for the linear systems that there are no analytical solutions for obtaining the response from nonlinear systems. Thus, an improved numerical integration method will be presented for the rolling element bearing-rotor system, which also works for high-speed and heavy-load operation condition. Firstly, the state variables x   x1

x2 … x10  , are introduced for transforming the Eq.(1) in T

generalized coordinates to equations in state space, which are obtained by Eq.(12) u   x1

x2

x3

x4

x5 

u   x6

x7

x8

x9

x10 

T

(12)

T

The acceleration vector u is denoted by nonlinear functions f   f1

f2 …

f10  as u  f (u, u) . T

For expressing the approach in a simpler way, the notations below are brought in: xi+5=Xi , fi= Xi , in which i=1,2,…,5 presents the five degrees of freedom for the system. Then the equations of motion containing 5 DOFs can be transformed into 10 first-order state equations: x(t )  X[x(t )]

(13)

The values of equilibrium points for each time step can be obtained through the numerical

integration methodology. For this kind of method, the finite difference approximations are introduced for replacing the first order of derivation term x(t ) . We introduce the small time increment t and the solutions for the Eq.(13) is expressed by the Taylor series: x  t  t   x(t )  tx  (t )  1

t 2 2 t N N t N 1 N 1 x (t )  ...  x (t )  x ( ) 2! N!  N  1!

(14)

In which the superscript inside {} is the derivative number with respect to time and the last term is the remainder with range t    t  t . The Eq.(14) can be rewritten as Eq.(15) while introducing the notation s for time discretization.

x  t  t   x(ts 1 )  x(s  1)

(15)

Thus the state vector x at time step s+1 is gotten by:

t N { j 1} x( s  1)=x( s)+ f  s   O( ) j 1 j ! N

(16)

For an arbitrary time step s in a dynamic process, the forces Fb and Fe which are instantaneous equilibrium is obtained by the state vector x(s). Physically, the varying compliance and the eccentricity of rolling element bearing will fluctuate Fb and Fe over ball pass cycle, which are the excitations for this dynamic model. The strong nonlinear relationships between the bearing force with the displacement vector will both consider the bearing equilibrium force Fb and the disturbance term Fb , which is presented by Eq.(17) in the integration process. Fb ( s)  u( s)  u 0 ( s)  k u u ( s ) Fb ( s)  u( s)  u( s  1)  k

(17)

u u ( s )

Thus the total bearing force at time step s is expressed as: Fb ( s)  u( s)  u0 ( s)  k uu ( s )  u( s)  u( s  1) k

The unbalance force terms

Fex  s   ms2e  s  sin   s 

u u ( s )

(18)

and Fey  s   ms2e  s  cos   s  are

presented by the system dynamic parameters e and γ as:

  s  =arctan

u y  s     cos(s ts) ux  s     sin(s ts)

e  s   ux  s   u y  s    2  2 ux  s  sin(s ts)  u y  s  cos(s ts)  2

2

(19)

(20)

For the time step s the unbalance force vector Fe, which is not only influenced by the static eccentricity ε and planar dynamic displacements ux and uy, is obtained by:

Fe  s   [ms2e  s  sin   s  ms2e  s  cos   s  0 0 0]

(21)

For this work, the fourth order Runge-Kutta algorithm is selected for the numerical integration. The process of computing the values of space vectors x at time step s through the algorithm of 4th-order

Runge-Kutta is showed as: x  s  1  x  s  

1 g1  s   2g 2  s   2g3  s   g 4  s   6

(22)

Where:

g1  s   tf  s  1   g 2  s   tf  x  s   g1  s   2   1   g 3  s   tf  x  s   g 2  s   2  

(23)

g 4  s   tf  x  s   g 3  s   The total response time tend  send t is necessary to be long enough to achieve the adequate steady state response. Additionally, the error for the numerical approximation of x  send  caused by time discretization also matters the convergence of the method. The convergence for the response can be expressed in Eq.(24): (24) lim x  send   x(tend ) t 0 s  send

From [29] the method that can be denoted as convergence is required to converge in a relative large range of all reasonable initial values. The initial vector for each step is needed which can be 0

1

s

2

written as x0 , x0 , x0 ,...., x0end

1

for this method. Thus, the initial values can be decided by Eq.(25): lim x0 (t )  x0 ,   0,1, 2,..., send  1

(25)

t 0

For the problem described in this work, which is existing in the practical engineering, the quasi-static results xq0 at the corresponding external loads, which means the velocity vector u is a null vector in mechanical system, are chosen as the initial values, which is showing as Eq.(26). xq 0

W  x  k x

Wy ky

Wz kz

Mx k x

My k y

 0 0 0 0 0 

T

(26)

The calculation flow chart for the vibration response of the rolling-element bearing is shown in Fig.4.

Fig.4 Calculation flow chart

5. Case Study One type of cylindrical roller bearing will be introduced in this chapter and the corresponding parameters are showing in Table 1. The mass of the rotor is set as m=2.14 kg. The internal clearance is chosen based on ISO5753-1:2009 [30] and set as 60 μm. For the cylindrical roller bearings without crown rollers, the external load is only applied on the Y direction and the plane motions are the main focus in this section for model simplification and less computation time consumption. Thus, the effects from the unbalance force and the internal clearances on the system dynamic plane motions are parametrically studied. Table 1. Cylindrical Roller Bearing Parameters (Based on SKF NU306 ECJ) Parameters

Values (Unit :mm)

Number of Rollers

11

Roller Length

13

Roller Diameter

11

Pitch Diameter

51.5

Bearing Width

19

Outer Diameter

72

Inner Diameter of Outer Raceway

62.5

Outer Diameter of Inner Raceway

40.5

Bore Diameter

30

The load deflection distributions of this kind of cylindrical roller bearing are captured through FE/CM mixed method in this work. Also, it can be obtained through other effective analytical or experimental techniques whichever are accurate enough to catch the nonlinear properties for the

REBs in static state.

5.1 Stability analysis on X direction Based on the description above, even the external load only applied on Y direction, the system also is disturbed on X direction caused by the asymmetric internal configuration inside of the bearing over ball pass cycles. The equation of the vibration response on X direction is recalled in Eq.(27).

mux  Fcx  Fbx  Fex

(27)

Which is an autonomous system can be rewritten in: ux  

1 1 Fcx (ux )   Fbx (ux )  Fex (u x ) m m

The functions h and g are denoting the velocity force and displacement force: 1 h(u x )  Fcx (u x ) m 1 g (u x )   Fbx (u x )  Fex (u x ) m

(28)

(29)

Hence for the dynamic process, the potential energy on X direction is presented by Eq.(30).

   g (ux )dux

(30)

The kinematic energy of the rotor on X direction is defined as: 

1 2 ux 2

(31)

The total energy only on X direction is denoted by: E       g (u x )du x 

1 2 ux 2

(32)

We can obtain the rate of change of energy E by: dE  g (u x )u x  u x u x dt  g (u x )u x   h(u x )  g (u x )  u x   h(u x )u x 

(33)

ux Fcx  u x  m

In this case, the Raleigh damping constants are chosen as α=10-2 and β=10-6; since the mean bearing stiffness on X direction is not active when no load applied on X direction for the REB with internal clearance, the average effective bearing stiffness inside of the relationship with damping force Fcx is the diagonal term relative to the X direction from the alternating bearing stiffness matrix without the tiny off-diagonal terms. Thus the term dE/dt can be rewritten in:

dE     ux2    k xx  dt m  

(34)

The alternating stiffness term k xx is keeping positive along the ball pass cycle and the rate of energy change term will keep negative. The total energy E is decreasing along every path of (ux , ux ) , and then the phase path of X direction goes to the equilibrium points. 10

 =0

Displacement u x / μm

0 -10

0

300

10

600

 =5 μm

0 -10

0

300

10

600

 =10 μm

0 -10

0

300

600

Ball pass cycle Fig.5 Vibration displacement on X direction, Wy=3 kN, ns=12000 rpm

Fig.5 shows the time history of the vibration displacement ux for different values of static eccentricity ε. For ε=0, the total energy is descending with the effect of viscous damping until the initial energy uses up, reaching to the equilibrium state which is (ux , ux )=  0,0  . For the cases with ε>0, when the energy from the initial position energy runs up and the energy from unbalance force is vibrating with time, the vibration behaviors will come to the convergent and steady state. The FFT frequency spectrum for the X direction vibration displacement with ε=10 μm is showed in Fig. 6. The major peaks of vibration response are observed at the rotation frequency fs and the combination of the ball pass frequency with the rotation frequency fb-fs, fb+fs, 2fb-fs and 2fb-fs,The largest peak amplitude occurs at fb+fs.

1.4

Amplitude / m

1.2

fb  f s

1 0.8 0.6 0.4 0.2

fs

fb  f s 2 fb  f s 2 fb  f s

0

1000

2000

3000

4000

5000

6000

Frequency /Hz Fig.6 FFT frequency spectrum for X vibration displacement, Wy=3 kN, ns=12000 rpm, ε=10 μm

5.2 Effect from speed fluctuation Peak-peak value of displacement u y / μm

100

Time-varying bearing stiffness model 50

0 500 20

1000

1500

2000

New bearing force model

2500

B

3000

3500

4000

Wy  1 kN Wy  3 kN Wy  5 kN

C

15

A 10 5 0 500

1000

1500

2000

2500

3000

3500

4000

Frequency /Hz Fig.7 Speed sweep plot for Y direction, ε=5 μm

The time-varying bearing stiffness model, which is widely applied to transmission systems with rolling bearings [9,21,31], is introduced in this section for the comparison with the new model built in this work. The bearing force Fb is decided by the bearing time-varying stiffness matrix Kˆ b . Neglecting the off-diagonal term, time-varying stiffness matrix is expressed as: ˆ  diag kˆ K b  xx

kˆyy

kˆzz

kˆ x x

kˆ y y 

(35)

The time varying bearing stiffness terms kˆxx , kˆyy , kˆzz , kˆ  , kˆ  inside of Kˆ b are calculated as the same x x

y y

way as the alternating stiffness showed in Eq.(9), but independent from the influence of the external loads. Thus, the bearing force is written as: ˆ u Fb   K b

(36)

In which  is the gap function for judging whether the rolling elements are in contact with raceways;  0  1



u  u0 u  u0

(37)

u 0 is the no-load displacement and also shows in Section 3.2.

Fig.7 shows the peak-peak value of displacement on Y direction with speed fluctuation from two different bearing force models. The rated rotation speed is 11000 rpm for this type of bearing, while the rotation speed range is from 8000 to 50000 rpm corresponding to the excitation frequency range from 576.7 Hz to 3604 Hz for science research purpose. For the initial speed, the initial points for the integration process are decided based on Eq.(26). For the subsequent speed condition, the steady state solutions will be the initial points for the previous speed circumstance. It is seen clearly from Fig.6 that the vibration peaks occurs at the natural frequency of Y direction zones and corresponding sub-harmonic areas. For the time-varying bearing stiffness model, the results from the speed-up sweep is very similar to the speed-down sweep results, thus only the vibration displacements results for the speed-up sweep process is presented by the dash lines in upper plot. The peak frequencies keep the same when the external load increases from 1kN to 5 kN and the peak-peak value of displacement of Y direction rises clearly with increasing load. But for the new bearing force model in this work, with the increasing loads, the peaks shift to larger frequency area. When the external load Wy=3 kN, the 1/2 subharmonic frequency component excites the largest peak-peak displacement value amplitude A. For the Wy=1 kN and 5 kN, the major peaks occur at the Y direction natural frequencies with the phenomena of jump amplitude. In addition, the bending-to-left peaks at zones B and C are obvious, which shows the strong soften nonlinear properties.

5.3 Effect from unbalance force Fig.8 shows the unbalance force varies visibly with rotor period cycle Ts both for different external loads. For the certain static eccentricity value and rotation speed, the unbalance force slightly increases when external load is largely enhancing, while the rotation speed is the main factor for the magnitude of unbalance force. The comparison for the results of bearing forces with and without unbalance force is plotted in Fig.9. The mean values of bearing force of models with unbalance force are larger than the results from models without unbalance force. The shapes of

bearing force vibrations are more complicated when considering unbalance force, while the magnitude of the rate of forces fluctuation are not sensitive to the unbalance force. 160 150

Unbalance force Fe / N

140 130 120 110 100 90 Ts

80 5

Wy =0.2kN

Ball pass cycle

/ 102

Wy =3kN

Wy =5kN

6

Fig.8 Unbalance force with varied external loads, ε=5 μm, ωs=11000 rpm with unbalance force without unbalance force

600 400

Wy  0.2 kN

Bearing force Fby / N

200 0 500 4000

515

Wy  3 kN

3000 2000 500 5200

515

5100

530

Wy  5 kN

5000 4900 500

530

515

530

Ball pass cycle Fig.9 Effect of unbalance force on bearing force vibration, ε=5 μm, ωs=11000 rpm

5.4 Effect from external load The vibration behaviors for the rotor-bearing system under rated speed ωs=11000 rpm are studied in this section. With larger load, the RMS value of displacement escalates while the fluctuation rules of the displacement are not related to the magnitude of loads. For the operation condition showed in Fig.10, the largest fluctuation is obtained by the case with Wy=3 kN, which reaches to 3.15 μm of peak-peak value. From the FFT frequency spectrums, for the low load occasion, the effect of rotation frequency of rotor is obvious on the fluctuation amplitude. As showing in the case of Wy=0.2 kN, rotating frequency of rotor fs, ball pass frequency fb and its multiples 2fb,3fb and the

combinations frequencies fb+fs, 2fb+fs dominate the displacement fluctuation amplitude. When the load grows, resonant peaks are observed at ball pass frequency fb and multiples 2fb, 3fb and 4fb for

32.5

fb

0.4

32

2 fb

2 fb  f s

0.2

31.5

3 fb

fs

31 3

5

7

40

Velocity u y / mm  s 1

Amplitude u y / μm

Displacement u y / μm

Wy=3 kN and for Wy=5 kN with additional obvious peak at rotation frequency fs.

0 0

fb  f s 2

1.5

4

10 0

Wy  0.2 kN

-10 31

32

33

50

3 fb

1 38

0 0.5

36 3

5

7

fb

2 fb

0 0

42

0.2

41.5

0.1

2

5

7

4

-50 36

fb

38

Wy  5 kN

0

0 0

3 fb

4 fb

2

4

-2

41.6

41.8

42

Displacement u y / μm

Frequency /103 Hz

Ball pass cycle / 102

40

2

2 fb

fs

41 3

Wy  3 kN

4 fb

500

0

5

7

4000 3000 2000

5

7

Frequency components amplitude / N

Bearing force Fby / N

Fig.10 Vibration behaviors of displacements with varied external loads, ε=5 μm, ωs=11000 rpm

2 fb

100

3 fb

50

fs 0 0

2 fb  f s 3

1000

3 fb

Wy  3 kN 2 fb

4 fb

0 0 40

5200

20 f s

3

3 fb

5

Ball pass cycle / 102

7

6

500

5400

5000

Wy  0.2 kN

fb

0 0

fb

2 fb

6

Wy  5 kN

4 fb 5 fb 6 fb

3 f b  f s3 Frequency /103 Hz

6

Fig.11 Vibration behaviors of bearing forces with varied external loads, ε=5 μm, ωs=11000 rpm

Fig.11 shows the vibration situation of bearing force Fby over the steady state. The larger external

loads lead to larger RMS value of bearing forces but have no obvious relationship with the peak-peak values of the bearing forces. As showing in Fig.10, the bearing force is fluctuating from around 2340 N to the values around 3930 N when Wy=3 kN, while the triple of ball pass frequency 3fb masters the force fluctuation amplitude. The same conclusions as the response of fluctuation of displacements that the response dominated by rotation frequency are more obvious with low load, also come to the force response, which is observed by the peaks at fs and 2fb+ fs for Wy=0.2 kN. For the high-load condition, the components of higher multiples of ball pass frequency are more noteworthy. While for Wy=5 kN, the frequency components of 6fb still plays an irreducible role in the fluctuation amplitude response of bearing forces.

6 Conclusion This work proposed an advanced high-efficient model for the rotor-bearing system dynamic behaviors. This model can well preserve and take advantages of the full 3-dimensional load-deflection distributions of rolling element bearings, which consider the bearing varying compliance, body compliance, radial internal clearance and the strong non-linear characteristics attached to varied external loads. Two different types of bearing stiffness, which are the mean stiffness and alternating stiffness, are classified physically and defined clearly. A new computation model including these two different bearing stiffness model of rolling element bearing-rotor system for dynamic behaviors is proposed. One kind of cylindrical roller bearing is introduced for case study and the dynamic behaviors of one direction applied load are studied deeply. Intricate and highly strong nonlinear vibration behaviors are observed: 1) For the system only with Y direction external load, the nonlinear behavior are observed on X direction because of the asymmetric configuration over ball pass cycles. While there is internal clearance, the bearing doesn’t transfer loads on X direction. Only the alternating stiffness term k xx is activated, which is depended on the value of bearing force on Y direction for the static model build upon FE/CM method. Based on this physical model, the X vibration is dominated by the ball pass frequency, rotating frequency and their combinations, but always goes to a static state with the effect of viscous damping. 2) The peak resonances are obtained when the exciting frequency passing the zone of natural frequency of Y direction and its sub-harmonic areas. The new model built in this work can capture the nonlinear dynamic behaviors caused by the change of the instantaneous equilibrium position and the disturbance displacements. The results show that resonant areas shift to the larger frequency when increasing external loads while the situation of vibration fluctuation is less related to the magnitude of external load. Remarkable jump phenomenon and soften bend-to-left peak properties are observed at certain loads. 3) The unbalance force is not sensitive to the variation of external loads, but firmly depended on the

rotation speeds. The mean values of vibration characteristics are decreasing when the unbalance force is irrespective. 4) For the lower loads, the resonances are influenced by rotation frequency more obviously. For the condition with larger loads, the multiples of ball pass frequency play more conspicuous role in the frequency spectrum of vibration characteristics.

Nomenclature O  XY Generalized coordinate O'  XY Coordinate attached to rotor m Rotor mass

I b Moment of inertia of the rotor I p Polar moment of inertia of rotor

 ,  Rayleigh damping constants

s Rotor angular speed ns Rotor rotational speed Wx ,Wy ,Wz External loads

M x , M y External torques Fex , Fey Bearing force on X-Y plane Fex , Fey Unbalance forces

Fcx Damping force on X direction e Dynamic eccentricity

 Static eccentricity

 Dynamic phase angle ux , u y , uz , x , y Dynamic displacements uˆx , uˆ y , uˆz ,ˆx ,ˆy Effective elastic displacements

k x , k x , k z , k x , k y Bearing mean stiffness on different directions k xx , k yy , k zz , k x x , k y y Diagonal terms in alternating stiffness matrix

k xy , k yz , k x x , k x y , k yx , k yz , k y x , k y y , k zx , k zy , k z x , k z y , k x x , k x y , k x z , k x y , k y x , k y y , k y z , k y x

bc Bearing internal clearance h(ux ) Velocity force g (u x ) Displacement force

 Potential energy

 Kinematic energy

Off-diagonal terms in alternating stiffness matrix

E Total energy

f b Ball pass frequency f s Rotor rotating frequency Ts Rotor period cycle

t Time t Small time increment

s Time step

 Small time interval in remainder u, u, u Displacement, velocity and acceleration vectors M Mass matrix

C Damping matrix G Gyroscopic matrix Fb Bearing support force vector Fe Unbalance force vector

W External loads vector K Average effective bearing stiffness matrix

u e Displacement at a certain instantaneous equilibrium point

Δu Small disturbance displacement around equilibrium point Fb Bearing disturbance force caused by disturbance displacement Fb Bearing support force at a certain instantaneous equilibrium point

u 0 No-load displacement

x, x State variables x q 0 Initial values for state equations X State equations

f Nonlinear functions between acceleration vector with displacement and velocity vectors

Acknowledgments This work is financially supported by Chongqing Science and Technology Project (Grant No.10262235320180019), National Natural Science Foundation of China (Grant No. 51575061) and China Scholarship Council. The authors also gratefully acknowledge for the help from Professor Robert Parker’s group in Virginia Tech and Dr. Vijayakar Sandeep for the FE/CM modelling.

Reference: [1] [2]

Cong, F., et al., Vibration model of rolling element bearings in a rotor-bearing system for fault diagnosis. Journal of Sound and Vibration, 2013. 332(8): p. 2081-2097. El-Thalji, I. and E. Jantunen, A summary of fault modelling and predictive health monitoring of rolling element bearings. Mechanical Systems and Signal Processing, 2015. 60-61: p.

[3]

[4] [5] [6] [7] [8] [9] [10]

[11]

[12]

[13] [14]

[15]

[16] [17]

[18]

[19]

252-272. Zhang, M., Z. Jiang, and K. Feng, Research on variational mode decomposition in rolling bearings fault diagnosis of the multistage centrifugal pump. Mechanical Systems and Signal Processing, 2017. 93: p. 460-493. Sharma, A., et al., Nonlinear dynamic investigations on rolling element bearings: A review. 2018. 10(3): p. 1687814018764148. Sunnersjö, C., Varying compliance vibrations of rolling bearings. Journal of sound and vibration, 1978. 58(3): p. 363-373. While, M., Rolling element bearing vibration transfer characteristics: effect of stiffness. Journal of applied mechanics, 1979. 46(3): p. 677-684. Fukata, S., et al., On the Radial Vibration of Ball Bearings : Computer Simulation. Bulletin of the Jsme, 1985. 28(239): p. 899-904. Lim, T.C. and R. Singh, Vibration transmission through rolling element bearings, part I: Bearing stiffness formulation. Journal of Sound & Vibration, 1990. 139(2): p. 179-199. Lim, T. and R. Singh, Vibration transmission through rolling element bearings, part II: system studies. Journal of Sound and Vibration, 1990. 139(2): p. 201-225. Tiwari, M., K. Gupta, and O. Prakash, Effect of radial internal clearance of a ball bearing on the dynamics of a balanced horizontal rotor. Journal of Sound and vibration, 2000. 238(5): p. 723-756. Tiwari, M., K. Gupta, and O. Prakash, DYNAMIC RESPONSE OF AN UNBALANCED ROTOR SUPPORTED ON BALL BEARINGS. Journal of Sound & Vibration, 2000. 238(5): p. 757-779. Upadhyay, S.H., S.P. Harsha, and S.C. Jain, Analysis of Nonlinear Phenomena in High Speed Ball Bearings due to Radial Clearance and Unbalanced Rotor Effects. Journal of Vibration & Control, 2010. 16(1): p. 65-88. Mevel, B. and J.L. Guyader, Experiments on routes to chaos in ball bearings. Journal of Sound & Vibration, 2008. 318(3): p. 549-564. Kappaganthu, K. and C. Nataraj, Nonlinear modeling and analysis of a rolling element bearing with a clearance. Communications in Nonlinear Science and Numerical Simulation, 2011. 16(10): p. 4134-4145. Zhang, Z., Y. Chen, and Q. Cao, Bifurcations and hysteresis of varying compliance vibrations in the primary parametric resonance for a ball bearing. Journal of Sound & Vibration, 2015. 350: p. 171-184. Li, Y., et al., A General Method for the Dynamic Modeling of Ball Bearing–Rotor Systems. Journal of Manufacturing Science and Engineering, 2015. 137(2): p. 021016. El-Saeidy, F.M.A. and F. Sticher, Dynamics of a Rigid Rotor Linear/Nonlinear Bearings System Subject to Rotating Unbalance and Base Excitations. Journal of Vibration & Control, 2010. 16(3): p. 403-438. Zhang, X., et al., Stability analysis of a rotor–bearing system with time-varying bearing stiffness due to finite number of balls and unbalanced force. Journal of Sound and Vibration, 2013. 332(25): p. 6768-6784. Sharma, A., M. Amarnath, and P. Kankar. Effect of unbalanced rotor on the dynamics of cylindrical roller bearings. in Proceedings of the 9th IFToMM International Conference on Rotor Dynamics. 2015. Springer.

[20]

[21] [22]

[23] [24] [25]

[26] [27] [28] [29] [30] [31]

Kurvinen, E., J. Sopanen, and A. Mikkola, Ball bearing model performance on various sized rotors with and without centrifugal and gyroscopic forces. Mechanism and Machine Theory, 2015. 90: p. 240-260. Sheng, X., et al., Calculation of ball bearing speed-varying stiffness. Mechanism and Machine Theory, 2014. 81: p. 166-180. Gunduz, A. and R. Singh, Stiffness matrix formulation for double row angular contact ball bearings: Analytical development and validation. Journal of Sound and Vibration, 2013. 332(22): p. 5898-5916. Gargiulo, E., A simple way to estimate bearing stiffness. Machine Design, 1980. 52(17): p. 107-110. Knaapen, R., L. Kodde, and d. Kraker, Experimental determination of rolling element bearing stiffness. 1997: Technische Universiteit Eindhoven. Gunduz, A., Multi-dimensional stiffness characteristics of double row angular contact ball bearings and their role in influencing vibration modes. Dissertations & Theses - Gradworks, 2012. Li, H.L., et al., Periodic response analysis of a misaligned rotor system by harmonic balance method with alternating frequency/time domain technique. 2016. 59(11): p. 1717. Fu, C., et al., Steady-state response analysis of cracked rotors with uncertain‑ but‑ bounded parameters using a polynomial surrogate method. 2019. 68: p. 240-256. Lu, K., et al., Review for order reduction based on proper orthogonal decomposition and outlooks of applications in mechanical systems. 2019. 123: p. 264-297. Deuflhard, P. and F. Bornemann, Scientific Computing with Ordinary Differential Equations. 2002. International Standard ISO 5753-1:2009, Rolling bearings - Internal clearance - Part 1: Radial internal clearance for radial bearings. International Standards Organizations, 2009. Liu, G., J. Hong, and R.G. Parker, Influence of simultaneous time-varying bearing and tooth mesh stiffness fluctuations on spur gear pair vibration. Nonlinear Dynamics, 2019: p. 1-22.

Declaration of interests 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.