René E El

Melhem Bat

Blaise Pascal

Insa : R Gourdon

M Sté

Daniele

M Gé

Scorletti
email: gerard.scorletti@ec-lyon.fr

Isabelle Mme

Von Bueltzingloewen

Zhu S C | Ré

Keywords: Elasto-hydrodynamic Lubrication (EHL), contact model in bearing matrices, bearing-rotor model, alternative bearing bench Elasto-hydrodynamic Lubrication (EHL), transient EHL, oil film stiffness and damping distributions Chapter Direct bearing model, roller bearing, slicing technique, interaction among slices, force and moment distributions, bearing stiffness Bearing test bench, inertial loading, topological optimization, electromagnetic vibrator, sliding mode control

Contribution to the modeling of rolling element bearing

Acknowledgment/Remerciements

I would like to express to Mr. Zhang my very deep gratitude to have provided me with advanced projects of National level, financial supports and technical advice, especially of brilliant ideas in train techniques, and for the most exciting encouragement he gave and for the constant interest that he showed during these five years.

I am extremely sensitive to the honor that Mr. Nelias has allowed me to follow him for over two years to study bearing techniques, completely grateful for technical documents he gave, very pleased by the nice discussions we had and always thankful to his constant kindness. I could never forget what I experienced these years in Lyon, thanks to him. I always feel indebted to my parents and other members in my family. I spent less time with them, until my grandfather (maternal) passed away 2 months ago, I suddenly find something in my life lost and some space in my heart empty. Now, I really miss him.

I would like to thank particularly Professor Shuying Gao and her husband-Professor Luming Wu, who passed away on the same day when I returned from Lyon back to China one year ago. I cannot forget him forever; he firstly gave me the hope in bearing research and gave me great help in my road abroad administratively.

The work for a Ph.D. candidate is to find a road in the maze. Even although a lot of exits can be in fine found, a lot of blind alleys can also be possibly in the front. It sometimes needs immediate retreats and timely helps from colleagues and friends. I would like to thank Professor Dongli Song, Professor Guiming Mei, and my colleagues Q. Xiong, L. Wei, H.L. Wu for their helps in China, and a lot of other colleagues, Y.C. Zhang, S. Chen, W.J. Gao, Eric, Damien, etc. for their support in France.

I would also like to thank my Chinese friend, W. Gong, headman of hardware test department in Huawei Company, Chengdu, who gave me a lot of suggestions in programming

Introduction

From contact model to bearing model, and from bearing-rotor model to bench reproduction technique, the four chapters deliver the author's understandings in bearing engineering application.

Chapter 1:

A new method is proposed to solve transient elasto-hydrodynamic lubrication (EHL) problem. First, original dynamic EHL behavior is estimated by steady-state EHL solution, where elastic deformation theory is specifically combined with oil film stiffness distribution equation instead of steady-state Reynolds equation. Second, subsequent dynamic EHL procedure develops from recursively using transient distributed oil film stiffness and damping, where each time-marching solution is iteratively searched by ensuring both oil film force Introduction | Zhu S.C.

Chapter 3:

A uniform derivation for overall stiffness and damping matrices of both tapered and ball bearings is conducted, using the pioneer work of Luc-Houpert (1997), with slice technique considered in distributed roller load calculation and varying contact angle considered in fluctuating ball transmission load calculation. Simplified EHL and Hertz relations are linearly included in contact force and in fine analytically facilitate formulated matrices. Nonlinear bearing-rotor system is modeled with a better linear mathematical equation and solved in simpler enhanced numerical manner of Adams. Meanwhile, nonlinear bearing-rotor system is also remodeled and resolved with analytical method of HB-AFT. It performs a faster efficiency and an easier system behavior-identifying capacity, especially when considering 3D condition. Besides modeling and solving nonlinear bearing-rotor system, some factors participated can also be qualitatively and quantitatively considered, including initial states, damping, and eccentricity, internal and external excitations.

Keywords:

Bearing damping and stiffness matrices, ball and roller bearings, slice technique, varying contact angle, bearing-rotor system, enhanced Adams method, HB-AFT.

Chapter 4:

This chapter shows a new design of bearing test bench. Two ideas are presented herein: 1 functionally, it reproduces dynamic loads acting on rotor by considering internal and external excitations. Internal excitations are applied, as most bearing manufacturers manipulate, in the manner of motor drive, while external excitations are specifically loaded herein by generating housing motions and corresponding rotor inertial loads, which requires programming an electromagnetic vibrator which supports entire system. It appeals to the idea the authors have ever delivered in work about nonlinear bearing-rotor system dynamics. Meanwhile, an alternative idea of uniquely reproducing bearing contact load is also performed for bearing study. 2 financially, to reduce costs for the electricity used in bearing tests (e.g. for train products) conducted in our first-generation bench, a second-generation bench is also proposed to increase usage of vibrator, where a sliding mode control idea is also used in optimized design. With the construction of bearing test platform, a series of research ideas on high-speed, over-loading and vibration conditions can be explored.

Pré face | Zhu S.C.

Préface

La technologie des roulements est un savoir connu, les questions traditionnelles sur les paramè tres interfaciaux dynamiques sont toujours populaires à la fois dans les roulements et é galement dans les applications adé quates de transmission, puisque les analyses de couplage pré cises et efficientes sur les relations dynamiques diffé rentes et les conceptions subsé quentes de la composition structurelle et des paramè tres sont toujours requis en fonction des exigences diverses. L'interfé rence des roulements caracté risé e par les matrices de coefficients et son association avec le systè me des rotors sont de primeure inté rê t à cette thè se.

Davantage de connaissance adéquate sur le contact à l'intérieur d'un roulement et son approche beaucoup plus simple d'interpréter les charges de contact du roulement sont très importants pour comprendre et dé crire de maniè re correspondante les erreurs que les paramè tres simplifié s de contact ont autrefois introduit dans les matrices de rigidité et d'amortissement des roulements. Inversement, comment interpréter les paramètres distribués de contact de roulement avec les paramè tres assemblé s est é galement une question à long terme d'intérêt et d'importance, puisque les paramètres assemblés de contact sont plus facilement impliqué s dans les matrices de transmission des roulements.

Cette thè se explique dans

Cette thè se explique é galement dans le chapitre 3:

Comment dé velopper des matrices de roulement plus pratiques et des mé thodes plus fiables de ré soudre des modè les de roulements-rotor non liné aires de maniè re uniforme et é conomique?

Si les mé thodes de fré quence-domaine repré sentent une meilleure capacité d'identification par assistance informatique plutôt que les mé thodes de temps-domaine ? L'analyse des domaines temps/fréquence sur le système dynamique du rotor présenté au chapitre 3 donne soit qualitativement soit quantitativement des informations utiles caché es entre le roulement et le rotor. Quelques vé rifications par des tests correspondants sont requis pour reconnaî tre les paramè tres des roulements ou les motions du rotor, de mê me que pour estimer les conditions de fatigue du roulement. Une premiè re gé né ration de banc d'essai de roulement proposé reproduit les comportements dynamiques du rotor avec diffé rentes caracté ristiques de roulements é mergeant sous diffé rentes tensions. Cependant, une vibration intense couplé e à une charge importante de contact imposé e de maniè re fluctuante à l'intérieur d'un roulement ne peut pas être reproduit facilement, puisque la plupart des vibreurs ne sont pas capables de faire vibrer intensément les bancs d'essais géants des roulements.

Cette thè se pré sente enfin dans le chapitre 4:

Comment construire un banc d'essai alternatif des roulements pour reproduire les comportements de contact à l'intérieur d'un roulement sous une charge importante et une vibration à haute fré quence.

Ré sumé | Zhu S.C.

Résumé

Quelques dé veloppements sont dé montré s dans cette thè se, du modè le de contact au modè le de roulement et du modè le de roulement-rotor à la technique de reproduction par bancs.

Tout d'abord, la mé thode de perturbation est é tendue afin de prouver que la rigidité et l'amortissement distribués en état stable EHL ne remplissent pas les conditions transitoires et que les coefficients inté gré s ne sont pas approprié s pour les analyses dynamiques。 La loi de distribution de la rigidité et l'amortissement du film d'huile de l'EHL est proposée. En combinant l'état initial d'une solution stable avec l'amortissement et la rigidité distribués de l'EHL transitoire, les comportements du film d'huile peuvent être prédits.

Par la suite, l'effet du contact en bordure entre le cylindre et le courant est pris en compte dans les transmissions de roulement. Les méthodes d'approximation première et seconde sont proposées pour interpréter numériquement cet effet. La méthode d'approximation première se concentre sur l'adoption de la méthode du coefficient d'influence pour reproduire l'effet de contact de même que l'utilisation du procédé de normalisation pour respecter les résultats empiriques ; alors que la méthode d'approximation seconde se concentre sur une maniè re analytique rapide de dé crire la pression Hertz au contact entre le cylindre et le courant.

Ensuite, les matrices de roulements pratiques sont utilisé es pour dé velopper un code de dynamiques 3D du roulement-rotor, basé sur des formules Hertz et EHL classiques. En tempsdomaine, la mé thode Adams quatre-ordres est amé lioré e pour ré soudre les é quations incré menté es des dynamiques non-liné aires. En fré quence-domaine, la mé thode analytique de HB-AFT est particuliè rement employé e pour combler des matrices de roulement en 3D et des solutions d'état stable 3D. De plus, la mé thode de continuation arc-longueur augmente la vitesse de calcul. En utilisant les deux mé thodes, certains facteurs peuvent ê tre considé ré s qualitativement et quantitativement, en incluant les états initiaux, l'amortissement, l'excentricité et les tensions internes et externes.

Introduction | Zhu S.C.

Introduction

A partir du modè le de contact au modè le de roulement, et à partir du modè le de roulement-rotor à la technique de reproduction par bancs, les quatre chapitres dé livrent notre compréhension dans l'application de l'ingénierie des roulements.

Chapitre 1:

Une nouvelle mé thode est proposé e pour ré soudre le problè me de lubrification transitoire é lasto-hydrodynamique (EHL). Tout d'abord, le comportement de dynamique EHL commence par une solution de EHL à l'état stable, où la thé orie de la dé formation é lastique est spé cifiquement combiné e à l'équation de la rigidité distribuée du film d'huile au lieu de l'équation de Reynolds à l'état stable. Ensuite, le comportement subsé quent de dynamique EHL se dé veloppe en utilisant ré cursivement la rigidité et l'amortissement distribués du film d'huile transitoire, où chaque solution transitoire est ité rativement recherché e par l'assurance que le chargement externe par échelon s'associe à la fois avec la croissance de la force du film d'huile ainsi que la ré vision de la dé formation é lastique. In fine, cette mé thode, non seulement augmente la vitesse de calcul, mais aussi dé montre la disponibilité de la rigidité et de l'amortissement distribué s de l'EHL en acqué rant des courbes idé ales de EHL transitoire. La mé thode rapide d'analyse du contact de film dynamique sur les roulements, les engrenages, etc. pourrait convenir techniquement.

Mots clé s:

Lubrification Elasto-Hydrodynamique (EHL), EHL transitoire, distribution de la rigidité et de l'amortissement du film d'huile.

Chapitre 2:

Ce chapitre renforce la technique partielle classique pour calculer les forces et les moments agissant sur les roulements ronds avec un modè le direct de roulement. Tout d'abord, les calculs de dé formation et de charge empirique classique sont revus.

Ensuite, les modè les d'approximation premier et second du contact du roulement sont dé veloppé s de faç on polyvalente pour prendre en compte la concentration de la pression dans le bord du roulement. Enfin, en considé rant la capacité lors de l'application, un modè le adé quat d'ingénierie est proposé , en utilisant le process de normalisation. Les Cela appuie l'idée des auteurs qui ont remis le rapport sur la dynamique du système des roulements-rotor non liné aires. Cependant, une idé e alternative en reproduisant uniquement la charge du contact du roulement est également présentée pour l'étude des roulements. 2financiè rement. Afin de réduire les coûts d'électricité utilisés dans les tests de roulements (par exemple les produits ferroviaires) mené s dans le banc de premiè re gé né ration, un banc de seconde gé né ration est proposé pour augmenter l'usage du vibrateur, où une idé e de contrôle en mode glissant est é galement utilisé e en conception optimisé e. Avec la construction d'une plateforme de banc d'essais, une sé rie d'idé es de recherches sur des conditions à grande vitesse, surcharge et vibrations peuvent ê tre exploré es.

Mots clé s : banc d'essais de roulements ; charge inertielle ; optimisation topologique ; vibrateur é lectromagné tique ; contrôle en mode glissant.

Chapter 1

Chapter 1 Calculation on Elasto-hydrodynamic Lubrication 1.1 Introduction

EHL contact qualities in transmitting components (bearing, gear, etc.) and in special dynamic systems with open connections (wheel-rail contact, etc.), can be measured by stiffness and damping parameters. In dynamic design, integrated (or lumped) parameters are often preferred instead of distributed values due to their feasibility in multi-body dynamic system, and it must be also recognized that steady-state parameters are used exclusively since transient values are technically hard to seize.

In the past decades, strenuous efforts have gone into improving the accuracy and efficiency of steady state EHL calculation, direct iterative method for light load in point contact was introduced by Hamrock and Dowson [START_REF] Hamrock | Isothermal EHD Lubrication of Points Contacts: Part I-Theoretical Formulation[END_REF], inverse method was also developed for heavy load, first employed in line contact by Dowson and Higginson [2], and later in point contact by Evans and Snidle [3]. Since Newton Raphson method was used, relevant researches in EHL problem have become more popular [4], [START_REF] Oh | Numerical solution of the point contact problem using the finite element method[END_REF], [START_REF] Houpert | Fast Approach for calculating Film Thicknesses and Pressures in Elasto-hydrodynamically Lubricated Contacts at High Loads[END_REF]. When multigrid method was introduced by Lubrecht [7], incredibly high efficiency in calculation was obtained.

Interest in numerical technique also led to relevant developments of derived methods, such as multilevel multi-integration method [8], FFT method [9], etc.

No matter for what method, it is still difficult to apply it to dynamic systems, from transient EHL problem to subsequent application in vibration system.

Apart from full numerical methods, a simplified method neglecting the Poiseuille term used by L.Chang [10], S. Messé and Lubrecht,A.A. [START_REF] Messe | Approximating EHL Film Thickness Profiles Under Transient Conditions‖[END_REF], and others has brought some efficiency in analysis on transient EHL problems.

As for integrated EHL behaviors, other efforts were made. A rough and direct integration of Reynolds equation was proposed and often used in parallel to Hertz relation in bearing dynamic calculation [START_REF] Walford | The sources of damping in rolling element bearings under oscillating conditions‖[END_REF]~ [START_REF] Nonato | On the non-linear dynamic behaviour of elastohydrodynamic lubricated point contact[END_REF]. Recently, another rough and direct integration method was also used, oil film parameters are integrated within different spaces in [START_REF] Wiegert | A simplified elasto-hydrodynamic contact model capturing the nonlinear vibration behaviour[END_REF], [START_REF] Wiegert | An analytical expression of the normal force of hydrodynamic line contacts under transient conditions[END_REF], but in series with Hertz stiffness. Both numerical method and curve fitting method were used to extract lumped parameters from distributed EHL parameters by Shriniwas P. Chippa and Mihir 2017.09. [START_REF] Wiegert | An analytical expression of the normal force of hydrodynamic line contacts under transient conditions[END_REF] Sarangi [START_REF] Sarangi | Stiffness and damping characteristics of lubricated ball bearings. Part 1: theoretical formulation‖[END_REF], [START_REF] Chippa | On the Dynamics of Lubricated Cylindrical Roller Bearings, Part I: Evaluation of Stiffness and Damping Characteristics‖[END_REF], [START_REF] Chippa | On the Dynamics of Lubricated Cylindrical Roller Bearings, Part ii: linear and nonlinear vibration analysis‖[END_REF], and until now, this could be the most complex way to get EHL damping and stiffness coefficients. Before getting transient lumped parameters, transient distributed parameters should be studied. This chapter is oriented towards the study of transient distributed EHL behaviors.

Besides the mentioned manners, we focus on improving Reynolds equation. With distributed damping and stiffness, distributed transient EHL behaviors will be fast described.

Decomposition of Transient EHL force

This paragraph discusses in depth interacting mechanism on EHL surface and introduces new principle of simplification.

Generally, oil film stiffness and damping co-exist on EHL surface. Distributed behaviors are reduced to three concentrated interactions. It has been well understood that zero effect at exit zone on total load is due to cavitation caused by generation of negative pressure, and oil film also performs indifferently at Hertz's contact zone since central approach is dominated by Hertz deformation, so that, by neglecting these two weak interactions, oil film's effect on contact load is described at entry zone by a rough and direct integration of Reynolds equation, in the form of parallel normal supporting load [12],

1.5 4 6.66 min, min 00 1.5 ,

R u L u R L y y z iso is y h o Q h   (1. 1)
Inlet equivalent relation and central hertz relation turn one EHL pair to damping and stiffness coefficients, having been widely used to perform bearing dynamic calculation or to formulate bearing stiffness or damping matrices in the past.

However, this simplification fails to consider the deformation outside Hertz area.

Materials compress at hertz area (primary stiffness), but stretch at entry zone (extra damping produced by -inciting wings‖). Grubin [START_REF] Grubin | Fundamentals of the Hydrodynamic Theory of Lubrication of Heavily Loaded Cylindrical Surfaces[END_REF] predicted inlet oil film thickness according to this assumption.

Meanwhile, preceding analytical expression well matches highly-loaded contacts, but fails to serve weakly-loaded contacts, thought by Benedikt Wiegert, etc. Meanwhile, structural force equal to hydrodynamic force was organized by kinematic constraint, so that their new force model was expressed as a function of approach value, approach velocity and steady-state central oil film thickness. However, both film and static structural contact coefficients were derived separately according to Reynolds equation or Hertz theory, more daring an assumption of introducing a linear damping constant into the elastic structure was also specially used in his paper.

Although two preceding methods are well performed, surface loading is distributed around entire deformed area, interaction between elastic structure and lubricant should be considered, that is why Grubin's prediction was incomplete. Until full numerical method is introduced into the coupling calculation, accurate distributed solution can be obtained in a time-consuming way. Chippa, S.P. and Sarangi, M. worked on this method to extract lumped parameters, but their perturbation method was still based on steady-state results. Wherein, a same oil film perturbation was used along entire contact surface, still leading to rough estimations of integrated contact coefficients by directly accumulating each parameter.

In this chapter, elastic deformation keeps original calculating strategy; but lubricant film stiffness and damping forces decompose from Reynolds equation. The mutual effect between oil film and elastic deformation under transient EHL can be maintained by transient damping and stiffness distributions. Theoretically, this simplification is more accurate than two preceding methods.

Reynolds equation:

3 12 12 h dp h h u x dx x t             (1. 2)
Based on boundary assumption

               (1. 3) 32 12 ( ) 3 ( ) ( ) ( ) ( ) ( ) ( ( ) 12 36 12 6 4 3 ( ) ( ) ( ) 12 3
) ( )

) () ( h x h x h x h x ee h x h x h x dk u h h u u u film x x x x dx h h h dc x film dx h e e e hx e                          (1. 4)
Under steady states, formulation (1. 4) can be transformed into

                            (1. 6) 3 2 3 3 ( ) ( ) 0 0 0 0 0 0 0 ( ) ( ) ( ) 12

h dP dP dP h h dP h d h d h h

h h h hh h h P P u hh x dx dx dx dx x dx dx x t                                          (1. 7)
                   (1. 8)
To boil down to formula (1.8) and to respect differences between original and simplified equations, some comparisons among internal terms in original equation (1.6) tell the approximating conditions.

                        
It is found that only from inlet to exit of EHL contact can get ( ).

Further to this, dynamic distribution of oil film calculated by steady-state distributed stiffness and damping can clarify its non-robustness.

Steady-state equation and fluctuation equation of transient

                                                               (1. 9)
and come from steady-state EHL solution, and serve fluctuation equation.

P h h h h x dx x dx dx A B C D                                     2 0 2 dh dx   (1. 12)

Distributed perturbed h

 and variable h are substituted by discrete forms.

1 () 2 1 1 2 ( ) ( ) 22 () 0 1 1 1 1 ( ) ( ) ( ) 00 11 1 2 1 2 2 ( 2 ) ( 2 ) ( 2 0 0 0 1 2 kk d h h h tt k dx x k k k k d h h h h h t t t t k dx x k k k h h h tt t k k k k k k d h h h h h h dh t t t t t t t k dx dx xt k k k k k k k h h h h h h h dh t t t t t t t dx                                                  12 ) 11 2 () kk hh tt k xt                      (1. 13)
Where, the superscript k denotes the film thickness at location x k and the subscript t the film thickness at time t.

Substituting (1. 13) into (1. 12), h  h , , , A B C D Chapter 1 1 1 1 1 2 1 2 1 2 ( ) ( ) ( 2 ) ( 2 ) ( 2 ) 0 0 0 0 0 1 1 1 1 1 12 2 1 1 ( ) ( 2 0 0 0 ( k k k k k k k k k k k k k h h h h h h h h h h h h h t t t t t t t t t t t t t A B k t xt k k k k k h h h h h k k t t t t t C h D A B C h D t t k k k t x x t x k k k h h h t t t AB t                                                         1 1 2 12 ( ) ( 2 ) ) 2 2 1 0 0 1 1 1 1 1 ) ( ) ( ) 2 ( ) ( k k k k k kk h h h h h h h A B D k t t t t t t t A B C h t k k k k t t x t x t x t x E F G B D A B kk hh tt k k k t x t x x t H                                              ) 12 D C k x I kk E F G h H h k tt h t I             (1. 14)
Assume that two first increments of film thickness are zero at each state. Dimensional results under both steady and transient states:

Solution of Transient EHL contact

Transient EHL procedure starts from steady-state solution. Newton-Raphson method has been explored explored by Houpert to solve steady-state EHL, as seen in Appendix B. To take advantage of proposed contact stiffness distribution, a new steady-state EHL algorithm is performed instead.

Film thickness equation is,

2 ( ) * 2 4 ( ) ln 0 k R b x h x h b s x p s ds k E s x                   
Based on deformation-matrix method,

          1 11 1 22 22 22 ln (s)l 2 ln ( ) 1 2 n 1 xx n n p s x s ds p x s ds p D k k i i k i ik EE i i xx X X X X i k k k k D ik i                            Film thickness equation turns into, 2 ( ) * 2 1 n pD i ik R x h x h i     (1. 15)
Newton iterative method and the following frame constitute the solver. Assume that pressure increment along contact line relative to previous moment is p  distribution. Thus, the updated film thickness and its increment are,

2 * 0 2 1 32 3 0 0 0 12 0 2 ( ) 0 n x A h h p D i ik R i h dP h dP h Bu x dx dx k k k C P P P h h b k Hertz D P h h N             Chapter 1
2 () 2 ( ) * 1 ( ) ( ) ( ) 1 n p p D i i i x h x h t i h x h x h x t k R t                  (1. 16)
Real-time stiffness and damping distributions are ( ), ( )

k x c x (equations (1. 4)) establish, ( ) ( ) ( ) ( ) ( ) h x k x h x c x p x    (1. 17)
Combining (1. 16) with (1. 17),

1 1 1 [ 2 1 * ( )]( ( ) ( )) () 1 2 t t t x h h x k x c x p n p p D i i ik i t i R             (1. 18)
Besides (1. 18), conservation of load increment should be also satisfied.

2 ( ) ( ) b Hertz p Q t t Q t i N       (1. 19)

2017.09.18

Combining (1. 18) with (1. 19), the increments of film thickness and contact pressure under transient conditions can be obtained within a few steps by Newton iterative method.

Please note, in the dynamic process,  Real-time update of contact grids  Zero distributed stiffness under pure squeezing condition, but non-zero under both rolling and squeezing conditions Newton iterative method and the following frame constitute the solver.

2 ( ) * 1 ( ) ( ) ( ) 1 2 () ( ) ( ) 2 t n E p p D i i ik R F k G p P h P h ih x h x h i h x h x h x tt b Hertz Qt h H t N p Q i t                     Using equation (1. 18), equations E, H, G can turn to I, where [] D D ik  .         2 * () 2 * 1 , p i pp ii c R film k film t h c CK x h D h D t ii                           2 [ ( ) ] ( ) ( , )
,: ,: , * () Using proposed semi-analytical method, pure squeezing effect is calculated below, and its behaviours are studied after a sudden halting at rotating speed of U. Chapter 1

( ) ( ) ( ) Hertz p inv I N CIJ CK h c CK CIJ i P Hertz CK D diag i DC i i i i i K I            
Considering engineering facts, only squeezing action remains when the contact media stop rolling. The method proposed herein cannot work for a long simulation time since squeezing effect is fading away with time going by. Fortunately, some traces have been found at the early stage of this simulation that Hertz pressure distribution will replace it eventually.

As for oil film growth at contact center, lubricants are explained to fail to escape immediately due to the sudden stop of rolling motion. Practically, motion change needs time.

The pure squeezing action rarely exists if it evolves from normal lubricating operation.  With entrainment velocity increasing, pressure tends to distribute close to contact center. Oil film thickness increases and takes more effect.

Oil film transition between successive states is simulated with oil film stiffness and damping fluctuating at the same contact area. But pressure and thickness distributions involved should be updated in time to reproduce each new transient state.

Please remember, it is still a dummy transition. Different from the simulation on pure squeeze effect, in this case, the growth in oil inlet can be well predicted, but the recovery in oil exit cannot. Due to rolling, the oil film stiffness and damping distributed near the oil outlet Chapter 1 should be partly abandoned. If cavitation can be well located during rolling, the force jump at the oil outlet will naturally disappear in simulation.

Conclusion

An extended use of perturbation method is made to prove that steady-state EHL contact damping and stiffness distributions should not be used to analyze the transient conditions, let alone acquire equivalent total contact damping and stiffness coefficients for dynamic analysis.

                                           (A.4) Chapter 1

Appendix B Solution of steady-state EHL problem

Newton-Raphson method (Houpert etc.)

x p b h X P H H b p R HH w h U W G E H K ER HH W                       
Dimensionless form of Reynolds equation is given by,

3 ( ) 0 dP H ee f H K H i i i i i dX i           (B.1)
Thickness is given by 2 0 2 1

N X i H H D P i ij j j      (B. 2
)
i  is given by, 9 0.6 10 1 9 1 1.7 10

PP Hi i PP Hi       (B.3)
The viscosity-pressure is given by Roeland,

(unit,Pa s) exp{ [(1 ) 1]} 0 PP z iL i zP L         (B.4)
Where,

                   Iterative style is ( ) ( ( )) ( ) ( ) ( ) ( ) , 1, , 0 () 0 2 N f f f o n o n o n o i i i H P H f i N e e j i H P H e e i j                   (B.5) Load balance condition is ( ) ( ) 2
N f f f f N N N N H W H P P H e e N cc N                                                                                      o   (B.6)
The numerical method of Gaussian elimination is used to get the solution.

Chapter 2

Chapter 2

Calculation on distributed forces and moments of bearing [START_REF] Hartnett | A general numerical solution for elastic body contact problems‖[END_REF] [START_REF] Hartnett | A general numerical solution for elastic body contact problems‖[END_REF]and Jan M. de Mul, Joost J. Kalker (1998) [29] and others. They solved it in a similar way. However, pressure distribution in each slice was differently assumed to be semi-elliptical, uniform and parabolic. Technically, the first needs a time-consuming numerical calculation; the second presented by an analytical formulation causes a gross error; the third is also formulated analytically, but needs more discussion on potential error of contact width. Following their steps, the author of this thesis also develops an improved calculation by using some transverse rectangular strips to symmetrically replace elliptically-distributed pressure shape.

Compared to second approximation method, first approximation method is less popular in accuracy, but more in efficiency. In this chapter, first approximation method of R.Teutsch Palmgren and Kunert empirical method, others analytical method; Palmgren [START_REF] Palmgren | Grundlagen der Walzlagertechnik‖[END_REF], Kunert [START_REF] Kunert | Spannungsverteilung im Halbraum bei elliptischer Flachenpressungsverteilung über einer rechteckigen Druckflache[END_REF],

Houpert [START_REF] Houpert | An Engineering Approach to Hertzian Contact Elasticity-Part I[END_REF] and Teutsch R. [START_REF] Teutsch | An Alternative Slicing Technique to Consider Pressure Concentrations in Non-Hertzian Contacts[END_REF]obtained explicit solutions, while Lundberg [START_REF] Houpert | An Engineering Approach to Hertzian Contact Elasticity-Part I[END_REF], Tripp [START_REF] De Mul | The Contact between Arbitrarily Curved Bodies of Finite Dimensions[END_REF],

Dinnik [START_REF] Rothbart | Mechanical Design and Systems Handbook[END_REF], Kowalsky [START_REF] Rothbart | Mechanical Design and Systems Handbook[END_REF] and others acquired implicit solutions.

Incorporated in this alternative technique, roller discretization (slice number) is also crucial to solution accuracy as found by Singh and Paul [START_REF] Teutsch | An Alternative Slicing Technique to Consider Pressure Concentrations in Non-Hertzian Contacts[END_REF]. An average technique (ever called as Method of redundant field points) was then adopted [START_REF] Teutsch | An Alternative Slicing Technique to Consider Pressure Concentrations in Non-Hertzian Contacts[END_REF]. A similar strategy of normalization by mean values of all weighting functions was used by R.Teutsch [START_REF] Teutsch | An Alternative Slicing Technique to Consider Pressure Concentrations in Non-Hertzian Contacts[END_REF] as well.

The accuracy of AST is skillfully discussed in this chapter. The corresponding strategies and another reduced roller contact algorithm are proposed to maintain its technical robustness and convenience to engineering problems.

Chapter 2

With bearing questions differently argued, plenty of works are also differently performed to consider from normal contact to tangential contact and from contact model to bearing model. However, effective lubrication functionally ensures little tangential force. Bearing force is transmitted by normal contact rather than by tangential contact considering order of magnitude of force. Although insufficient and uneven oil film forces could result in roller slip and roller skewing respectively, tangential oil film damping, roller/rib resistance and cage constraint would ensure assumption on the normal contact in direct bearing model.

In order to reduce repeated calculations for each roller balance, as well as to easily solve bearing stiffness matrix, direct bearing model (Stribeck, Harris, Luc-Houpert 1997 [START_REF] Houpert | A Uniform Analytical Approach for Ball and Roller Bearing[END_REF]) is technically convenient to use, compared to equilibrium loop model (Liu, J. 1989 [START_REF] Liu | Equilibrium and Associated Load Distribution in Ball and Roller Bearings Loaded in Five Degrees of Freedom While Neglecting Friction-Part I: General Theory and Application to Ball Bearings[END_REF],

D.Né lias 2003 [START_REF] Bercea | A unified and simplified treatment of the nonlinear equilibrium problem of double-row rolling bearings -Part 1: rolling bearing model[END_REF], etc.) and dynamic model (Gupta [START_REF] Gupta | On the dynamics of a Tapered Roller bearing‖[END_REF], etc.). And, it is wise to choose relative motions between races as model inputs, in which case, by being part of dynamic system, races are balanced dynamically instead of statically.

In this chapter, the works to better predict bearing stiffness matrix include, 1. First and second approximation methods are versatilely used to develop alternative slicing technique (AST) and reduced roller contact algorithm (RRCA), with the error caused by AST explained.

Improvements of second approximation method;

 Semi-analytical method of flexible rectangular pressure strips inside slice is versatilely used to bridge an analytical and fast way to calculation;

 Semi-empirical method of enhanced load-deflection relation along slice thickness is used to stably obtain roller moment distributions; 2017.09.18

Deformation distribution in direct bearing model

With five known displacements applied on the bearing as showed in Figure2. 1, initial and general relative positions of cone and cup are shown, and the rollers occupied in loading area undergo elastic deformation which will be calculated in this part.

To develop general codes, initial geometric relations in a bearing are prescribed.

 When dy and dz are nil, central lines of cone and cup coincide with each other.

 dx equal to zero corresponds to both raceways just touching rolling elements.

Please note, the same -zero‖ convention still works when describing a double-row rolling element bearing.

Figure2. 1 Geometric relation in taper roller bearing (Houpert) Using Euler angles, relative rotations between raceways are described. Under the provision that rotation is successively around z axis and updated y axis, additional displacement of roller/cup contact position sin '/cos cos ( '/cos sin )cos ( '/cos sin )sin

x D x y R x o z Rx o                        
caused by tilting motions between races is given.

d d d d d d x x y z y z y x d y d d y y z z d d d d d z z d y z y z y z                                                        (2.1)
Deflection calculation follows the assumption made in the model, so that different manners in which deformation distributes in the roller should be distinguished between different models to clarify the possibly induced error.

Chapter 2

In both equilibrium loop model and dynamic model, deflection is measured between geometric centers of mating bodies, while, in direct bearing model, displacement of contact point at mating surface is directly used without considering the shape of contact bodies.

Herein, the effect of contact body shapes on deflection distribution in direct bearing model is considered. Roller/cup interference along roller length at azimuth angle  is suggested as

( , ') ( ) sin {( ) cos ( ) sin ... ...[( ( '/cos sin ( ) tan ))(1 cos ( '))]}cos x d d d d d d x x y y z R x d o x x z d x                           (2. 2)
It is specifically proposed that in case of tilting motions, the roller needs to deform curved part in cup. This part overcome by additional loading will recover with time going by, as shown by 2

A initial in Figure2.2. Proposed deformation consists of (a) displacement of contact point at azimuth angle  (penetrating into a wedge with conical degree  ), and (b) curved part before penetrating into a wedge.

As shown in Figure2.3, contact curve C 1 C 2 in the cup incorporating contact point ( , ') x  will axially move due to ,

dd xx  rather than , , , d d d d y y z z   .
Where, '

x is abscissa along the roller length.

d d d d y y z z x R x d d o x x              (2. 3) Figure2.2
2017.09.18

Calculating strategies for Bearing force and moment

In one roller, load applied on contact point x will be calculated herein by alternative slicing technique (AST, used by R.Teutsch) instead of classical slicing technique (CST) to investigate interaction among N slices in roller of length L.

x=-L/N+0.5L/N: L/N: L/N+0.5L/N

(2. 4)

Classical Roller Contact Calculations (Classical Slicing Technique)

In general use, deflection-load relation is described analytically and explicitly.

L i Q        (N, mm) (2. 6) (Palmgren, 206 , 0.3 1 2 1 2 E E Gpa      ) In slicing technique, p Q N 
is exercised, and sin( ) 0.99 sin( )

Q of i c Q o i f       
is used to revise the mechanical relation in a tapered roller (Figure2.5).

Figure2.5 Loading in a tapered roller

Slice/unit force is given by [START_REF] Harris | Rolling Bearing Analysis[END_REF], using a linearized fraction of total roller stiffness calculated from Palmgren's equation to describe each slice in case of misalignment.

    1.1 1.1 1.1 0.1 ( ) 1.1 p x k k x k x centre centre centre centre lc lc lc                 (2. 8) 0.1 2 ( ) (0) 1.1 centre M p x p x k lc x       (2. 9) 1.1 2 () 2 centre L p p x dx k L lc L     (2. 10) 2017.09.18 0.1 2 1.1 2 () 2 2 centre L M M x dx k L lc L     (2. 11)
Summation method can also provide the load and moment distributions.

(2. 12)

Equivalent forces and moments at all roller centers are transmitted to tan

R i
 point in cone central line. Within scope of theoretical mechanics, bearing loads are calculated. no

0 0 0 N Q F x x N Q F y y N Q F z z                                (2.
N Fp xi N Fp yi N Fp z i i i i                          ( 1) 2 1 ( 1 sin cos 1 cos cos 1) 2 1 L N N Mp yi i N i N N M i N i p z L i                            tan 0 tan 0 N M R Q M y i z y N M R Q M z i y z                              Cette thèse
                                                     (2.
                                (2. 15)
Where,

              (2. 16)
Hence, total load that cup applies on cone is derived according to Figure2. 6.

F F F x lx rx F F F y ly ry F F F z lz rz         ( ) ( ) ( ) (
             (2. 17)
It is exposed that load peaks (summation method) will be weakened compared to integrated results when bearing rotates, for higher order items in -generalized binomial expansion‖ vanish in integration method. Summation method performs more accurately in the case of gross tilting motions (see Figure2.8). CST enables each slice to share linearized fraction of single roller stiffness, but not to interact with other slices. Differently, approximation methods focus on analytical and direct interacting relations between distributed load and deformation. Second approximation gives an accurate 2D relation, while first approximation describes an inaccurate 1D relation.

To compensate for this inaccuracy in 1D relation, distributed load and mean load is bridged indirectly. Average slice load is described by both linearized fraction of single roller load and weighted average value of each slice load.

In the former description, linearized fraction of single roller load is, In the latter description, weighted average value of slice load vector is, Where, basic weighting function is generated newly,

    () 1 px N NN e ee Sq w      (2. 20) 2 2 0.9 0.39[4(1 )/ 4(1 )/ ] 0.9 0. 12 12 8 Q L EE         206 , 0.3 1 2 1 2 E E Gpa      Chapter 2
  1 ( ) 1 pp coef i i i r ij i pp j i j            (2. 21)
With influence coefficient reduced to,

1 () 1 () 1 , 4 jk coef j coef r jk w jk jk LN j                      (New) (2. 22) 1 , 4 jk r jk w jk jk LN coef coef                    (R. Teutsch) (2. 23)
When j=k, the influence is considered between slice edge and slice mid-plane for two half slices.

Different from the half-space theory, a force at slice j contributing to the deformation in another slice k was ever interpreted by first approximation suggested by Singh and Paul [START_REF] Teutsch | An Alternative Slicing Technique to Consider Pressure Concentrations in Non-Hertzian Contacts[END_REF].

Contribution of pressure at slice j to the deflection in another slice k decreases by 1/r, where r is the distance between slice j and k.

Using different deflection-load relation coefficients, influence function between slice load and average slice load can be predicted differently.

E.g. Using Palmgren relation, the governing equation is,

    1 = 1 , 1. ww jk jn k lc N we j k jk ww nk nn n q        
The alternative strategy of influence coefficients in bearing calculation includes another three aspects: the manner of roller balance, the selection of slice number and contact relation 2017.09.18 incorporated. This comprehensive consideration is responsible for a clear bearing prediction on different applications.

This chapter removes rollers from the equilibrium loop, and uses direct model to perform load calculation. Influenced deformation on both roller sides are adjusted by a most stable manner, in which each slice is balanced in a same way. E.g., in cylindrical roller bearing, a pair of equal loads is applied on both slice sides; in tapered roller bearing, a same ratio c of load pair takes place on each slice. The manner ensures entire roller balanced in the safest way.

Based on AST, the simplified line-contact system involved in tapered roller bearing is architectured in matrix form.

        11 11 1 . 1 1.1 1 1 k c c k lc lc q S q S i w i o NN N N N N ee wo e e e e NN ee                     1.1 1 () 1 1. 1 1 1 1.1 (1 ) io qS o k lc N NN e ee w c N e                  (2. 24)
Load intensity in one roller is observed below with different slice number considered. The influence of slice number on roller load is also investigated.

Figure2.10 gives slice number-roller loads relation curves. It illustrates that each load is nearly linear to the common (decadic) logarithm of slice number. So it could be objected that adequate slices are more accurate in first approximate use than a small number of slices.

To explain the phenomenon, a compared governing equation is supplemented.

    1 1 0 1, 1 = 1 0 .1 , 1 ww w jk jn kk k lc N e ww nk nn w n k nk q                  (Compared) (2. 25)
This compared governing equation and Figure2.10 show that the weighted average function in first approximation is the cause to this natural error.

An accurate prediction on integrated roller load is very difficult, but a reasonable depiction of distributed roller contact load tendency can be accessible and compensated in use more or less.

A simplest way in engineering is to apply load normalization method by using experimental value exp

Q . () ( ) exp () q predict i q Q correct i q predict i  
More complicated works can be also performed to confirm slice number. Load concentration is reflected in computed roller force and moment. The experimental manner can 2017.09.18

help to find a good agreement with a required tendency of load distribution, by comparing either measured moment or force with numerical result with different slice number used. This work should be numerically feasible, as Figure2.10 tells that each load is monotonically related to the common logarithm of slice number, and slice number can be searched at its law.

Compared to slice number, the error comes more from fundamental deflection-load relation. When the curvatures and thicknesses of contact bodies are considered, the formula in Appendix C can be used. E.g. in railways axle-box bearing, the Palmgren formula performed in this chapter can be adopted, since the corresponding structures are all very thick.

Alternative roller contact calculation is in essence based on first approximation. With

Palmgren formula used, it is a semi-empirical method. Besides AST, a semi-analytical method is also proposed in next two paragraphs.

Complex Roller Contact Calculations (Second Approximation Model)

Similarly, a second approximation used in slice technique can also fabricate a group of non-linear equations.

According to Boussinesq-cretu theory [START_REF] Boussinesq | Application des Potentiels à l'Etude de l'Equilibre et du Mouvement des Solides Elastiques[END_REF], the deformation at any point is given by integrating the pressure field

                 (2. 26) 1 1 1 12 () 1 2 12 K GG     1 1 2 1 2 12 () 2 2 2 2 12 K GG    
It can be solved numerically by deformation matrices and equivalent pressure in each finite and adequately small rectangular contact patch.

( , ) ( , ) , 12 (0,0) (0,0) * () D D f ij i N j N N N K P K P x y jj ii jj        (2. 27) Chapter 2 1 2 2 ( ) ( ) ( , ) ( , ) ( , ) ( , )
xa y jj j dxdy xa y

x x y y jj j i i

F x x a y y F x x a y y i j i j i j i j j j F x x a y y F x x a y y i j i j i j i

D j j ij j                                   * * * ( )cos ( )sin 2 2 ( ) ( ) ( , * * ) ( , ) ( , ) ( , )
xa y jj j x x y y i i dxdy x x y y xa y i i jj j F x x a y y F x x a y y i j i j i j i j j j F x x a y y F x x a y y i j i j i j i j j

j j D i                                       
Where,

1 2 2 2 2 ( , ) ln( ) ln( ) 22 F X Y dXdY X Y X Y Y X X Y XY          cos sin * 2 2 1 2 2 1 ( , ) [ ln( ) tan ( )]cos [ ln( ) tan ( )]sin 22 X Y Y X F X Y dXdY Y X Y X X X Y Y XY XY             
P j , as an equivalent normal load, occurs at the area of 22 a j , whose center is ( , )

xy j j .
Assuming that 1. Roller is divided into N slices (Figure2. 11) with N indent values discretized. These values are more determined by distributed pressure along roller length than that along roller width.

2. Longitudinal distribution of roller load is more related to tilting than skewing.

3. Both mating bodies of the same material,

        2 2 2 2 1 ln max, ' 2 2 a j y y y y x i j i j x dx Ea j a p ij j y y y y x i j j i j                               (2. 28)
Another versatile method simplified by the analytical equations uses a uniform pressure strip (Figure2.13) to replace Hertz distribution inside each slice. Herein, an alternative model is developed based on them.

Simplified model

i i j             (2. 29) 2 2 ( ) ( ) ( , ) ( , ) ( , ) ( , )
a y j j dxdy a y x y y j j i F a y y F a y y i j i j j j F a y y F a y y i j i j j

D ij j                             Figure2.
     is one resultant element in   1 j N e   .

Alternative model

Compared to empirical results (Figure2.17), enhanced results reveal that mean Hertz pressure performed in calculation led to overestimated results. Besides uniform, semielliptical and parabolic pressure distributions involved in the preceding calculations by other researchers, several rectangular pressure strips symmetrically superposed on each slice is also proposed in this chapter to replace semi-elliptically distributed Hertz pressure. An evident advantage is that distributed loads and deformation will be bridged analytically.

In this alternative method, rectangular pressure strips can be organized differently.

Below tentatively gives three situations (Figure2.14~2.16). A better use of proposed method can be claimed. The upshots show that the latter two methods of   A versatile idea can be employed in the following way. It is relevant to think that it is significant to have a more detailed description for the contact points of large deformation, but As such, the condensed matrix will be modified with a product factor c j . The product factor is used to describe the importance (strong or weak) of load distribution to each slice, in order to reduce the calculating scale of elastic deformation in equilibrium loop. If the load in j th slice contributes less to others' deformations, fewer strips can be used for slice j, it means that each column of {D ij } can be specially selected.

    1 1 2 2 max 1 1 1 n k c k c D D D sj sj sj k c p sj N e sj NN sjn n N e j N N
N N e e e e e e

                           (2. 35)
Where, sj donates j strips used in the deformation matrix, and j varies from 1 to n.

, kk s sj are constant coefficients corresponding to different strips.

0 1 0 T c j slice j         
is the product factor of selecting j th column of influence coefficients matrix.

Equilibrium loop method

Besides interaction among slices, unequal loads generally acting on two sides of each slice should be also comprehensively considered in the equilibrium loop model.

Chapter 2

Using discretized pressure   mean, 1 p j N e  distributed along contact line as variable in Hertz line contact relation, deflections distributed at both inner and outer races are calculated.

        2 ' 1 1 mean, m 2 en ' , 1 1 a E N N NN layers e e ee E N N NN layers e e e i i i k D p j j o o o k D p j j e                               (2. 36)
In tapered roller bearing, different contact geometries at inner and outer races participate in calculation. ' max, 2

pE s p jR x   p s R x  max, max, i p o j Rc x o i p R x j   , ( i p s c o p s  ) (2. 37) 2 ' max, R x ap j E j  max, Rp x j  i a i j Rc x o o a R x j   (2.
N i D N ee     is calculated explicitly.         ' ) 2 1 1 ' 2 1 1 ( ( ) E N N NN layers e e ee E N i i i p inv k D j j o o o p inv k D N NN layers e e ee j j                               (2. 39)
In equilibrium loop, geometrical relation is Moment equilibrium around roller end/flange contact point is written by

     
( sin ) ( ( ) sin ) 22 F LL h D h w M F M i i o o        (2. 42)
Where,

22 22 1 1 1 1 ii pa jj j i oo pa jj j N N e Fp ii N N e Fp oo i                          1 1 ( 1) 2 1 ( 1) 2 1 L i N i L i N i N N Mp ii N N Mp oo                     (2. 43)
Roller end/flange contact position is located by distance h. In this method, absence of force distribution along slice thickness is the leading error in moment calculation. The more slices are used, the more accurate moments are. When the number of slice is sufficient to fully describe the force distribution, the moments (M y, M z ) will stay stable. Theoretically speaking, the series {Lj} affect the synchronous stability of moment M and force Q calculations. {Δp• Lj} decreases to{0} in the unique case of adequate slices in calculating stably distributed {p}. But too many slices in bearing calculation are not realistic.

Semi-empirical method of stabilizing roller moment

2017.09.18

A semi-empirical method is suggested to improve the stability of moment calculation.

Herein, a constructing function is introduced to simulate pressure distribution along slice thickness. The deformation along the contact length is linear to contact coordinate. As such, the constructing function is linearly related to coordinate.

0 1.1 (( ) ) 0 p x a new x x ii xb i        (2. 45)
Please note, unit of

() p x newi is Nm, not 2 m N .
Interestingly, any constructing function with a power of 1or 2, etc. cannot ensure the computed moments to stay stable with more slices involved in calculation, expect when a power of 1.1 is chosen. Below gives the constructing idea in depth.

1.1 0.1 1.1 0.1 0.1 1.1 1.1 1.1 11 1 1 1 2 p k k n n n n n n n n n lc lc                                 
Unequal (non-average) interaction among slices is considered to significantly develop new distributed stiffness along each slice length.

                          (2. 46)
Stiffness for each slice is no more constant, dependent of enhanced pressure values as well as slice deflection distribution. The constructing function in fine proves the correctness of the power 1.1 after integration of (2. 46).

Besides power value 1.1 in constructing function, the integral coefficients , a i b i should also be calculated by both load continuity and load conservation inside each slice.

( ) ( ) 0 _ 1 1_ 22
( ) ( ) 0 _ 1 1_ 22 LL p x p x x newi i slice center newi i slice center NN LL p x p x x newi i slice center newi i slice center NN                 (2. 47) Chapter 2 _ 2 ( ) 2 2 _ , 2 L x i slice center N p x dx a newi i L x i slice c p mea en i t N n er         (2. 48)
Using a simple matrix operation, the coefficients of each slice are calculated recursively and successively, and resultant loads are established. Figure2.23 presents updated moments, to find that a small number of slices are sufficient to obtain stable results. The introduction of slice stiffness recalculation is in fine proved to be feasible in bearing engineering. Please note, in roller bearing stiffness and damping matrices, such a relation can replace uniform stiffness formula to perform a better prediction in system vibrating characteristics.

L x i slice center N N M y p x xdx newi L x i slice center N L x i slice center N N M z p x xdx newi L x i slice ce N i i nter                                  (2.

Reduced Roller Contact Calculations (1 st and 2 nd Mixed Approximation Models)

The first approximation (1/r) suggested by Singh and Paul should be used carefully, since its interacting distance is wide. A more precise work is efficiently performed below.

In Boussinesq work [START_REF] Boussinesq | Application des Potentiels à l'Etude de l'Equilibre et du Mouvement des Solides Elastiques[END_REF], deformation at any point is given by integrating the pressure field ( , ) p x y over the effective contact area ( Sc ), 2017.09.18 Using Hertz line-contact theory, Based on this reduced algorithm, the simplified line-contact system involved in tapered roller bearing is architectured in matrix form.

                   (2.
        ( (c ) , , , , ) 1 1 1 1 i o i o K w q w w w i o o j k j k j k q K q N N N N
                                   (2. 52)
s-approaching value decides interacting distance incorporated in w ii calculation, and determines availability of new influencing coefficient.

LN LN j k n LN LN w j a L N j k n N r aL r                                 (2. 53) , 1 , , , , ( , ) ( , ) 1 ( 0
                                        (2. 54)
Where, It is investigated that a) a very low coupling value n is sufficient to obtain accurate load intensity along roller length; b) in the case of large slice number, use of more second approximation influencing coefficients w j,k (a slightly higher n) will make reduced algorithm more stable.

1 2 2 2 2 ( , ) ln( ) ln( ) 22
F X Y dXdY X Y X Y Y X X Y XY          -0.
F a F a F a F a j a a L N L N w jk q q a a L N L N L N L N L N L N L N L N j j j b F a F a F j j                                       0.571 , ) ( 0. 0.5 / 0. 571 , ) 5/ a F a j j L N L N              (2.
L N a L N a L N a L N w jk q r a L N a L N a L N a L N                                          ( 2

Analysis on Bearing Stiffness Matrix

The influence of roller edge effect on bearing transmission lies in angular motion items in stiffness matrix. Under tilting motion, asymmetric pressure concentration will cause considerable moment acting on roller center.

2017.09.18

Within two rotations, a tapered roller bearing is calculated, using slice number 31 and setting five inputs: Figure2. [START_REF] Mevel | Accuracy in Tapered Roller Bearing Equilibrium‖[END_REF] Comparisons between CST and AST stiffness in rotation Figure2.30 reveals that phase (distribution of rollers in bearing chamber) will cause more evident periodical fluctuations of stiffness using AST. CST results build up a symmetric stiffness matrix, while AST results produce an asymmetric matrix, but the amount of asymmetry does not deviate very far. Besides, coupled angular stiffness has a more apparent change using AST.

K K K K K xx xy xz x x yz K K K K K yx yx yx y y yz K K K K K zx zy zz z z yz K K K K K x y z y y y y y y z K K K K K x y z z z z z y z z                             ( , , , , ) ( , , , , ) () ( , , , , ) ( , , , , ) () ( , , , , ) ( , , , , ) ()
                                       0 2 4 3.5
To terminate this chapter, nonlinearity of bearing stiffness coefficients is considered into resonant frequency.

Described under zero phases, nonlinearity of some stiffness terms to inputs: In some special applications (like axle-box bearing in railways), bearing inputs are of high frequencies but low amplitudes. So that, system resonant frequency f can be predicted to understand how these bearings support the excitations, using linearized static stiffness.

Where, m and J are the suspended mass and moment of inertia.

is a damping factor, on the amount of 0.05 for example. Figure2.36 still pictures that rotation causes a fluctuation, but plays less importantly than displacements. The results also predict that resonant frequency changes much more steeply at lighter preload than at heavier preload, so that adequate preload is beneficial to improving stability.

Lots of frequencies exist simultaneously in motion. Some components are unexplainable at times, now the ones resulted from five inputs are more accurately predicted.

Conclusion

Bearing transmission is better predicted by considering roller edge effect. R.Teutsch's and new developed ASTs are introduced based on first approximation method with their errors explained. Meanwhile, reduced roller contact calculation is also developed by mixing first and second approximation methods. At last, bearing stiffness components and relevant vibrating characteristics are discussed. This paper in fine provides a numerically stable and engineering-adequate bearing model of suitable complexity. Such convenient calculations will be possible to better predict the performances of bearing and its adjacent components by implementing in any non-linear F.E.A package or any complex M.B.S package.

Nomenclature

Appendix C Deflection-load relationships for inner and outer race-to-roller contact

Palmgren [START_REF] Palmgren | Grundlagen der Walzlagertechnik‖[END_REF] 0.9 5 3.85 10 0.8

Q io L      
Kunert [START_REF] Kunert | Spannungsverteilung im Halbraum bei elliptischer Flachenpressungsverteilung über einer rechteckigen Druckflache[END_REF] 0.925 5 4.05 10 0.85

Q io L       Houpert[36] 2 1) 0.074 1/1.074 ] 0.2723 2 1)
1/1.078 0.27835

( [ 0.078 ( [] 1/ d m E Q i Q L L t o Dd m E                  
Teutsch [

RR E b E b Rt E b E b Q i L Q o L                           Dinink[39] 22 1 1 2 1 1 2 1 1 2 2 ( ln ) ( ln )] 2 ( ) [ 2 () 33 12 22 1 1 2 1 1 1 1 2 ( ln ) ( ln [ )] 33 12 RR E b E b Rt E Q i L Q o b L Eb                            Kowalsky[40] 22 1 2 1 2 1 1 2 2 (0.407 ln ) (0.407 ln )] 2 ( ) [ 2 () 12 22 1 2 1 1 1 2 (0.407 ln ) (0.407 l [ n )] 12 RR E b E b Rt E b E Q b i L Q o L                         
Chapter 3

Chapter 3

Calculation on nonlinear bearing -rotor system

Introduction

Bearing technology is a kind of old knowledge, the traditional questions-dynamic interfacial parameters are always welcome in both bearing and relevant transmission applications since accurate and efficient coupling analysis of dynamic relations and subsequent design of structural composition and parameters are required according to diverse demands. Bearing interference characterized by coefficient matrix and its coupling approach with rotor system are of prime interest in this chapter.

From bearing stiffness matrix …

Numerous study efforts have come to bearing stiffness matrix. It has ever been performed by Jones From modeling methods of integrated rotor-bearing system … Dynamic performances of rotor system depend on different configurations, bearing model can be assumed to be rigid, linear [START_REF] Lund | Stability and Damped Critical Speeds of a Flexible Rotor in Fluid Film Bearings[END_REF] or nonlinear, and rotor can be also modeled as rigid system, lumped mass system (with massless beam) or continuous mass system (Euler-Bernoulli beam element, Rayleigh model, Timoshenko beam element [69], etc.). At the same time, some effects of stiffness, damping, rotary inertia, gyroscopic moments, or shear deformation of different unit on rotor system have also been studied to different extent. In this chapter, nonlinearity of bearing is emphasized, while rigid rotor is employed.

To mathematically connect sub-models composing rotor-bearing system, transfer matrix method (TMM) and finite element method (FEM) are two mainstream methods. The marching idea in TMM makes the incorporated matrix present at a fixed dimension, beneficial for programming development. FEM should be the earliest method, succeeds in accuracy, but fails in efficiency since the selection of unit needs a careful discussion.

As far as TMM is concerned, since it was first employed in linear bearing-rotor system by Prohl, it has experienced a long-term development, but it still needs a discussion when nonlinear bearing is considered.

Chapter 3

… to solving methods of integrated rotor-bearing system

Normally, both time and frequency methods are developed to solve the equations.

Compared to frequency methods, not too many arguments are about how to solve nonlinear bearing-rotor equations in a more efficient time-domain way, but higher focus is on parametric influence on system responses as well as discovery of rich nonlinear phenomena in it.

In time integration method [70], periodic, aperiodic, and quasi-periodic solutions can be obtained together with transient responses. However, steady state responses can only be acquired after transient responses. But, it still belongs to a qualitative method by watching the clock of periodic solution's appearance. To perform stability analysis quantitatively, shooting method is often used to get periodic solution, but it requires strict initial values and integration step.

In frequency domain methods, the transient responses are deliberately skipped so that steady state responses of desired frequencies can be directly calculated in a fast way. However, the complexity of non-linear items included in non-linear differential equations determines accessible treatments. Some treatments, like perturbation method, multi-scale method, classical HB method, etc. require application of Taylor expansion or curve-fitting methods, where truncation of high order series does introduce errors and limit periodic solution only obtained around Taylor expansion point. More efficient methods [71], [START_REF] Detroux | The Harmonic Balance Method for Advanced Analysis and Design of Nonlinear Mechanical Systems[END_REF] for non-linear questions-IHB, HB-AFT, and GHB were later developed. Different from IHB and GHB methods, HB-AFT method performs better technically in nonlinear bearing-rotor system by simultaneously interpreting unknown response and bearing force as combined harmonics.

This chapter undertakes to develop more practical bearing matrices and more reliable model-solving methods in a uniform and economic way.

Practically, we can find a wide usage of analytical model (Lim and Singh, Houpert, etc.) in gear box or other machines, which has been regarded as a fast and enough accurate use. But theoretically, in roller bearing, the integrated load for one roller ignored the high order load term after performing Taylor expansion of Plamgren's roller load relation [59], it would produce errors when slightly gross tilting angle of roller occurs, so that slice technique is intentionally used instead. And in ball bearing, varying contact angle affecting both 2017.09.18 deformation value and load transmission is also considered. The two points distinguish and also supplement main diversities in uniform matrices of ball and roller bearings.

Meanwhile, application of roller bearing matrices in time/frequency domain analysis on rotor dynamics system is performed respectively. In time domain analysis, solving method is improved to ensure a convergent result. While in frequency domain analysis, a practically accepted way to periodic solution with less programming effort is manually programmed, using the current popular HB-AFT (Harmonic Balance-Alternating Frequency Time) algorithm to form a uniform 3D code.

Nonlinear Models of Roller Contact Damping and Stiffness

The dynamic parameters for each lubricated contact in both roller and ball bearings are primarily extracted here. Different from ball bearing, description of both EHL and elastic contact requires an extra usage of slice technique in the roller bearing for analyzing load configuration.

Characteristics of Lubricated Contact

Figure3. 1 Equlivalent stiffness and damping of a lubricated contact

Equivalent contact stiffness at each contact can be simplified by

dF dF K K H EHL k d d dh K K EHL H       (3. 1)
According to Harris' relation [START_REF] Harris | Essential Concepts of Bearing Technology[END_REF]

i t t R e                              .
So, individual contact force including damping and elastic loads is

( ) F k i c       (3. 2)

Linear Model of damping and stiffness matrices

Finite differential method [60] has been used to approach stiffness and damping matrices by comparing bearing load vectors under different similar displacements with only tiny fluctuation for designated differential value. We can see the entire dynamic effect, but not internal transient parametric details.

Herein, stiffness and damping matrices are assembled by transient dynamic parameters from all available contact sources.

In such case, overall stiffness and damping matrices of a roller bearing are derived below,

{ } { } { } [ ] [ ] [ [ ( )]] { } { } { } 1 1 T Z N F Q ij ij K T T T d d i j ij               (3. 3) { } { } {} [ ] [ ] [ [ ( )]] { } { } { } 1 1 T Z N F Q ij ij C T T T d d i j ij               (3. 4)
And for a ball bearing

1 N  , they are { } { } { } [ ] [ ] [ [ ]] { } { } { } 1 T Z Q F i i K T T T d d i i             { } { } { } [ ] [ ] [ [ ]] { } { } { } 1 T Z Q F i i C T T T d d i i             Chapter 3
To fulfill the above-mentioned matrices, each contact load as well as equivalent contact parameters for one rolling element should be understood at first. In this case, resultant equivalent frequency-dependent contact parameters for a pair of roller slice's lubricated contacts can be evaluated by means of complex equation (Dietl 1997),

1 1 ( ) s k i c k i c k i c i i o o          (3. 5)
Where,  is the angular frequency, generally of the inner race, outer race, cage or ball pass frequencies, or combinations of two or more of these [77].

Load Distributions in Bearing

With five known displacements applied on the bearing as showed in Figure3.3 and For easy analysis, some initial geometric contact relations are prescribed in advance for the analytical model of tapered rolling bearing.

 When dy and dz are nil, central lines of both races coincide with each other.

 dx equal to zero corresponds to both raceways just touching rolling elements.

Please note, the same -zero‖ convention still applies when describing a ball bearing.

Geometric relation in tapered roller bearing

2017.09.18

                                                       (3. 6)
Arbitrary roller contact position is sin '/cos cos ( '/cos sin )cos ( '/cos sin )sin

x D x y R x o z R x o                                     (3. 7) With 1 ' ( 1) 2 N l x j N     
Total approaching value along contact line at azimuth angle is proposed as

( , ') ( )sin {( )cos ( )sin }cos x d d d d d d x x y y z z             (3. 8)

Geometric relation in ball bearing

The ball diameter is D, the race curvature radii are r i and r o with curvature centers i and o at a distance of R i ' , R o ' from rotation axis respectively. Inner and outer race centers are I and O.

the initial contact angle is α 0 .

Chapter 3

Figure3.4 geometrical parameters in ball bearing (Houpert) Inner and outer race curvature centers i, o are respectively located at

    sin 'cos 'sin 0 T T x y z r R R I I I i i i         ( )sin 'cos 'sin 0 T T x y z D r R R O O O o o o     
Different from roller bearing, change of ball-race contact angle brings about different deformation-load relation, while flexible bearing inputs lead to variation of the loaded contact angle α, affecting the paths of load transmission and local elastic deflection.

To pursue a uniform analytical approach, relative displacement between race curvatures centers i and o is calculated.

Similarly, tilting angles of outer race relative to inner race cause additional relative displacements between race curvature centers and update contact angles for different balls.

                                                           (3. 9)
According to a derivation in Appendix D, under varying contact angle, resultant interference of ball at azimuth angle owns a uniform expression with tapered type,

( ) {( ) sin [( ) cos ( ) sin ] cos } 0 0 d d d d d d x x y y z z                 (3. 10)
And, variation of contact angle is

Q o f i s Q o i f        
For tapered roller contact

k k ei eo K U G U G EHLC o o i i      

For ball type

Once load is obtained, each contact parameter in Part 3.2 will be acquired.

Formulations of Bearing Stiffness and Damping Matrices Load Transmission through tapered roller Bearing

Chapter 3

The force along the roller with N slices will be obtained. The total transmission load at roller contact center that the cup applies on at azimuth angle  is given by sin cos cos si c

1 1 s 1 no N Fp xi N Fp yi N Fp z i i i i                          ( 1) 2 1 ( 1 sin cos 1 cos cos 1) 2 1 L N N Mp yi i N i N N M i N i p z L i                            (3. 13)
Where,

0 n i k p ii  
acting on outer race,

0 p ii n k i   
 acting on inner race.

Along roller length, the force acting on a roller will be replaced by an equivalent force and moment at contact center (Houpert, T.C. Lim (ignores additional moment)).

Load transfer to roller contact center helps to concentrate all the forces at tan R i  point of inner race's central line, and additional roller moments in rollers can also be transplanted to tan R i  point of inner race's central line directly within the scope of theoretical mechanics.

Hence, total load that the outer race applies on the inner race assembly is given by 0 0 0

N Q F x x N Q F y y N Q F z z                                tan 0 tan 0 N M R Q M y i z y N M R Q M z i y z                              (3. 14)
The deformation values are far smaller than R i , and can be neglected when getting total moments.

Organized into

T ij d x ij d y Z N F ij ij K N l j d ij z i j N R i ij R i d y ij d z                                                                                                        (3. 15)
According to geometric and motion relations,

2 2 2 2 2 2 ( ) ( ) 2 2 2 ( ) ( ) F s k k w c k k w c ij i o o o i i k ij s k k s c c w o i o i              sin ij d x      cos cos ij d y       sin cos ij d z      sin cos
                                                                           cos( ) sin
                                                                           
The same derivation is also employed for damping matrix,

Chapter 3 sin 0 cos cos 0 1 [ ] [ { sin cos 0 [ 1 (
1) sin sin sin cos 2 cos sin cos cos

T ij d x ij d y F ij ij C N l j d z N R i ij R i d y ij d z                                                                                                      ]}] 1 1 Z N ij i j     (3. 16) 2 2 2 2 2 2 ( ) ( ) 2 2 2 ( ) ( ) F s c k w c c k w c ij i o o o i i c ij s k k s c c w o i o i              sin ij d x      cos cos ij d y      sin cos ij d z      Due to sin cos cos ( ) ( ) sin cos T d x ij d y dd y z y z d z                       d dd x xx d d d d d y y y y z dd yz dd d zz z                                        d d x x

Load Transmission through ball Bearing

According to Appendix E, contact angle affects load transmission and deformation fluctuation.

N Q x i N Q y i N Q z F F i i F F i i F i F i i                                             (3.
FF ii FF i NN M R R r y i i i NN M R R r i z i i i                                       (3. 18)
Organized into Overall stiffness and damping matrices cos cos sin sin sin sin ( cos cos (

sin cos 0 0 cos sin 0 0 [ ] [ { cos sin 0 0 ( ) sin cos ) 0 0 sin cos ) 0 0 R R r i i i R R r i i i d x i d y F i i i K d F z i i i d y i i d z                                                                                            }] 1 T Z i                     (3.
) 0 0 R R r i i i R R r i i i d x i d y F i i i C F d i z i i d y d i i z                                                                                }] 1 T Z i                                   (3. 20)
All related partial differential terms have the same formats as those in tapered roller bearing, except for

( ) ( ) ( ) FF ii F i i i i                 ( ) ( ) ( ) ( ) ( ) ( ) d dx dy dz d y z dx dy dz d d i i i i y i z i                                                 
Where,

                                           also for ( ) ( ) ( ) FF ii F i i i i               

Application in Nonlinear Bearing-Rigid rotor dynamics

Apart from the model configurations and modeling methods referred in the introduction, another very important part concerns how to solve models as well as to discover the hidden information between solutions and models.

Nonlinearity of bearing unit included in mathematical model used to be in the form of nonlinear force (3.21), mainly characterized by Hertz deformation-force relation, at times, also considering bearing clearance and element waviness.

( ) ( ) ( ) ( ) ( ( ), ( )) MX t CX t KX t f t f X t X t b rotor      (3. 21)
To boil down to format of linear equation, it can be also found in some publications that transient bearing stiffness and damping matrices participate in bearing-rotor dynamic calculation, instead of varying mean bearing matrices of coefficients (3. 22), using finite differential method.

( ) ( ( ( ), ( ))) ( ) ( ( ( ), ( ))) ( ) ( ) MX t C C X t X t X t K K X t X t X t f t bearing bearing      (3. 22) ( ( ), ( )) ( ( ), ( )) ( ( ), ( )) ( ( ), ( )) f X t X t f X t X t bearing b repl rotor K X t X t K X t X t bearing bearing XX ace        
( ) ( ) ( ) ( ( ( )) ( ) ( ( )) ( )) ( ) 1 1 2 2 M X t C X t K X t C X t X t K X t X t f t               (3. 23) ( ) ( ) ( 1) X t X t X t     (3. 24)
Where, C 1 , K 1 are constant coefficients, C 2 , K 2 are nonlinear bearing coefficients.

Better solutions with less programming effort are explored to understand how bearing interacts with rotor, using a most simplified rotor system described in Figure3.5, where a rigid rotor is symmetrically supported by a bearing. The matrices included are given here.

2 [ cos( )] () 2 [ cos( )] Q x m e t mg Q y ft m e t Q z M y M z                      m m Mm J d J d          , 0 0 0 1 0 0 C J p J p              , 0 1 K  , () T X t x y z yz     Jp =0.24505kg• m 2 L=1.1133m R=65mm Jd=11.98kg• m 2 m=116kgρ=7850kg/m 3
Chapter 3

In this 5-DOF system, the matrices M, C 1 are on behalf of mass and gyroscopic matrices;

Study on rotor responses considering nonlinear bearing i n time domain

Using Four-order Adams Method given in Appendix F, in the manner of bearing matrices, the governing equation (3.23) can be easily solved to recognize incremental variables () Xt 

. Four-order Adams method performs perfectly in calculation efficiency and convergence of () Xt  . However, the error from () Xt  will be continuously accumulated for () Xt, its fragile accuracy and convergence ask for a comprehensive solver.

A new solver of Incremental nonlinear dynamic equation

The error from () Xt  must be seriously considered after a lot of simulation failures, using current numerical methods (like Four-order Adams Method). Incremental governing equation ), should also consider each parameter-varying coefficient's fluctuation and its robustness at each integration step.

This chapter specifically proposes a new constructing function.

1 () 0 2 1 i i K K x x      is suggested to be replaced by 1 1 [ ( ) ( )] 00 2 2 11 ii i K K x x K x x         . It is qualitatively explained that generation of x i  considers effects of preceding states ,, 3 2 1 x x x i i i  
 and x i  itself on dynamic equation. This effect mainly reflects on varying bearing coefficients, and adoption of mean value between

1 () 0 2 1 i i K K x x      and 1 () 0 2 1 i i K K x x     
in the solver is one essential manner considering this effect.

Chapter 3 movement of red dash line. It is equal to the following equation.

( ) ( ) ( ) ( ) ( ( ( )) ( ) ( ( )) ( )) 1 1 2 2 ... ( ) ( ) ( ) ( ) ( ( ( )) ( ) ( ( )) ( )) 0 1 1 2 2 M X t C X t K X t f t C X t X t K X t X t actual actual actual M X t C X t K X t f t C X t X t K X t X t actual actual actual                           (3. 26)
No convergent drawback is in fine overcome by bridging to an equivalent question of the same convergence as static bearing mechanic balance question, where 0.5 is proved to be the right coefficient. And, obviously, other setting of unequal coefficients fails to get rid of the influence of the error inaccurate () Xt introduces to

( ) ( ) ( ) MX t CX t KX t
(demonstrated by red dash line). Meanwhile, it also well corresponds to simulation findings.

( ( )) ( ( )) ( ( )) ( ( )) 2 2 2 2 ( ) ( ) ( ) ( ) ( ) ( ) 0 11 22 C X t C X t K X t K X t M X t C X t K X t f t X t X t                 (3. 27)
These nonlinear variables are in fine predicted implicitly instead of explicitly, since it becomes a combination of dynamic marching procedure and static balance procedure. Now, by assuming boundary conditions, solutions are obtained.

=-2 - 0 ( ) ( ) ( )=0 t t t X t X t X t      , , , 0 Q Q Q M M x y z y z

  

Fundamental results on nonlinear bearing-rotor system

The different manners in which undamped or damped bearing and balanced or unbalanced rotor perform in healthy rotor vibration responses are studied. Although much focus has been put on numerical solver, more explanations about possible simulation results are also given to make proposed solver persuasive.

Undamped conditions

) ( ) MX t KX t G F preload    (3.28) ( ) sin( ) GF K preload X t A t KM M      (3.29) fgdgsd
Similarly, in nonlinear conservative system, it is found as well that bearing vibrations in Figure3.8 and Figure3.9 are also influenced by initial conditions ( ), ( ), ( )

X t X t X t ( = -2 - 0 t t t  , , ).
In Figure3. How influences of untransferable and transferable potential energies are distributed in responses is not clear using a random initial setting, however, in application, (0) X can be given by a non-zero value to represent the initial kinetic energy, which, in essence, is an equivalent potential energy. To further confirm them, preload and initial kinetic energy can be carefully calculated using first two groups of state equations. Here, no more discussion is performed.

Below give two situations uniquely changing initial potential energies. Combined with the findings on rotor dynamic characteristics from Figure3.9, it is again evidenced that initial transferable potential energy affects vibration amplitude and natural resonant frequencies, while preload affects vibration balance point and natural resonant frequencies. Furthermore, the impact of initial transferable potential energy will even become reduced due to damping, so, in practical engineering, it seldom takes many difficulties of getting rid of initial states' influence in identifying system as well as external excitations.

In fact, if the damping is considered, the initial potential energy and even some small excitations (like VC vibration) will disappear.

Damped conditions

Figure3.13 Responses of unbalanced rotor using constant damping (C yy =C yz =C zz =C zy =200 N/m•s)

Figure3.14 Responses of unbalanced rotor using proposed damping (Simplified steady-state damping coefficients in an integrated manner)

It is noticed that damping dispersive behavior does suppress some vibration. In In conclusion, two points rarely emphasized are specifically studied in this chapter.

1. Use of incremental differential equation delivers a uniform but also persuasive engineering interpretation in case of nonlinear stiffness/damping coefficients.

2. The Runge-Kutta method used to perform a better convergence in nonlinear stiff equations, since the integration step in this one-step method can be changed as accuracy requires, however it also naturally produces an unpredictable convergence time. The four-order Adams method used in this chapter is a multi-step method with unchanged integration step, but easily diverges in nonlinear questions. The success herein is mainly profitable with the specific treatment for varying simulation error shared by inertial loads and gyroscopic moments (for which part, no safe method in the past, like Jacobian matrix in static equation, can help to converge it, unless stepvarying control of dynamic accuracy). Although several iterative times are taken in static calculation included at each integration step, less integration steps out of a max motion period are wasted in the whole convergence. Only when damping exists, some more periods are needed. The described enhancement can be observed from all preceding results.

2017.09.18

The proposed time method has a short convergent period, only a small number of periods are needed to observe system behavior. However, initial states still need a discussion on distribution of initial transferable potential and kinetic energies, as well as constant preloads.

The authors in this thesis also believe that a combination of the Runge-Kutta method and proposed constructing function can be a good choice to settle it.

Study on Rotor responses considering nonlinear bearing in frequency domain

The dynamic relationship between local parameters with global rotor working behaviors is the core of designs and the destination of simulation predictions, and inversely important to identify and control each element in bearing-rotor system.

More than what we have observed from the preceding pictures, global working behaviors reflected in transient and steady state responses are much richer in reality. However, it is sure that periodic solutions and stability are technically important for a nonlinear system.

Again, local parameters consist of those describing each internal element (disk, supporting or suspending bearings, etc.), and those depicting external excitations (rotor gravity, rotor unbalance, environmental vibrations, etc.), so that parametric influence on systematic stability is a diverse question. Accordingly, acquisition of periodic and nonasymptotic solution needs a more direct path in frequency domain.

HB-AFT applied in nonlinear bearing-rotor system using uniform matrix format

Using truncated Fourier series to describe unknown displacements, bearing forces and other external loads acting on rotor,

[ cos( ) sin( )] 0 1 K X A A jwt B jwt jj j      (3. 30) [ cos( ) sin( )] 0 1 K F C C jwt D jwt jj j      (3. 31) [ cos( ) sin( )] 0 1 K F E E jwt H jwt jj external j      (3.
T T T T T T T A A B A B A B P KK T T T T T T T Q C C D C D C D KK         Substituting [ sin( ) cos( )] 1 2 2 ( ) [ cos( ) sin( )] 2 1 K d X jw A jwt B jwt jj dt j K dX jw A jwt B jwt jj dt k                   
into dynamic equations, the same harmonic compositions in Q and P are balanced with a Galerkin projection. The resultant format in frequency domain is yielded in general way.

0 j   0 0 0 0 0 0 0 0 T A C mg E      1 j  2 ( ) M C 2 C ( ) M A C E jw jw j j j G B D H j j j jw jw G                              
() 1 dynamic FFT Q F X P Q Q i bearing FFT i i     
Figure3.15 Iterative paths of fixed point P included in frequency method

The straight lines in the pictures describe the relation between P and Q hidden in systematic dynamic relation as linear relation. (Remember, when k=0, the elements in Q have nothing to do with those in P. This linearity will not be influenced). This finding reflects that a function as rebuilt in time domain is not necessary. In strategy 1, direct iterative procedure apparently fails to get a converged P vector according to the trends of blue arrows in Figure3. [START_REF] Wensing | On the dynamics of ball bearings[END_REF]

                                                                                              (3. 33)
In our study, strategy 2 is performed in advance, it is found to be very time-wasting, and given different conditions, iterative step should be changed frequently. Also, whether such a direct strategy will work is uncertain, because not all situations are considered to prove the correctness of the preceding exponential function relation. ( , ) 0

f Q P  ( , ) 1 '( , ) f Q P PP ii f Q P   ( , ) ( , ) '( , ) f Q P f Q P Q f Q P P Q P        Q P  
, as main work, should be derived based on bearing mechanics relation and alternative time-frequency relation of P. Considering that the complexity will increase with the numbers of harmonics used and degrees of freedom involved increasing, a uniform matrix relation is suggested in this chapter.

The acknowledgment of all Fourier coefficients can directly tell an arbitrary-time output using equations (3.30)~(3.32). However, reverse question should be answered differently. A simple idea below performs an inverse calculation of each hidden Fourier coefficient using time signals within a time period. [] 0 0

             22 1 2 [ ]cos( ) ( )cos( ) { ( cos(( ) sin ))cos( )} 2 0 0 00 N k ww in N F F t iwt dt C jwt D jwt iwt dt C C j j i i n N w n j               cos( 22
N CF t n N n     (3. 34) 1 2 2 [ ]cos( ) 0 N in CF t n i N N n      (3. 35) 1 2 2 [ ]sin( ) 0 N in DF t n i N N n      (3. 36)
Integration result can be analytical, for as each harmonic wave involved performs a multi-period integral.

               
Consequently, each sub matrix included in the overall Jacobian matrix is obtained.

0 0 0 0 0 0 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 0 1 1 0~0 1 1 0 1 1 C C C C C A A B A B KK C C C C C A A B A B KK D D D D D Q A A B A B KK K P C C C C C K K K K K A A B A B KK D D D D D K K K K K A A B A B KK                                                                     
Where,

NN F C t X in in i n K t n A N N X A N N nn NN F C t X in in in K jwt t n A N X A N N N jj nn N F C t X in in i n K jwt t n B N N X B N j j n                                 1 ) 0 N N n   
NN F D t X in in i n K t n A N N X A N N nn NN F D t X in in in K jwt t n A N X A N N N jj nn N F D t X in i n K jwt t n B N N X B N j j n                                     1 2 ) 0 N in N n    
Where,

2 n t wN   , 2 / Nt w   .
Some applications of frequency domain method in bearing system

Parametric influence on system stability

Even for the steady state responses, not all periodic responses included are stable, so that, some laws of stability have ever been studied, among which, Floquet theory is the most popular one.

As found in Part 3.5.1 (time method), initial states, damping, unbalance, and others strongly influence steady-state responses, it is very visible but needs post-treatment (like error, period identification, etc.) in quantitative analysis, and acquisition of more information also needs to burden expensive computations.

In the frequency domain method, direct output of system responses irrelevant to time bridges a more direct access to stability analysis, using the Floquet multipliers. Normally, monodromic matrix () At is extracted from state-space equation.

(

( ) () ( ) ( ) X t X t At X t X t              )
Here, due to complex nonlinearity, partial differential relation is used.

j NN n k k j kk t k C A t dt k k k t k N C kk k T At T N At Nj                             (3. 38)
Proposed frequency method makes direct use of time domain matrix ( ), C ( ) 22 K t t , so that matrix  is much more easily determined to search the instability regions by calculating its eigenvalues.

Under gravity, a response at y direction is obtained using 3D cylindrical bearing to support a rotor rotating at speed of 1000 rad/s. The corresponding Floquet multipliers can be calculated. It is found that the system is critically stable. Although now no constraint (preload) is set at x direction to increase stability of system, a slight disturbance still will not easily destroy the stability of system at speed of 1000 rad/s.

Speed is inevitably the main inspected parameter to control running conditions. Now, take the impact of speed on stability as an example. In equation (3.38),

f f P f v f f P v T P v P v T P P v                                              (3. 39)
Then, solution track follows a predictor-corrector order.

( , ) 1 1 ( ) ( ) 11
ff f P v ii PPPv ii T Pv T vv Pv ii P P v v ii                                             (3.42)
Where, Referring to [80], [START_REF] Zhang | Harmonic balance method with alternating frequency/time domain technique for nonlinear dynamical system with fractional exponential‖[END_REF], the tangential vector [ , ] Pv    can be calculated as follows, using equations (3. 39).

    ( 1) ( , ) det[ ( , )] 1 ( 1) ( , ) det[ ( , )] ( 1)
P i C P v DH P v i v N C P v DH P v N                     
Deleting i th column from ( , ) [ , ] ff DH P v Pv    ,   ( , ) ( ( , )(:,1, 1), ( , )(:, 1, 1))

DH P v DH P v i DH P v i N i      and Chapter 3   1 ( 1) sgn det[ ] 0 ( , ) 2 1 (det[ ( , )]) 1 f N P C P v N DH P v i i                
2. Around turning regions, even singularity of Jacobian matrix is mathematically avoided; it is still physically near singularity and fails to globally converge towards local minima. We think, it has a solution of same relative error ( / PP i actual



), but of bigger absolute error ( P i  ), the convergence turns bad, especially when a lot of variables are involved. To improve convergence rate and calculation stability when approaching its regime, the criteria used in convergence of Newton iterative method and selection of step size should be understood.

 Either absolute and relative accuracies of independent variable or absolute accuracies of independent and dependent variables should be satisfied.

1 ( ) ( ) 1 q q q q J f q i i i i i i         , Residual ( , ) ( ) 1
q i q v J q f q i i i i q i        
 On the other hand, iterative times counted at last integration step can be also used to consider whether reduced step size should be used or even iterative process should be stopped.

 Unfortunately, either accuracy control or step-size reduction still doesn't work. The solution asks for more efficient method of linear search technique for Newton's Method, developed in [83] and performed in [START_REF] Guo | Dynamic Analysis of Wind Turbine Planetary Gears Using an Extended Harmonic Balance Approach‖[END_REF].

If it is necessary to locate precisely the dangerous points, supplementary equations are sometimes needed to locate system singularity.

Results

Results of using 3D cylindrical bearing to support a rotor are plotted at a bigger range of speeds.

Figure3.17 Multi-Harmonic waves appearing at y direction It can be found from Figure3. [START_REF] Sarangi | Stiffness and damping characteristics of lubricated ball bearings. Part 1: theoretical formulation‖[END_REF] that the stability analysis is more difficult for 3D

condition and shows a more irregular distribution. Since the model used is cylindrical bearing-rotor model and no horizontal preload is imposed, all situations present slight instability, but some tendencies can still imply that the ranges of 380m/s~600m/s and 910m/s~1000m/s look much tenderer.

Further works on bifurcation, like the factors of bifurcation, period doubling, hysteresis, and the routes to chaos, can be found in the book [79].

Parametric influence on system responses

System responses include vibration amplitude, frequency as well as balance position etc. at output stage. Using frequency domain method, some information hidden in responses caused by internal and external excitations can be discovered. Here, influences of both internal and external excitations are briefly discussed.

 Internal excitation is still caused by varying rotating speed.

 External excitation acting on bearing-rotor dynamic system is described in an alternative but simpler manner.

( This is conveniently analyzed not only in the programming manner (in Chapter 3), but also in the experimental manner (It appeals to the exploitation of First generation bearing test bench proposed in Chapter 4). In the latter case, the desired structural parameters or dynamic relations in nonlinear bearing-rotor system can all be found.

( ) ( ) ( ) ( ( ), ( )) MX t CX t KX t f t f X t X t b rotor      ( ) ( ) ( ) X t X t )

Results

Since each composition in alternative external load is The balance position is influenced by both internal and external excitations.

For internal excitation, the change of balance position is the main reason that dominates bearing stiffness and influences rotor resonant regions and unstable regions.

If preload is largely imposed, resonant point can be predicted, using linear vibration theory. Even resonant frequency is also found to be time-varying, but larger the preload is applied, more accurate the prediction is. For same external excitation, speed rise will decrease vibrating amplitude in Figure3.26.

For same exciting frequency, different amplitude used has a similar distribution of resonant regions; the difference is that unstable area is enlarged at bigger amplitude (Blue dash lines in Figure3.20 and Figure3.21). For same exciting amplitude, higher frequency imposed benefits to stabilize balance position (Green dash lines in Figure3. 21 and Figure3.22). This once again evidences that the anti-correlation exists between system stability and system vibrating intensity. There should be a more accurate analysis on resonant regions and system bifurcation, but the accessible speed scope incorporated still needs to be increased by a more advanced and efficient mathematical solver. Currently, arc-length continuation method, accuracy control, step-size control, even simple linear search technique for Newton's method are used, but some problems of convergence in 3D question still occur, with more different singularity phenomena taking place at certain speed for more degrees of freedom, even they are proved to work very well in 2D questions.

Conclusion

Proposed bearing model of uniform format is engineering-convenient. It will therefore be possible to do such calculations to predict easily the performances of bearing and its adjacent components by implementing in any non-linear F.E.A package or any complex M.B.S package.

Both the new time-domain solver and the uniform frequency-domain method HB-AFT for calculating nonlinear bearing-rotor system can work to find some hidden information and to guide bearing engineering design and application.

Nomenclature

' cos 0 R r R i i i    - ' ( )cos 0 R D r R o o i     - ( )/
                                                                                   ) cos ( )sin ]cos } 0 {( )sin [( ) cos
                        
Where, sin sin cos 00

       cos cos sin 00        tan ( ) cos
                           0 0 0 N Q F x x N Q F y y N Q F z z                                ( ' tan sin ) 0
( ' tan sin ) sin cos 0 0 ' sin sin sin sin cos 0 00 ' sin ( sin sin ( 0 00 sin cos ) cos sin ) 0 0 0 0

(' 0 M R r Q y i i z N Rr ii F i FF ii FF i NN Rr ii NN Rr ii N i Rr ii                                   cos ) sin cos sin ) 0 0 0 0 2
FF ii FF ii FF ii N Rr ii NN R r R r r i i i i i NN R R r i i i                                          Similarly, ( ' tan sin ) 0 ( ' tan sin )cos cos 0 0 cos cos ( sin cos ) 00 00 M R r Q z i i y N Rr ii NN R R r i i i F i FF ii                              2017.09.18
Chapter 4

An alternative bearing test bench reproducing contact loads 4.1 Research challenge for vibrating bearing -Powerful bench

As running speed increases, the dynamic interaction between environment and transmission system will deteriorate bearings incorporated in rotor system. These bearings will experience a long-term dynamic external loading. Once they are broken, the system will lose stability and safety. Normally little different dynamic response outside bearings will be aroused only due to initial bearing failures, since other suspension elements installed outside bearing box used to isolate some influences of bearing vibration on other components. Thus, it is more important to detect bearing states than others at most times, to avoid further bearing degradation and to reduce financial cost. As a practical strategy, a powerful bench is also required either to verify computer-aided theory calculation or to better guide relevant structure to escape from some intense vibration regions.

Some resonance mechanism can be tested, using the first generation bearing test bench proposed in this chapter. However, bearing failures are always accelerated due to more intense relative vibrations and larger contact loads between races. Now, it is more elusive to make clear of failure mechanism: resonance, abnormal failure, or other fatigue failures with different conditions mixed. How to use the bench to reproduce them has powerful meaning.

However, searching from most of current techniques, only simulation and entire vehicle test can replay extreme conditions. But reliability and cost are still what we are concerned about.

To deal with it, an alternative bearing test bench is proposed as the second generation one, which directly reproduces the resonant results or other similar vibration conditions of races, but does not ensure the environments' reproduction of rotor or bearing box. The fatigue scenarios under different vibrations inside a bearing are expected to be discovered before reaching its life end or failing to work properly.

Due to

 All simulation results are not compatible to be reproduced by test bench. Some single stable solution can be found, but some double-stable solutions or unstable solutions cannot be found.

 Elements like rollers, cage are more easily ruined, so that if the dynamic behaviors of axle box and rotor are not very significant, in study, we can skip reproducing their environments, but directly to copy the dynamic bearing contact loads (or reproduce 2017.09.18 absolute motions of races), and then observe bearing non-dynamic parameters' change.

 Using the first-generation bench, the maximum dynamic load scope the vibrator can reach and the max motion acceleration reproduced in races are constrained due to neighboring rotor and supporting frame synchronically vibrating and consuming energy.

Step 1: Since external loads acting on train bearing are very complex and include a wide range of frequencies and amplitudes, which makes it of considerable impact to reproduce dynamic vibrations of the train axle box bearings to investigate failure mechanism under different vibration conditions. The bearing test bench proposed in this chapter meets this industrial demand. The maximum mass that the electromagnetic vibrator can support is 20 Tons, meaning that the electromagnetic vibrator can generate a maximum thrust force of 200,000 N, which will be responsible for all static and dynamic forces occurring above vibrator. Given a maximum vibrator power, the mass above the vibrator will constrain the max acceleration (frequency or amplitude), so that there should be a compromise among the platform mass and max acceleration. Given a minimum platform mass, a maximum acceleration must be strictly designed, considering other static and dynamic loads generated by moving coil and expansion

First generation bearing test bench

2017.09.18

The variable mass is connected rigidly to the tested bearing, while it is elastically connected to the ground with pre-compressed spring used to ensure an approximate range of static force.

Since high vibration acceleration and low vibration amplitude are performed in test, the dynamic bearing contact load is approximated.

Assuming that: Bearing contact force It can be found that both dynamic bearing force and motions of rolling elements inside bearing can be well copied in the alternative method.

Input signals of train bearing bench test

Bearing test bench acquires motions from vibrator.

Development on 2 nd bearing test bench conception

The bench test elaborated in Part 4.3 is reachable theoretically, but there are two technical problems. First, the spring used to load static force will have a slight deformation during bearing vibration. Second, supporting point in the lever moves on the raceway of selfaligning bearing, bringing error to leverage amplification ratio. To solve them, an optimized design is used with sliding mode control method adopted. along with rotating shaft, this design deliberately settles a spring of small stiffness at the bottom of left part, which matters little to system load, but helps to guide frame motion.

In this method, it results in that the load transmitted from the left to the right will be further reduced, since bearing static loads acting on both bearings at the left side are almost canceled. It means that more energy will be saved and higher vibrator ratio can be obtained.

Following the topic, the error caused by supporting point is also mathematically explained.

The dynamic bearing contact force is expected, The best approximation a ˆ meets continuous control law 0 s . In operation, an appropriate k is selected in real time by detecting unknowns required in the equation. a incorporating k is in fine input into vibrator controller, using sensor-feedback system.

The test system is asymptotically stable owing to its quick response to simple signals and its robustness under conditions of unstable system parameters and external disturbances.

However, the drawbacks of sliding mode control are still this asymptote or slight fluctuation of stability, as well as the necessity of understanding boundaries of uncertain parameters. 2017.09.18

Conclusion

This chapter proposes an alternative bearing test bench, from the first generation bench to the second generation bench, due to the increasing demands on high speed and heavy load, This thesis also explains in Chapter 3:

How to develop more practical bearing matrices and more reliable methods of solving nonlinear bearing-rotor model in a uniform and economic way? Whether frequency-domain method performs a better computer-aided identification capacity than time-domain method?

A uniform derivation for overall stiffness and damping matrices of both roller and ball bearings is conducted, using the pioneer work of Luc-Houpert (1997). Slicing technique and varying contact angle respectively considered in roller and ball bearings can lead to a more accurate prediction. Bearing matrices formulated with simplified EHL and Hertz relations can be analytical, however, the contact descriptions incorporated are still in doubt.

Nonlinear bearing-rotor system is modeled with a better linear mathematical equation and solved in simpler enhanced numerical manner of Adams. Meanwhile, nonlinear bearingrotor system is also remodeled and resolved with analytical method of HB-AFT. It performs a faster efficiency and an easier system behavior-identifying capacity, especially when considering 3D condition. After modeling and solving nonlinear bearing-rotor system, some factors participated can be well considered qualitatively and quantitatively, including initial states, damping, and eccentricity, internal and external excitations.

This chapter in fine provides a good interface between raceways, which can be involved in any non-linear F.E.A package or any complex M.B.S package. The solving methods can also be extended and then applied in more complex bearing-rotor system.  Linear superposition of the lubricant film stiffness, damping forces and Hertz force is a bolder hypothesis, since distributed lubricant force nonlinearly reconstructs Hertz deformation and force distributions. So that it will be more persuasive to perform that in other more complicated manners.

This thesis at last introduces in

 The proposed method for line EHL problem can be applied in Point EHL problem.

 The transition from line contact to point contact can be well considered when bearing is applied in rotor system. The selection of slice number can be further confirmed by comparing numerical bearing moment with experimental bearing moment.

 More information hidden in bearing-rotor system can be studied by considering edge effect, using linear vibration theory or non-linear rotor theory.

 The proposed slice load can be further introduced into bearing damping and stiffness matrices before performing dynamic analysis.

 More complex bearing-rotor system can be modeled using TTM or FEM. And a versatile code can be developed to solve TTM model of complex bearing-rotor system, using frequency domain method.

 The corresponding experimental results will be compared with theoretical solutions, with the construction of bearing test bench.

| Zhu S.C. First, perturbation method is extended to prove that steady-state EHL distributed damping and stiffness do not satisfy the transient conditions, and steady-state EHL lumped coefficients are not suitable for dynamic analysis.

Distribution law of EHL oil film damping and stiffness is attemptly proposed. Combining initial state of steady-state solution with transient EHL distributed stiffness and damping, dynamic behaviors of oil film are predicted.

Second, edge effect of roller/race contact is considered in bearing transmissions. First and second approximation methods are proposed to numerically interpret this effect. First approximation method concentrates on adopting influence coefficient method to replay edge effect as well as using normalization process to respect empirical results, while second approximation method focuses on a fast-analytical manner of describing Hertz pressure inside roller slice.

Third, practical bearing matrices are used to develop 3D bearing-rotor dynamics code, based on classical Hertz and EHL formulas. In time domain, four-order Adams method is improved to solve incremental nonlinear dynamics equations;

In frequency domain, analytical method of HB-AFT is particularly employed to bridge 3D bearing matrices and 3D steadystate solutions, besides, arc-length continuation method increases calculation speed. Using both methods, some factors can be considered qualitatively and quantitatively, including initial states, damping, eccentricity, internal and external excitations.

Last, an alternative technique is proposed to reproduce contact behaviors inside bearings under heavy loading and high-frequency vibrations, when working parameters are beyond bench capacity.

  Ph.D. work is conducted in both TPL (Traction Power National Laboratory) of Southwest Jiaotong University and LaMCoS (Laboratoire de Mé canique des Contacts et des Structures) of INSA-Lyon, CNRS UMR5259, F69621, directed by Professor Weihua ZHANG and Professor Daniel NELIAS.

le chapitre 1 et 2 :

 2 Comment les coefficients d'amortissement et de rigidité distribués de film d'huile correspondent à l'épaisseur distribuée du film d'huile sous des états stables et transitoires? (Chapitre 1) Si l'effet du contact en bordure entre le cylindre et le courant affecte les transmissions de charge et de motion? Comment interpré ter numé riquement cet effet de maniè re plus simple ? (Chapitre 2) Les calculs complexes ré alisé s dans les chapitres 1 et 2 essaient de dé crire les paramè tres distribué s de contact des cylindres en axes rotatif et longitudinal, et d'expliquer les caracté ristiques ré sultantes des vibrations du rotor, etc. Cependant, dû aux coûts des calculs, cette complexité est davantage pré fé ré e dans le calcul des é tats stables que dans celui des é tats transitoires. Sous des é tats stables plutôt que dans les é tats transitoires, les caracté ristiques Zhu S.C. | Pré face Contribution to the modeling of rolling element bearingrotor system for railway application 2017.09.18 dominantes de vibration du systè me devraient ê tre plus né cessairement distingué es et ê tre plus pré cisé ment compré hensibles.

  forces et moments du roulement peuvent ê tre obtenus d'une maniè re plus rapide et plus pré cise, ce qui Introduction | Zhu S.C. requiert la programmation d'un vibrateur é lectromagné tique qui soutient le systè me entier.

[ 17 ]

 17 . He proposed a new force model to apply in oscillatory system with multiple EHL contacts by interpreting each contact as hydrodynamic force in series with Hertz force. In his study, hydrodynamic Chapter 1 force was reduced to a complex function with regard to the cavitation position and its height, so that he reported another study[START_REF] Wiegert | An analytical expression of the normal force of hydrodynamic line contacts under transient conditions[END_REF] to get an approximate relation for the normal force of hydrodynamic line contacts under transient conditions by means of numerical results;

(Figure1. 1 :

 1 Figure1. 1: Transient oil film thickness distributions generated by steady-state solution

Figure1. 2 :

 2 Figure1. 2: Reduced algorithm for steady-state EHL solution The same advantage of proposed contact stiffness and damping distributions in solving transient EHL problem is also discovered.

Figure1. 3 :

 3 Figure1. 3: Reduced algorithm for transient EHL solution

Case 1 :Figure1. 4 :Case 2 : 18 Figure1. 5 :Figure1. 5 

 1421855 Figure1. 4: Comparison of steady-state and transient EHL contact pressures and oil film thicknesses During rapid acceleration and deceleration, the film behaviors are governed by two different mechanisms-fluid entrainment and film squeeze.

However, according to

  transient phenomena of oil film tested under ball longitudinal oscillation (Figure1. 6), Romeo P. Glovnea and Hugh A. Spikes [23] gave the same reason that during rapid deceleration, the formation of fluid entrapment dominated the film thickness over the contact; meanwhile, the depth and shape of this entrapment is dependent of the rate of speed variation. Zero halting time assumed in this chapter is an extreme situation, but also tells a same significant truth of existence of fluid entrapment.

Figure1. 6 : 18 Figure1. 7 :

 6187 Figure1. 6: Film thickness profile during rapid deceleration in longitudinal oscillatory motion ( Romeo P. Glovnea and Hugh A. Spikes, 2005)

[ 26 ]

 26 will be developed from a new perspective. Involved in this alternative technique, the fundamental load-deflection relation is particularly important. Hertz gave a width relation in line contact, but did not publish its deflection relation. The unknown issue was later discussed over a century. The relevant results are listed in Appendix C. N.Ahmadi and M.R. Hoeprich used numerical method,

3 .

 3 Direct bearing model is used to calculate loads (three forces F x ,F y ,F z and two moments M y ,M z ) with five relative displacements acknowledged (three linear displacements d x ,d y ,d z and two angular displacements dӨ y ,dӨ z ); 4. Bearing-rotor vibrating characteristics are predicted by bearing stiffness matrix calculated by finite differential method [44].

  azimuth angle due to deformation recovery of cup edge.

Figure2. 3 Figure2. 4 Figure 2 .

 342 Figure2.3

13 )

 13 General codes are also developed for double-row tapered roller bearing. Due to symmetry of the rollers at double rows, integrated model can be decomposed into uniform sub-models which accept inputs and transmit loads in different manners. Both coordinates transformation and load reconfiguration used improve general calculation in input and output stages.

Figure2. 6

 6 Figure2. 6 Load transmissions in double row tapered roller bearing

Mz

  est accessible à l'adresse : http://theses.insa-lyon.fr/publication/2017LYSEI136/these.pdf © [S. Zhu], [2017], INSA Lyon, tous droits réservés According to equation (2.1), line displacements of each row center caused by five displacements of entire bearing are given by,

,

   (in  ) for right row.Angular displacements of each row are also revised according to relative change of right and left coordinate systems.

dx=0. 8 /

 8 1000;dy=0.8/1000;dz=0.8/1000; w=5/180×π; dӨz =0.1/180×π; dӨy =0.1/180×π;

Figure2. 7

 7 Figure2.7 Load results in rotation

Figure2. 8

 8 Figure2.8 Load results in rotation 2.3.2 Alternative Roller Contact Calculation s (First Approximation Model) Uniform Harris' relation used to be applied in independent slice calculation. Relevant deflection relations in Appendix C are also useful either empirically or analytically.

Figure2. 9 Figure2. 9

 99 Figure2. 9 Load intensity of one roller with different slice number used

Figure2. 10

 10 Figure2. 10 Influence of slice number on roller calculation

Figure2. 11 2017

 11 Figure2. 11

Figure2. 12 Figure 2 .

 122 Figure2. 12 Figure 2.11~2.12 Slice division scheme (Bruno Mevel)

Figure2. 14 2 Figure2. 17 Figure 2 .

 142172 Figure2.14 Two strips, uneven strip height

  use, but more time will be consumed if more strips are chosen.

  2017.09.18it is still acceptable to have a less detailed description for the ones of small deformation, since loaded area is much smaller and error caused by transverse pressure distribution is smaller. So the way of involving different number of strips in each slice is given below.

  contains two geometrical unknowns (tilting angle t  and central deformation r  ).

  Figure2. 21 Geometric relations of Roller end/Flange Contact (D. Nelias)

  this case, the developed equilibrium loop model can recognize relative tilting motion of roller to cone as well as reciprocating motion of rib/roller contact position. Strictly speaking, tilting equilibrium of roller mainly depends on relative displacements between races. The resultant t  and r  will in turn slightly change overall distributed deformation   1 j N e   between each slice and interacted races in direct bearing model. However, the unique impact on additional deformation proposed in Paragraph 2.2 can be ignored.Discussion on Moment Accuracy of Second Approximation MethodIn calculation, compared to the forces (F x, F y, F z ), the moments (M y, M z ) are more easily affected by slice number. A more detailed picture depicts the relation between slice number and bearing loads.

  Figure2. 22 Slice number-bearing loads relation

  Figure2.23 COMPARISONS between original and enhanced moments

Figure2. 24 1 -

 1 Figure2.[START_REF] Reusner | The Logarithmic Roller Profile-the Key to Superior Performance of Cylindrical and Taper Roller Bearings[END_REF] 1-D Reduced interacting relation dependent of interacting distanceAccording to 1-D reduced interacting relation,

Figure2. 25 ( 2 . 53 )

 25253 Figure2.[START_REF] Harris | Rolling element analysis[END_REF] Application of 1D reduced contact algorithm in roller load calculationIt is found that wii is dangerous to calculation and needs a better approximation.

  Figure2. 26 Strip Scheme of replacing elliptically distributed pressure

. 56 )Figure2. 27 Figure2. 28

 562728 Figure2. 27 Distributed load intensity under different strip schemesIt is observed that two-strips scheme is sufficiently accurate to perform roller contact calculation.

Figure2. 29

 29 Figure2. 29 Comparisons between CST(red) and AST(blue) results

M M dx dx

  dy dz d y d z M dx dy dz d y d

dy=0. 8 /

 8 Figure2. 31 Comparisons between CST and AST stiffness with d x changing

Figure2. 33

 33 Figure2. 33 Double-row tapered roller bearing supporting a shaft at the centre Since angular displacements bring about more different nonlinear stiffness, it will be necessary to see how this difference evolves under more varying angular displacements.Below give two waterfall plots of preload(dx)-motion(dӨy)-resonant frequencies and one plot of rotation(t)-motion(dӨy)-resonant frequencies.

Figure 2 .

 2 Figure 2.34~2.36: dx=0.8/1000;dy=0.8/1000; dz=0.8/1000; dӨz=0.1/180×π;,N=31

Figure2. 34 Figure2. 35 Figure 2 . 2 Figure2.

 343522 Figure2.[START_REF] Palmgren | Grundlagen der Walzlagertechnik‖[END_REF]

[ 52 ]

 52 , Harris [53], Palmgren [54] and others, where either ideal boundary assumption or ignorance of some certain degree of freedom was employed, and failed to simulate the coupled translational and rotational vibrations. Bearing transmission load is obtained based on a discrete summation method ever used by Lim and Singh [55]~[57] or based on an integration method adopted by Hernot [58]before formulated into stiffness matrix. Houpert [59] provided a uniform engineering model for both tapered roller and ball bearings in 1997. Later in recent years [60]~[61], a transition from point contact to line contact as load increases has been considered, it chose displacements as model inputs, and a finite differential method was also referred in his work to get stiffness matrix. Like his interest, the authors of this thesis also developed a fast method of calculating bearing stiffness matrix considering roller edge effect in the work [59]. However, possible difficulties in finite difference formulation are the selection of differential step size as well as the need of replicated load calculation for several times. Besides these efforts above, what cannot be also forgotten is a FEM (Finite element method) model used by Guo and parler [62] and an experimental approach used by Knaapen 2017.09.18 R [63]. Although much would be cost, they still well stand for current engineering conceptions. … to bearing damping matrix On the other hand, bearing applications often occur in high-speed machineries vibrating at high frequency. Strict application conditions also lead to accurate analysis in bearing damping matrix. Sarangi, M. et al. [64]~[66] have proposed two models -Linear spring-damper model of Lubricated Contacts‖ and -Nonlinear spring-damper model of Lubricated Contacts‖ for both ball bearing in 2004, 2005 and for roller bearing in 2013, but they were for 2D questions, also based on their former works, in which, they got curve-fitting results by a transition from distributed EHL parameters to lumped parameters, but it was still founded on steady EHL. The authors of this chapter also show a different understanding of transient EHL contact stiffness and damping distributions in the first chapter.

Figure3. 2

 2 Figure3. 2 Equlivalent stiffness and damping of one rolling element Equivalent contact parameters for one rolling element are derived here by ignoring roller's mass and assuming linear vibration.

Figure3. 4 ,

 4 Figure3.4, initial and general relative positions of inner and outer races are shown, and the rollers occupied in loading area undergo elastic deformation.

Figure3. 3

 3 Figure3.3 Deformation distribution in roller bearing(Houpert) Euler angles are used to demonstrate the relative rotations between raceways. The rotations around z axis and updated y axis are considered successively. Tilting angles of the cone relative to the cup cause additional displacements at any spatial point.



  In roller bearing, slicing technique used meets with uneven distributions of EHL and Hertz behaviors in each misaligned roller. In ball bearing, variation of contact angle described agrees with essential fluctuations of deformation and load transmission. Bearing transmission matrices are in fine derived with five relative displacements ( ) between raceways acknowledged.

2017.09. 18 …

 18 Using incremental linear equation to preform nonlinear relation, it will be more persuasive in application. Differential format (3. 23), (3. 24) of equation (3. 21) is employed.

.

Figure3. 5 A

 5 Figure3.5 A symmetric rigid rotor-nonlinear bearing System

2 C

 2 y z stand for preloads or varying external excitations acting on the rotor center equivalently; and e represents eccentricity value. Parameter-varying bearing coefficients 2 K are derived in Part 3.4. Below, two novel solvers are respectively developed for both time and frequency domains analysis.

Figure3. 8

 8 Figure3.8 Responses of balanced rotor starting at Zero point

Figure3. 9

 9 Figure3.9 Responses of balanced rotor starting at Non-zero point

(

3 InFigure3. 9 ,,

 39 different vibration scope and different balance point are found, it is because the setting of ( ), ( ), ( ) X t X t X t at = contains more information than initial transferable potential energy, preloads hidden inside represent untransferable potential energy and constrain the motions at = , different initial untransferable (preload) or transferable potential energy (the extent of compressing bearing -spring‖) result in different vibration frequencies.

Figure3. 11 (

 11 Figure3.11 Responses of balanced rotor starting at Zero point

Figure3. 13 ,

 13 Figure3.13, vibration amplitude is reduced compared to Figure3.10, using constant damping values. In Figure3.14, VC (varying compliance) vibration is nearly totally suppressed, when adopting simplified steady-state damping coefficients in an integrated manner.

  materials and lubrication. More accurate damping coefficients expect more work. At the stage of present knowledge, a most complex lubricating damping function is the curve-fitted relation developed by Chippa and Sarangi [64]. The inlet force prying up the rollers will experience fluctuations, once contact vibration happens. This fluctuation forms an additional damping force. It is admitted that formulations proposed by Chippa and Sarangi are numerically excellent, but  The perturbation method used is based on a synchronous perturbed value of film thickness along contact width. Synchronous perturbation of film thickness distribution, in essence, trusts robustness of steady-state EHL under transient conditions. A consideration of asynchronous film thickness variation along contact width is taken by the authors in the first chapter.  Linear superposition of the lubricant film stiffness, damping forces and Hertz force is still a bolder hypothesis, since distributed lubricant force nonlinearly reconstructs Hertz deformation and force distributions.

  recent literatures, systematic dynamics relation and bearing mechanics relation are mixed together, where assumed preceding exponential function relation is replaced by a new one using HB-AFT.



  So that, C i and D i are described using[] F t n .2017.09.18 1 1

4 10 , 5 Nn

 5  is assumed. Chapter 3 Normally, T VC -1 solutions can be searched between 1 w vc and 2 w vc , therein, 10 rad/s as a stepsize, afterwards, 1 rad/s as a finer step-size around unstable region. However, a full-parameter tracking simulation as demanded above encounters a robustness question of frequency-domain method. 1. Around unstable regions, multi-solutions make initial iterative position and iterative step-size more elusive compared to general regions (even for single solution, it is also hard to locate initial position around accurate solution). To deal with the problem, especially around the turning points, solution experiences a dramatic change. The arc-length continuation method which delegates a continuous version reformulating from the homotopy [79] can be well embedded.

Figure3. 16

 16 Figure3. 16 Arc-length continuation method

  Turning points 2017.09.18 At predictor stage, it is calculated on basis of an accurate neighboring point ( , ) 00 Pv . The direction perpendicular to tangential direction of previous solution determines next initial guess. The curvature of solution line also determines the size of iterative step [80].

( 3 . 41 )

 341 At corrector stage, the global Newton-Raphson iterative method is employed to calculate the extended system (3.42), where, both 0 f  and 0 g  should be satisfied. In the global Jacobian matrix, f v   can be approximated by finite difference method, a way different from that performed in [81], where other forces (bearing forces partly dependent of rotating speed) are not taken into consideration except inertial forces.

1

 1

  accurately perform next initial guess and avoid singularity of Jacobian matrix.

  which has been well recognized in [81][85][86] as a parallel equation, and also a replacement of 0 g  combined with 0 f  to directly identify bifurcation points.2017.09.18Since  is defined as the eigenvector corresponding to eigenvalue an initial iterative vector of bifurcation point-identifying algorithm.For a good parameter-varying track, det f P to inform algorithm to use reduced step size.

Figure3. 18 3 Figure3. 19

 18319 Figure3.18 Multi-Harmonic waves appearing at z direction

  Figure3.23 Multi-Harmonic waves appearing at z direction 1, 0.001 ka k 



  For high-amplitude bearing resonance resulted from external excitation, normally no prior estimation for approximate regions can be made, since balance point takes less effect on entire bearing stiffness. Only full track simulation can help to predict the resonant regions.

  Figure3.26.

  effective radius of curvature in the principal x-plane Cette thèse est accessible à l'adresse : http://theses.insa-lyon.fr/publication/2017LYSEI136/these.pdf © [S. Zhu], [2017], INSA Lyon, tous droits réservés Equivalent radius in rolling direction ( for principal x or y plane or sum radius)

  Different bearing responses caused by different external loads do send some hidden alarming signals of bearing life-span and update period. These warnings may come from (a) experimental descriptions of bearings, like shorten time of temperature increase, haste of relative vibration intensities between races, louder noise, and so on (they are useful) or (b) from analytical prediction of undetectable element (when relative motions cannot be detected directly, absolute motion of outer race can be detected instead), since rotor motion can be estimated by varying vibration intensities of bearing box. For detectable element, not all lower vibrations can ensure no failure. E.g. zero vibration may represent infinite bearing stiffness, which means that some part associated to rotation is broken and cannot move; Moreover, not all higher vibrations necessarily trigger off higher possibility of failure, it mainly depends on the factors arousing relative motions.On one hand, bearing entire vibration does relate to bearing failure. Relative vibrations between races that, the thesis thinks, have direct relation with failures, should be caught at different degrees of freedom and under different excitations.On the other hand, complex failure conditions always mix together, but in general, cage deformation, roller or race materials loss, and radial diameter expansion of raceway due to temperature rise, etc. all influence dynamic relations in original design. It sounds clear to distinguish them with design results.Chapter 4As a modern strategy, using a nonlinear mathematics relation to characterize it is the most economical way to promote future design, perform fault diagnosis and maintenance.One major manner in which bearing feels the vibrations is so called resonance, it points out that the balance of rotor is broken, vibration mode reset, and relative vibration intensities changed once different external excitations or failure excitations are imposed. The contact force between rotor and housing dominated by their relative displacement will change due to absolute housing motions or failure excitations introducing rotor inertial forces (3. 43).



  Outer race motion In the alternative bearing bench, effective dynamic loads are recreated, and the motions of tested bearing are accelerated, while ineffective dynamic loads are reduced. Test benches from different bearing manufacturers [89][90][91]differ in main structure, loading technique, testing technique and control technology. Test benches in this chapter are enhanced alternatively in loading technique and control technique. With tests conducted on the electromagnetic vibrator, this scheme can not only reproduce heavy loads, but it can also reproduce high-frequency vibration, which is rarely involved in bearing researches in the history and relevant for train bearing tests.

  As far as train bearing is concerned, using multi-body dynamic software SIMPACK, loads and motions of axle box bearings in the train model can be differently predicted according to different bearing model participated. Although normally bearing is rigidly modeled, it still has little effect on predicting dynamic bearing contact loads, due to low correspondingly unchanged dynamic loads of neighboring components (except axle box). However, the evaluation on dynamic bearing motions would be more complicated. Only precise bearing model used can find it.Using online test, both race motions and contact loads of bearing can be acquired. The motions of components in a bearing differ from each other. Motions of the weakest elements, like roller, cage, etc. can be practiced by setting similar motion boundaries of bearing in test bench with boundary motions set under specific frequency and amplitude.

Figure4. 4

 4 Figure4.4 Second Generation Bearing Test BenchIn order to ensure no relative motion between the downside of loading spring and housing, this chapter also uses the design in Figure4.4. Another bearing is added on the rotating shaft so that two bearings on the left side interact with each other, using same radial static spring force. Meanwhile, to ensure bearings on the left side to synchronically work

F

  . In improved design, a simple variable structure control (sliding control[START_REF] Jean- | Applied Nonlinear Control[END_REF],[START_REF] Liu Jinkun | Research and development on theory and algorithms of sliding mode control[END_REF]) is performed to predetermine a -sliding mode‖.Since system structure parameter changes with state changing, the system has to input a differently in accordance with the state trajectory of the predetermined -sliding mode‖. The uncertainty of system caused by b is reduced below. Bearing acceleration is obtained by vibrator through the -lever‖. Considering sliding range of fulcrum on the raceway of selfaligning bearing, acceleration coefficient b varies within a certain range.To meet bearing acceleration value x  with the designed dy x  , the initial value should fit, interference (neglected), and b has a range of, i.e., a involved in governing equation is a multiplier, gain b can be estimated, using the geometric mean value of the above-mentioned boundaries. must remain on the sliding surface St ( )

  as well as inducing enhanced requirements on vibrator capacity and more expensive test cost, original test on nonlinear bearing-rotor dynamic behaviors is replaced by an alternative test on bearing contact behaviors by reproducing dynamic bearing contact load and dynamic roller motion of some main components, like resonant components. It does not reproduce rotor behaviors, but economically considers the endurance or failure questions of weakest contact bodies inside a bearing. Compared to current bearing tests, this method is advanced, since the conditions performed not only include reproduction of dynamic loads, but also reproduction of absolute motions of both races, even under heavy load situation. This reproduction can consider the impact of internal and external excitations on bearing contact load and roller motion or roller inertia load. Even although roller inertia load is low, roller motion does also impact its contact behaviors at unloaded areas and the ways of temperature rise, etc.In conclusion, to understand more contact information inside a bearing under different application situations, the resultant change needs to be discovered numerically and experimentally. Some simple results can be mentioned in Part 3.5; working parameters at a range of as wide as possible can be considered in a very simple manner in this bench. And meanwhile, fatigue or so-called longevity test can be also economically performed by accelerating vibrations.General conclusions and prospects | ZhuS.C.  An engineering-adequate way of bearing calculation by combining normalization process and advantages of both approximation models;  The influence of new algorithm on bearing stiffness components and relevant vibrating characteristics;This chapter in fine provides a better description of longitudinal contact between roller and raceways, which might be involved in any non-linear F.E.A package or any complex M.B.S package with a very simple code.

Chapter 4 :

 4 | Zhu S.C.  The perturbation method used is based on a synchronously perturbed value of film thickness along contact width. Synchronous perturbation of film thickness distribution, in essence, trusts robustness of steady-state EHL solution under transient conditions. A consideration of asynchronous film thickness variation along contact width is taken by the authors in the first chapter. But compared to distributed parameters, the corresponding lumped parameters are still unknown.

  'UNIVERSITE DE LYON OPEREE AU SEIN DE L'INSA LYON NOM : Shaocheng DATE de SOUTENANCE : le xx novembre 2017 Prénom: Zhu TITRE: Contribution to the modeling of rolling element bearingrotor system for railway application NATURE : Doctorat Numéro d'ordre : 2017LYSEI136 Ecole doctorale : MECANIQUE, ENERGETIQUE, GENIE CIVIL, ACOUSTIQUE Spécialité : Mécanique -Génie Mécanique -Génie Civil RESUME: Some developments on bearing applications are emphasized in this thesis, from contact model to bearing model, and from bearing-rotor model to bench reproduction technique.

Appendix A Perturbation Formulation for stiffness and damping distribution In

  Reynolds equation, main variables are perturbed with respect to its steady-state values, so the perturbed film thickness terms and corresponding pressure changes are given by,

							0 0 h h P P P h P h h hh            	(A.1)
	By regarding the steady-state solution of EHL as initial state, and continuously acquiring Where,
	EHL stiffness and damping distributions, dynamic mechanical behaviors of oil film are
	predicted. The method proposed herein in fine provides a possible insight to understand the , PP PP hh h h     (A.2)
	distributions of EHL oil film damping and stiffness, which can be extended for transient point
	EHL problem as a fast calculating technique. The idea proposed can be easily delivered to any After substituting the perturbed terms (A.1) into Reynolds equation:
	engineer.	x  	   	(	0 h	3 ) ( h d P P h P h 0 hh dx      	)	(   0 h h x   ) 12 12 u    	(   0 h h t 	)	(A.3)
	Where, terms containing	,( h h n m nm     are neglected. 1)
	According to what is referred in [19], three equations are acquired,
								3 0 h dP 0	12	u	0 h	for( , ) h h	0
							x	dx				x
								32 3 0 0 h dP h dP h 0 0for( ) 1 h
							x	dx				dx
								3 0 h dP h	1 12 0for( ) h
							x	dx
												2017.09.18

  is used to calculate it.

	Reynolds equation,	( ) f p		d dx	3 h dp dx  () 12 d h u dx      		0
	For compressible fluids,	( ) f p	12    dx dp	u	(1 	2 ee h h  	)		0
									h
	By Newton-Raphson method, it is changed to	( ) ( f p k 	p k	1 		p k	)	() df p k		0
											dp
	Dimensionless variables are made below,
			,			,		2 22 , 0	0	,	, 00	,
			,		,			,		3 2 24 ,	2

  It will be solved with load balance condition and boundary conditions.

	w	1 x e x  	pdx	or	1 x e x 	Pdx		2 		and	0 P x  , 1	P x e x  		0
	By setting iterative increments (	 	H e e	) ( , n	P 	j	, ) ( n	H 	0	)	n	,
				(		H	) ( n		H	) ( o	H	)	n
						e e					e e	e e
				n o P P	() P	n
													j
					jj
				n 00 o H H			() 0 H	n
	1 (unit, GPa )  (ln P z   0	9.67)	,	L P			0.196 GPa  , z 0.5 0.98
	L
														2017.09.18

2.1 Introduction

	Poor roller end design as well as misalignment in bearing-shaft motion gives rise to
	pressure concentration and potential damage. In practice, SRBs (Spherical roller bearing) or
	TRBs (Taper roller bearing) frequently accept multidirectional motions or loads, so that
	misalignment turns into one of the main failure factors. Some researches [24] point out that
	the logarithmic roller profile produces uniform pressure distribution along roller length under
	uniform loading, but regular tilting motions (uneven loading) still cause frequent pressure
	concentration. So it requires an accurate calculation on force and moment distributions in a
	bearing for predicting bearing stiffness matrix, which is essential to demonstrate how well
	bearing accepts misalignments in transmission system.
	How lamina load relates to local deflection determines roller position and load
	distribution in bearing equilibrium. Slicing technique is a popular way to describe complex
	load distribution in a bearing, but it is also understood that it fails to consider interaction
	among laminas. The question that arouses our interest is the choice between classical slicing
	technique [25] and essential contact algorithm [26]~[29], when it refers to bearing
	transmission analysis under heavy-loaded conditions, under frequent tilting motions, and
	others.
	Bruno Mevel [

30]~[31] recently

	Although CG-
	(2013,2016) addresses how the roller equilibrium is
	influenced by the calculation method. His work happens to meet with our interest in
	improving bearing stiffness. Since tilting motion takes responsibility for load concentration in
	roller end, this effect will be on significant components of bearing stiffness matrix, affecting
	prediction on stability of system when applied to dynamic system. More concerns about how
	to find a better solution in bearing transmission should be taken, both in calculation efficiency
	and in accuracy.
	Linear classical lamina load is applied along roller length, as Harris formula explained
	by Taylor expansion. This failure of missing edge effect requires better codes development in
	simulation.
	2017.09.18

mlms [32] or CG-FFT [33] methods

  greatly extend our insight into knowledge about contact problem including roughness, indents etc., most engineering problems still need time-saving solution. This spectacular performance is not necessary for bearing transmitting calculation. Interest is to find time-saving solution of existing theories.

	Based on Bousinesq work, second approximation method was typically performed by
	Reusner H.(1960) [27],

  dx=0.8/1000;dy=0.8/1000;dz=0.8/1000;dӨz=0/180×π;dӨy=0/180×π;

	b) dx=0.8/1000;dy=0.8/1000;dz=0.8/1000;dӨz=0.1/180×π;dӨy=0.1/180×π;
	c) dx=0.8/1000;dy=0. 8/1000;dz=0.8/1000;dӨz=0 /180×π;dӨy=0.4~0.7/180×π;
	d) dx=0.8/1000;dy=0.8/1000;dz=0.08/1000;dӨz=0.1/180×π;dӨy=0.1/180×π;
		-3.73	x 10 6	Fx		9.55	x 10 6	Fz				-1.605	x 10 5	My			-3.74	x 10 6		Fx	9.9	x 10 6	Fz	-1.54	x 10 5	My
		-3.74													-1.61						-3.75				9.88	-1.56
	N	-3.75				N		9.5						N*m	-1.615					N	-3.76				N	9.84 9.86	N*m	-1.6 -1.58
		-3.76													-1.62						-3.77				9.82	-1.62
		-3.77	0	0.1	0.2	9.45	0	0.1	0.2		-1.625	0	0.1		0.2		-3.78	0			0.1	0.2	9.8	0	0.1	0.2	-1.64	0	0.1	0.2
					t [s]						t [s]						t [s]							t [s]	t [s]	t [s]
		-4	x 10 6	Fx			-4.08 1.094	x 10 6 x 10 7	Fx Fz				-4.14 1.128 x 10 6 x 10 7 -1.65 x 10 5	Fx	Fz My		-4.22 1.162 x 10 6 x 10 7 -1.7 x 10 5	Fx	Fz My	1.198 -1.7 x 10 7 x 10 5	Fz My	-1.7	x 10 5	My
		-4.01 -4.02				N	-4.09 1.092 -4.1 1.09						-4.16 1.126 -1.7 1.124 -4.18 N -1.75				-4.26 -4.24 1.16 -1.75 1.158 N -1.8			N	1.194 1.196 -1.8	-1.8 -1.9
		-4.04 -4.03					-4.12 1.086 -4.11 1.088						-4.2 1.12 -1.8 1.122				-4.28 1.156 -1.85			1.192 -1.9	-2
		-4.05	0 0.05 0.1 0.15 0.2		0 0 1.084 -4.13	0.1 0.1		0.2 0.2		0 1.118 -4.22 0 0.05 0.1 0.15 0.2 0.1 0.2 0 0.1 0.2 -1.85	0 1.154 -4.3 -1.9 0	0	0.1	0.1 0.1	0.2 0.2 0.2	0 -2 1.19	0	0.1 0.1	0.2 0.2	-2.1	0	0.1	0.2
											t [s]						t [s]							t [s]	t [s] t [s]
											-2.72	x 10 6		Fx			1.365	x 10 6		Fz
																		1.36
											-2.73						1.355
											N
											-2.74						1.35
																		1.345
											-2.75	0	0.1	0.2		1.34	0	0.1		0.2
															t [s]						t [s]

a)

  With given amplitude and initial phase, a harmonic vibration can be described.

								K	H	1.5( 		'2 2 3 ) (  3 k E RF e 	1 3 )	1.5 	F		1 	For ball contact
			K	EHL		2 N 0.27733 ' ER y	0.7 0.6 U G W [ d Ri e 1.13 1 ]   []  7.69240 dt   	Fh c	1 	For roller slice contact
		K Where,	(t) EHL  0	2 cos( [ 0.18023 0 ' ER y U 2 2 0.67 0.53 1.067 ) [( i )e ] 0 1 2 0.73 (1 0.61e G W    1 2	k e	1  )]	14.9254 	Fh c	1 	For ball contact
	Based on Sommerfeld boundary conditions and simplified integration of Reynolds 2 arctan 0 1
	equation, loads [76], [77] supported by inlet oil for roller and ball contact are
								ef F			(	4 u 0 h 	R l x x		3  0	1.5 1.5 R l x	u ) z	1 N	For roller slice contact
												c	h
												c
								F ef			h	4 u 0 2 ( / x x R a 1) k k  		6  2 0 h R	1.5 1.5 x	a	u	z	For ball contact
												c	f c	c
	Central film thicknesses are given by
												h c	4 3  1.6 U G W 0.7 0.6	0.13 	For line contact
				h c		2.69 U	0.67 0.53 G W	0.067 	(1 0.61e 	0.73 	k e	)	For ball contact
	Steady squeeze effect at inlet area can pry up the roller, similarly, the transient squeeze
	effect caused by transient term in Reynolds' equation should also have some effects, resulted
	from u	z , but relevant to formation of h c .
	Equivalent damping coefficient
				cC ef 		d dF u ef		3  0	1.5 1.5 R l x	1 N	For roller slice contact
												z	h
												c
				ef c C 			u dF d ef z		1.5 1.5 a x R h 6 0 2 	For ball contact
												c
	Inlet squeeze velocity is nearly equal to approaching velocity of two bodies,
								0.9 0.8 2.1054 ( ') E l KF 0.1 ( ) 1.2037 N dt N F H   z u d  		1 	For roller slice contact
	Chapter 3										2017.09.18

for line contact, Brewe and Homrock's simplified solution [74] for point contact, Dowson and Higginson's line EHL formula as well as Homrock and Dowson's point EHL formula[75], dynamic contact parameters are derived as (SI units)

Overall stiffness and damping matrices

	sin	0
	cos cos	0				1
	[ ] [ { sin cos 1 sin sin	0 sin cos	[ 1	(	2	1	1)	]}]
	cos sin	cos cos
								2017.09.18

  are all set as zero, and vibration highest peak is located at zero point. It points out that if no gravity and preload are considered,

	8, initial states at = -2 t , the rotor is well balanced at = -2 -t t  , -0 t t t  , , and has zero initial transferable potential 0
	energy.

  

	Strategy 2:
	Where,	ww cage 	is the fundamental angular frequency.
	If external forces participate, the Fourier coefficients of all harmonic items involved
	should be written into the vectors 0 E (on behalf of preload),	E	j and	H	j respectively.
	Now, in order to search a constant vector P constituting periodic solution, how to solve
	such a group of nonlinear algebraic equations is of another interest in this chapter. Two direct
	iterative strategies are commonly used for general nonlinear algebraic question.
	Strategy 1:	P X F i bearing  	FFT 	() P FFT i Q	dynamic 	P i	1 
								2017.09.18

  (left one). While in strategy 2, iterative relation seems to have a convergence. But we doubt that: 1, whether it is right that the bearing mechanics relation also performs an

	exponential function relation with exponent greater than 1(near 1.1 or 1.5) in frequency
	domain; 2, each calculation in strategy 2 contains lots of time-consuming bearing static
	balance operations, where iterative step involved is very sensitive to different load
	combinations. The step-size control coefficient in (3. 33) is very elusive even although it just
	needs to be programmed once for all.
								Fx	Fx
	x			x		x
								Fy	Fy
	y			y		y
	z			z	( ( coeff inv K	z	))(	Fz	Fz		)
								MM
	y			y		y
								yy
	z	i	1	z	i	z	i	MM zz bearingi	bearingi	1

Appendix D Ball deformation and Contact angle

						x O	d	x	d		x	x	I	x O	x	I
		( )	A	( )	0 A	y O	d	y	d	y	y	I	y O	y	I
						z O	d	z	d	z	z	I	z O	z	I
	E	2 d x x I 22 12 u u u  {( x d O x  0 (1 )/ (1 )/ 12 )sin [( 12 EE     2 ... {( '  )sin [( 0 x x y I O I	y d O y	Entraining surface velocity, and 1,2 for two contact bodies ) cos ( )sin ]cos }... d y z d d z y I O z z I Lubricant viscosity in the inlet area Reduced Young's modulus
				l											roller length
				F											Normal applied load
				e											Rotor eccentricity
				m											Rotor mass
			, JJ											Rotor inertial moment
			dp

  First generation bearing test bench is composed of electromagnetic vibrator, moving coil, expansion table, and bench platform. Bench platform design is the main job in developing dynamic bearing test technology.

reproduce  Cette thèse est accessible à l'adresse : http://theses.insa-lyon.fr/publication/2017LYSEI136/these.pdf © [S. Zhu], [2017], INSA Lyon, tous droits réservés Chapter 4 Figure4. 1 Entire Test Job Figure4. 2 First Generation Bearing Test Bench

table ,

 ,

	and vice versa.
		M	moving coil		M	expansion table		M	bearing platform	 200000 a  	N	(4. 1)
	M	moving coil		M	expansion table	510 	kg M ,	bearing platform	500 	kg	(4. 2)

Cette thèse est accessible à l'adresse : http://theses.insa-lyon.fr/publication/2017LYSEI136/these.pdf © [S. Zhu], [2017], INSA Lyon, tous droits réservés

Kanglez/anglez Cette thèse est accessible à l'adresse : http://theses.insa-lyon.fr/publication/2017LYSEI136/these.pdf © [S. Zhu], [2017], INSA Lyon, tous droits réservés

Acknowledgment/Remerciements | Zhu S.C. and my French friend A. Domissy, who helped me a lot in French learning, and gave his precious help in correcting my written French of abstract. At last, I must acknowledge the financial support provided by China Scholarship Council (CSC), and the conferences funding from National natural fund Committee of China and State key laboratory of Traction power, SWJTU, China, as well as relevant courses funding given by LaMCoS, INSA Lyon. Many thanks. 2017.09.18

2017.09.18

However, the story does not end yet, since how the relationship is quantitatively established is still not clear. Interestingly, all other settings of coefficients before To eliminate this error inaccurate () Xt introduces to

governing equation is rebuilt.

It is obvious that the system stability is parameter-varying. To further perform it, the discrete state transition matrix developed by Hsu [78] is used.

To reproduce engineering facts, excitation harmonic load here can be explained as

Displacement excitation of outer race is in essence an inertial load excitation of rotor, which will interfere the balanced rotor and reset new balance. In such a case, the influence of external excitation can be easily mathematically investigated.

Appendix E Moment of Ball bearing with consideration of varying Contact angle

Due to sin sin cos 00 cos cos sin 00

Appendix F Four-order linear multistep numerical Method (Four-order Adams Method)

Using polynomial of third degree, these nonlinear variables are implicitly predicted,

Participating coefficients are accordingly calculated based on accelerations selected at continuous simulating points in a recursive manner.

Introducing relevant variables into governing equation, acceleration is in fine derived explicitly.

  The calculation shows that if the bench platform weighs over 500kg, the vibrator will generate an acceleration of less than 20g. In practice, the accessible or reachable acceleration of vibrator is very sensitive to mass, since acceleration decrease/increase performs a faster step than mass increase/decrease. Therefore, the test bench mass should be carefully designed.

In first generation bench design, light material and simplified structure are used if they matter little to test requirements.

Using first generation bearing test bench, rotor behaviors caused by nonlinear bearing can be more focused. E.g.,

 Amplitude-frequency curve can be detected to investigate the resonant regions. How different factor from either internal or external excitation arouse hysteresis phenomenon, jumping phenomenon, or stable and unstable regions?

 The resonance behaviors can be further analyzed by recognizing the harmonics involved can be using FFT technique.

Conception of 2 nd generation bearing test bench

Bearing test bench reproduces bearing forces and bearing motions. This equipment (Figure4. However, in application, these forces are shared by both effective and ineffective masses.

Considering a great amount of electricity needed to drive bearing test bench, ineffective mass or force shall be eliminated to improve usage of vibrator. In the first-generation bench, ineffective mass has been reduced as much as possible; now this part looks for another way.

Mark spring force as F 0 , vibration acceleration as a, mass of bearings, suspension boxes, and rotor as m, the effective force can be approximated: F dy = F 0 + ma.

Chapter 4 In application of train bearings, even for axle box, the acceleration tested on the scene reaches 20g~50g, far beyond the capacity of the present vibrator. Using proposed bench, important components of bearing contact force can be versatilely reproduced to perform alternative tests.

Below gives a more economic bench, especially aimed for heavy-bearing load reproduction case, where the power shared by ineffective masses is reduced.

General conclusions and prospects

This thesis explains in Chapter 1:

How the distributed damping and stiffness coefficients of oil film correspond to the distributed oil film thickness under steady and transient states?

Perturbation method is extended to prove that steady-state EHL contact damping and stiffness distributions should not be used to analyze the transient conditions, let alone acquire equivalent total contact damping and stiffness coefficients for dynamic analysis.

This thesis explains in Chapter 2:

Whether edge effect of roller/race contact affects load and motion transmissions? How to numerically interpret this effect in a simpler manner?

Bearing, as a key part of vibration system, participates in vibration transmission. A better prediction of bear transmission is introduced by taking account of pressure concentration.

To interpret pressure concentration easily, below are performed: This chapter proposes an alternative bearing test bench. From the first-generation bench to the second-generation bench, due to the increasing demands on high speed and heavy load, as well as increasing requirements on vibrator capacity and increasing test cost, original test on nonlinear bearing-rotor dynamic behaviors is replaced by an alternative test on bearing contact behaviors by reproducing dynamic bearing contact load and dynamic roller motion of some main components, like resonant components. This reproduction can consider the impact of internal and external excitations on bearing contact load and roller motion (or roller inertia load). Even although roller inertia load is small, roller motion does also impact its contact (impact) behaviors at unloaded areas and the ways of contact manner, temperature rise, and others at loaded areas.

To understand more contact information inside a bearing under different application situations, the differences can be mathematically or experimentally. Besides mathematical manner in Chapter 3, which is economic; the experimental method is also developed in Chapter 4, furthermore, working parameters at a range of as large as possible can be also replayed in a very simple manner in this bench. And meanwhile, fatigue or so-called longevity test can be also economically performed by accelerating vibrations.

Prospects:

From contact model to bearing model, and from bearing-rotor model to bench reproduction technique, the four chapters deliver the author's understandings in bearing engineering application. However, a lot of relevant ideas are still there and waiting for some clearer solutions.

 Complex calculations performed in Chapter 1 and Chapter 2 attempt to describe distributed roller contact parameters in rolling and roller length directions, and to explain some resultant vibration characteristics of rotor, and so on. However, due to calculation cost, a simple and lumped parameter is more preferred in dynamic calculation of multi-EHL contacts. But, the extraction of accurate Transient EHL parameters is still very difficult.