Nathanael Thibault

Adeline Bourdon

Didier Remond

Damien Lecouvreur

Dynamic model of a deep grooves ball bearing dedicated to the study of instantaneous angular speed of rotating assemblies

Keywords:

The present paper describes a ball bearing model for transmission monitoring and adapted to non-stationary operating conditions. The model considers axial forces in bearing as well as centrifugal forces applied on rolling elements making it suitable for high speed and axially loaded systems. A spall is also introduced to study its detectability. The model is based on an "angular" approach which preserves the complete main rotation of the shaft and is not limited to the torsional deflections. Results obtained by modeling a simple system under different operating conditions with healthy and damaged bearings are discussed to highlight the bearing behavior dependence on axial load and nominal speed. It is showed that the increase of nominal speed mainly affects the repartition of the power on the harmonics of the bearing signature while the increase of axial load leads to lower influence of the bearing on the Instantaneous Angular Speed (IAS). Simulations performed including spall defect on outer ring pointed out the impact of the damage on the IAS and its detectability.

Nomenclature

Mass matrix of the complete system Damping matrix of the complete system Stiffness matrix of the complete system Complete degrees of freedom vector External forces applied on the system Connecting forces generated by bearings , , Translation degrees of freedom of node i , , Rotation degrees of freedom of node i Angular position of ball j Angular position of cage Angular velocity of cage Angular velocity of inner ring , Axial and radial position of the center of curvature of inner race at ball j , Axial and radial position of the center of curvature of outer race at ball j , Axial and radial position of the center of ball j , Radius and diameter of balls , Inner and outer ring radii , Inner and outer race curvature radii Radial clearance , Normal contacting forces on inner and outer ring , Angle of contact on inner and outer ring Centrifugal force ! * Reduced Young modulus # Ellipticity factor Σ% Sum of curvatures &#', (&#' Elliptic integrals of first and second kinds ) Sum of forces applied on a ball * + Jacobian matrix of G , , , Tangential contact forces on inner and outer race for ball j -Deportation distance of normal contacting force due to contact deflection . Ball moment of inertia / , / Primitive radius and diameter 0 Ball mass / Mean rotation angle 1 Spall depth , Δϕ 4 Spall angular position and width around the rotation axis , Δ Spall angular position and width in the race 5 , Radial and axial forces / Driving torque 5 Resistive torque 67 . Ball Pass Frequency on Inner ring 67 8 Ball Pass Frequency on Outer ring 69 Ball Pass Spin Frequency , Fundamental Train Frequency .

Introduction

Bearing dynamic models enable a better understanding of bearing intern behavior, which contributes to the design of diagnostic and predictive maintenance tools. These models are also a good alternative to time consuming and expensive experiments.

The study of the Instantaneous Angular Speed (IAS) originally derives from the measure of Transmission Error in automotive gearbox using optical encoder [START_REF] Remond | Practical performances of high-speed measurement of gear transmission error or torsional vibrations with optical encoders[END_REF]. The IAS has been proven to be relevant in case of non-stationary conditions due to its acquisition method that induces a natural angular sampling instead of a time sampling. This is more convenient for the behavior study of transmission elements such as bearings or gears whose periodicity depends on angle [START_REF] André | On the use of the Instantaneous Angular Speed measurement in non-stationary mechanism monitoring[END_REF], [START_REF] Remond | From Transmission Error Measurements to Angular Sampling in Rotating Machines with Discrete Geometry[END_REF]. Moreover, the IAS conveys information through transmission elements such as gears and enables the monitoring of a whole transmission using a single sensor what makes it very convenient in case of compact systems. Previous works have used the IAS as a monitoring tool. A first model have been elaborated by Bourdon & al [START_REF] Bourdon | Introducing angularly periodic disturbances in dynamic models of rotating systems under non-stationary conditions[END_REF] and focused on the shaft model, introducing the bearing fault as an angular periodic disturbance. A more recent work of Gomez & al. [START_REF] Gomez | Angular modeling of a rotating machine in non-stationary conditions: Application to monitoring bearing defects of wind turbines with instantaneous angular speed[END_REF] has adopted a more phenomenological approach and used a more detailed model of deep groove ball bearing. This work has also studied the benefits of non-stationary conditions.

Bearing dynamic modelling has been well studied throughout the last decades. One of the first work on the subject was carried by Gupta [START_REF] Gupta | Dynamics of rolling-element bearings[END_REF]. Later, Fukata & al. [START_REF] Fukata | On the radial vibration of ball bearings : computer simulation[END_REF] proposed a 2 degrees of freedom (DoF) model to estimate the full force vector transmitted from one ring to the other by using Hertz contact theory. Feng & al. [8] improved the model by adding two DoFs to include frame deflections. Their work also provides a first sliding approach by adding random noise in the angular positions of the rolling elements as well as unbalance effects and localized damages. Impact of fault on the vibration response has been widely modeled. Some authors used forces excitations given by arbitrary laws such as Tandon & al. [START_REF] Tandon | An analytical model for the prediction of the vibration response of rolling element bearing due to a localized defect[END_REF] and Behzad & al. [START_REF] Behzad | A New Model for Estimating Vibrations Generated in the Defective Rolling Element Bearings[END_REF]. Others followed a more phenomenological approach as Sopanen & al [START_REF] Sopanen | Dynamic model of a deep-groove ball bearing including localized and distributed defects[END_REF] and Liu & al. [12]. A more elaborate five DoFs model was later proposed by Sawalhi and Randall [START_REF] Sawalhi | Simulating gear and bearing interactions in the presence of faults[END_REF] introducing radial clearance and damages on rings or balls. This model assumes that rolling elements masses are negligible, that they are in quasi-static equilibrium at any time and that the rotating speed is constant. Gomez [START_REF] Gomez | Modelling deep groove ball bearing localized defects inducing instantaneous angular speed variations[END_REF] adapted their work to non-constant rotating speed and added rolling resistance at every ring-ball contact, introducing by the way a coupling between resistive torque and radial load. Authors showed how a ball bearing creates IAS variation in healthy conditions. They also studied the influence of a spall on the resistive torque applied on the inner ring and the angular speed. This model is however limited to a radial description of the bearing. Other authors such as Hamrock & al. [START_REF] Hamrock | Rolling-Element Bearings[END_REF] designed more complex models of bearings taking into account axial and centrifugal loads without considering rotating speed variations or fault. IAS variations due to rollingelement bearing are of very small magnitude but are nevertheless perceptible with optical or magnetic encoders when sample are long enough. Usual habits associated with classical vibrations study are therefore to be avoided for a pertinent extraction of information provided by IAS.

The model of ball bearing introduced in the first part of this paper is an extension of the work of Gomez [START_REF] Gomez | Modelling deep groove ball bearing localized defects inducing instantaneous angular speed variations[END_REF] to high speed and axial loading. The main improvement is the computation of axial and centrifugal loads on balls in addition to radial loads in order to improve the estimation of normal contact forces, total force vector and characteristic frequencies of the bearing. The contact resolution on spall is also refined by taking account of hybrid contact for which the contact zone covers both a part of the spall and a part of the healthy raceway. It allows the implementation of localized defects that does not cover the whole raceway width. The model is then used in a second part to simulate the behavior of a simple academic system including ball bearings under different operating conditions. This paper is focused on the impact of angular speed and axial load on the bearing behavior. Results from a set of simulations are also presented to demonstrate the limitations and performance of instantaneous angular speed for different loading conditions. Finally, the influence of a spall on the IAS is further studied.

I -Deep-groove ball bearing model description I.1 -Modeling principle and integration into a global system

The modeling principle of a complete rotating system is similar to the one described by Gomez [START_REF] Gomez | Modelling deep groove ball bearing localized defects inducing instantaneous angular speed variations[END_REF]. Shafts are introduced as Timoshenko beam elements. Their global behavior is therefore described by a differential system as bellow:

; + = + = &?' (1)
Where M, C and K are respectively the mass, damping and stiffness matrices of the shaft and they all are constant. C is obtained by a classic modal approach. X stands for the generalized displacement vector. Every node of the shaft is associated to six DoFs corresponding to three translations and three rotations in space. Thus = @ A , B , … , D E where = F , , , , , G is the displacement vector of node i. At this point it is important to remember that no assumption is made on the rotating speed. Hence, is the complete principal rotation of the shaft and does not only account for torsional deformations. Finally, &?' contains the external forces applied at the different nodes. A support or a casing could be introduced by completing the mass and stiffness matrix. A simple approach would be to add some extra nodes as suggested by Gomez. It is also possible to go further and introduce a casing as a super element in order to take account for its eigenmodes properly.

Bearing and other connecting elements are introduced as connecting force vectors. For instance, bearings are assimilated to connecting forces applied on two nodes associated to the inner ring (IR) and the outer ring (OR) whose positions are superposed at rest. The connecting forces depend more on the translations and rotations of the system than on time. This modeling method leads to a strong coupling between the different connecting elements, avoiding any linearization around an operating point and is then more relevant in case of non-stationary conditions. The following general differential system describes the behavior of a simple system such as the one used in part II where only bearings are added to the shaft.

; + = + = &?' + I , = , ; J (2)
Matlab ODE solver enables time integration of the behavior equations. In order to reduce computation time and improve precision, the equations are transposed into the angular domain as Bourdon et al. [START_REF] André | On the use of the Instantaneous Angular Speed measurement in non-stationary mechanism monitoring[END_REF] have suggested.

I.2 -Details of bearing model

The model introduced in this part is developed to estimate the dynamic behavior of a ball bearing integrated in a complex transmission that includes several bearings and gears and is not dedicated to the design of bearings. This is the reason why it is necessary to make some simplifying assumptions in order to limit the model mathematical complexity and the computation time. The following assumptions where retained.

A1:

The cage keeps the balls at a constant distance from their neighbor so that the knowledge of the angular position of the cage is sufficient to determine the position of every ball around the principal rotation axis. A2: At every ball-ring contact a non-slip rolling is imposed. A3: At any moment, each ball is in quasi-static equilibrium in its own radial plan &K LLL⃗, K 5 LLL⃗'.

A4:

The contact forces on the cages are negligible compared to the normal ball-ring contact forces.

A5: Balls and rings are rigid except in a small area around contact points.

The model is a function that takes as an input the displacements of the nodes assimilated to the IR I , , , , , J and the OR I N , N , N , N , N , N J and returns the connecting forces vector. Figure 6 provides the flowchart of the proposed model and each step will be detailed in the following sections. . The calculation of the Fundamental Train Frequency (FTF) is discussed latter.

I.2.1 -Contact deflection

The positions in the radial plan associated to each ball of the centers of curvature of the IR & , ' and OR & , ', introduced in Figure 1 are calculated as described in equations ( 3)- [START_REF] Gupta | Dynamics of rolling-element bearings[END_REF]. These positions are completely determined using the displacements of IR and OR nodes and assumption A5. Let be the radius of the ball, and the radii of IR and OR respectively and and the radii of curvature of the IR race and OR race. The following expressions are obtained geometrically.

= N cosI J + N sinI J + - (3)
= cosI J + sinI J + + (4)

= N + & - 'IsinI J N -cosI J N J (5) = + & + 'IsYZI J -cosI J J (6)
The contact deflection is then calculated depending on the diametral clearance as follow.

1 = - + -4 + \ I - J B + I -J B 1 = - + -4 + \ I -J B + I -J B (7)
Finally, the contact angles are calculated as:

sin& ' = - \ I - J B + I -J B sin& ' = - \ I -J B + I -J B (8)
Figure 1 : Efforts applied on ball j and positions of races curvature centers in the radial plane & , '.

The effect of press fitting of rolling bearings, the centrifugal expansion of the spindle and inner ring, and the thermal expansion of rolling bearing assembly elements are not considered in the proposed model. However, they can be introduced through a modification of clearance used in the computation of the contact deflections. Since is imposed as a parameter, the use of an external modeling tool or experimental data is required to compute it. The addition of clearance variation with working condition calculation could be the subject of future works.

I.2.2 -Normal contact forces and quasi-static equilibrium

Hertz contact theory [START_REF] Hertz | The Contact of Elastic Solids[END_REF] gives contact forces as described in eq ( 9) and [START_REF] Behzad | A New Model for Estimating Vibrations Generated in the Defective Rolling Element Bearings[END_REF].

] LLLLLL⃗ = N . 1 B _ -sinI J -cosI J ` a] LLLLL⃗ = . 1 B _ sinI J cosI J ` (9)
Where:

N = 2√2 3 e. #. ! * f (I# N J I# N J ^Σ% N = 2√2 3 e. #. ! * f (I# J I# J ^Σ% (10)
In these expressions, ! * stand for the reduced young modulus while # and Σ% are the ellipticity factor and the sum of curvatures respectively and depend on the geometry of the bodies in contact. The first and second order elliptic integrals &#' et (&#' calculation is complex and requires an iterative algorithm. A precise approximation method is however given by Antoine [START_REF] Antoine | Approximate Analytical Model for Hertzian Elliptical Contact Problems[END_REF].

In addition to the contact forces, the ball is subjected to a centrifugal force = 0. / . B , where 0 is the mass of the ball and / is the primitive radius. Several authors such as Harris [START_REF] Harris | Rolling Bearing Analysis[END_REF] introduce a gyroscopic momentum, but it is of second importance when computing the ball equilibrium and hardly affects the normal contact forces. Moreover, the tangential contact forces generated by the gyroscopic moment do not generate any torque on the principal rotation axis and therefore do not affect the IAS. The calculation of the gyroscopic moment is complex since it requires the proper rotation angle of the ball. It also generates hyperstatism in the ball equilibrium [START_REF] Bourdon | A new way of writing motion equation in rotating machines by translation into the angular domain[END_REF]. Consequently, the gyroscopic moment has not been involved in the model presented here.

Once all the efforts are known and assuming A3, it is possible to determine the equilibrium position & , ' of the center of each ball j. The Newton-Raphson method is used where the initialization is done with an arbitrary position & ,g , ,g ' for the ball. At each step, the next position is calculated as follow:

_ ,DhA ,DhA `= _ ,D ,D `-* + OA I ,D , ,D J. )I ,D , ,D J (11)
Where G is the sum of forces applied on the ball and * + is its Jacobian matrix whose calculation is performed analytically but is not of prime interest. For clarity purpose, it is not detailed is this paper. Iterations continue until the position of the ball center converges to a zero of G.

I.2.3 -Rolling resistance and tangential contact forces

In order to link load distribution on the balls and the IAS, the rolling resistance has been introduced. Since the model does not aim for a precise calculation of the contact, a macroscopic model describing the rolling resistance is adopted. Under load, the ball and the ring wrap leading to a deportation by a distance b of the application point of normal force from the center of the contact. Figure 2 illustrates this phenomenon. Calculation of b is done through rolling resistance coefficient: -= . 5 . For ball-raceway contact, typical value of 5 is 0.0015. Forces balance on the ball in the axial plane described in Figure 3 gives the following expression of tangential forces [START_REF] Gomez | Modelling deep groove ball bearing localized defects inducing instantaneous angular speed variations[END_REF].

i , = 1 2 I -+ -J + 1 2 & . 0. / . k + .. k ' , = 1 2 I -+ -J + 1 2 & . 0. / . k -.. k ' (12)
Here, k stands for the angular acceleration of the IR while . stands for the inertia of the ball and 0 for its mass. k and k are the cage and ball spin acceleration respectively. As a consequence of assumption A1, k and k are the same for every ball.

The calculation of k and k are difficult because they require the knowledge of the proper axis of rotation of the ball. Complete expressions can be found in Harris' book [START_REF] Harris | Rolling Bearing Analysis[END_REF]. In order to determine the orientation of the proper rotation axis of the balls, Jones [START_REF] Jones | A general theory for elastically constrained ball and radial roller bearings under arbitrary load and speed conditions[END_REF] proposed a "control" method based on the assumption that the motion of the ball is controlled by either one of the rings or the other. A drawback of this method is its binary behavior and the need of a posteriori validation. Changan et al. [START_REF] Changan | Raceway control assumption and the determination of rolling ellement attitude angle[END_REF] developed a more precise and continuous method by applying the dynamic momentum theorem to the ball. Finally, Noel [START_REF] Noel | Complete Analytical Expression of the Stiffness Matrix of Angular Contact Ball Bearings[END_REF] introduced a simpler expression based on the assumption of equal spinto-roll ratio on both ball-ring contacts. In the present model, the expressions implemented for the calculation of characteristic frequencies are simplified and have been established under assumption A2. The angle defining the spin rotation axis of the ball j is computed as suggested by Noel [START_REF] Noel | Complete Analytical Expression of the Stiffness Matrix of Angular Contact Ball Bearings[END_REF]. The balls are assumed to contribute to the cage motion proportionally to their normal contact force magnitude on inner ring. A mean spin angle / is then obtained by taking the weighted average of and is used in the computation of ball and cage rotating acceleration. The spin angle and the mean rotation / change during operation as the contact forces and angle on the rings change.

= + 2 / = ∑ I . m mJ ∑ m m (13) k = / 2 _1 -n / cos / o B `&k N -k ' k = 1 2 pk + k N + / cos / &k N -k 'q (14)
These simplifications have been made because the aim of the model is to give an estimation of characteristic frequencies of the bearing and of their variations under different operating conditions but not to compute them with exactitude. Experience indeed demonstrated that real bearing characteristic frequencies usually differ from the simulated ones because of sliding that occurs at ball-ring contacts. Furthermore, the inertia term in equation [START_REF] Liu | Vibration analysis of ball bearings with a localized defect applying piecewise response function[END_REF] has been noted to be of rather small compared to the normal contact forces term.

I.2.4 -Introduction of damages

In order to compare the behavior of a damaged bearing and a healthy one, it is necessary to implement damages (spall) in the present model. A spall is added as a rectangular pit defined by five parameters &1 , Δ , Δ , , ' described in Figure 4. They correspond to its depth, its angular width, its angular length, its position around the principal rotation axis and its position on the race. For the sake of clarity, this figure is not scaled but a real spall depth is typically close to 100 μm in depth. It is supposed that the presence of a spall only affects the normal contact force magnitude. The contact angle as well as the distance b and the resistance coefficient are supposed to be unchanged. In real cases, the shape of the spall is much more complicated and the pressure profiles are smoothed due to wear and plasticity. This does not affect our study because the IAS variations mainly depend on the integral over time of the loss of contact pressure on the spall. Hence, the use of realistic dimensions for the spall is important. Since the present model is dedicated to spall detection and not lifetime estimation, the propagation of the spall is not considered here.

Figure 4 : Spall definition parameters

Because of the ball and ring deflections, the shape of the contact area on healthy ring is an ellipse. Depending on the mating force and on the geometry of contacting bodies, the ellipse may cover partially or totally a spall. Consequently, the Hertz contact model can't be directly used. A more advanced resolution of contact has been proposed and implemented. For each ball, the shape of the contact area and the pressure distribution on healthy races are calculated according to Antoine & al. [START_REF] Antoine | Approximate Analytical Model for Hertzian Elliptical Contact Problems[END_REF]. Here, r and s are the normal contact forces calculated for a contact on a healthy race and a contact in a wide spall respectively.

is the total normal force and C is the contact coefficient introduced in I.2.2.

r = . 1 B tZu s = . &1 -1 ' A ^ (15) v w & , ' = 3 2 w e. t w . -w f 1 -_ t w `B -_ -w `B k ∈ @ℎ, zE (16)
where If the contact ellipse and the spall intersect, the contact pressure in the corresponding area is replaced by the contact pressure that would be obtained for a deflection reduced by the depth of the spall. Figure 5 shows the pressure distribution obtained with this method for a shallow spall. The total normal contact force is then computed by integrating the mixed contact pressure obtained. The contact resolution method presented here is not as precise as a complete resolution using advanced methods like finite elements or Boussinesq theory [START_REF] Boussinesq | Application des potentiels à l'étude de l'équilibre et du mouvement des solides élastiques[END_REF]. It does not account for the pressure concentration on the spall edges. It allows however computation times that are compatible with dynamic simulations and provides a correct estimation of the contact force. It should be noted that a real spall is not with a rectangular shape. Thus, a perfect resolution of the contact would not be relevant anyway.

t

I.3 -Validation in quasi-static conditions

In order to validate the simulation of the ball equilibrium position and the global connecting force vector, the present model has been compared to RBS2, a module of RBSDYN, which is a numerical tool developed by the CNES and available online. The geometry of the bearing used for this validation corresponds to a real bearing from the power transmission of an helicopter engine produced by Safran. The nominal working conditions of this bearing are known. Computations have been carried for different loading conditions. For each computation, the OR was built in while a displacement composed of translations and rotations was imposed on the IR. The magnitudes of the imposed displacements were chosen so that the resulting efforts were close to the nominal operating conditions of the bearing. The validation focused on the normal contact forces and contact angles which are the parameters involved in the computation of global connecting efforts, rolling resistance and resistive torque. Only one azimuthal location for the ball was studied during this validation. An exhaustive description of the validation is out of the scope of this paper. However, Figure 7 present the results obtained by simultaneously imposing axial and radial displacements, and an angular misalignment on the IR. In this figure, the normal contact forces and the contact angles on IR and OR are plotted as functions of the ball index which vary from one to eleven (the number of balls in the bearing in this example). Weak differences have been observed for the least loaded balls for which the centrifugal forces computed by RBS2 were significantly different from the ones obtained by the present model. This phenomenon is a direct consequence of the assumption A1 and of the simplifications in the calculation of the rotating speed . Yet, the differences are negligible and only impact the least loaded balls, which contribute the less to the dynamics of the bearing. The present model has therefore been validated in quasi-static conditions.

It must be noticed that in both RBS2 and the model introduced here, the contact ellipse is assumed to be included in the groove. Hence, the contact angle is not mathematically restricted and Hertz's contact law is assumed valid. This can lead to unphysical contact angle on inner race as for ball 6. The value of this angle is due to the very light normal IR load. The centrifugal force is then of primary importance and keep the ball deep in the OR which induce high contact angle on IR. Since this phenomenon only appears for light loaded balls, a more complex contact model is not necessary required.

II -Model exploitation

The bearing model presented in the first part of this paper has been integrated into a simple shaft model is order to highlight changes on the bearing behavior under different operating conditions.

II.1 -Simulated system description and reference run

In order to limit the complexity of the generated signals, the simulated system, presented Figure 8, has been designed as simple as possible. It is composed of a constant section shaft modeled by four Timoshenko beam elements of equal length. The nodes are displayed in red while the bearings are sketched by two balls only. A rigid disc has been added at the middle of the shaft to ensure a realistic inertia. A constant driving torque / and a viscous resistive torque are set on two nodes invisible in Figure 8 and related to the shaft extremities through elastic couplings. The torsional stiffness of elastic coupling is set to 1300 Nm.rad -1 . The resistive torque is defined as 5 = -}. where } is computed so that the resistive torque equals the driving torque for a given nominal rotation speed. Finally, an axial force and a radial force 5 are imposed on the central node. The bearing geometry involved in the simulations comes from a real bearing and parameters are given in Table 1. The number of balls has been reduced on the bearing located on the right side of the shaft to dissociate the bearing signatures. Hence the BPFOs of the bearings are 4.63 and 4.26 events per revolution respectively. In order to be as close as possible to a classic assembly, all the DoF of the right bearing OR are blocked while the axial translation of the left bearing OR is free. Hence, both bearings are radially loaded but only one of them supports an axial load. Furthermore, no mounting error is added to the system. However, due to the shaft flexibility, both bearing IRs are slightly tilted. For simplicity and confidentiality purposes, all the components of the models are associated to the same material. Several simulations were performed for different operating conditions summarized in Table 2 : Operating condition for each simulation. A first simulation was performed to stand as a reference. The operating conditions were calculated to be close to real nominal operating conditions of the bearing. A nominal rotating speed of 6 000 rpm is imposed as well as a radial force of 600 N and an axial force of 200 N. Every result will be plotted against the angular position of the shaft central node instead of time. In the same way, the spectrum will be presented versus angular frequencies instead of temporal frequencies. The associated unit is "event per revolution" (ev/rev) by analogy to the Hertz unit. The model presented here enable the study of both classic acceleration signals and IAS. Figure 9 shows complete spectra of shaft acceleration and IAS obtained for this reference run as examples. Acceleration of OR (or casing) could be easily obtained by integrating them into the modeled system but they are not of prime interest here. As mentioned earlier, in the rest of this paper, we will focus on the IAS and more precisely on the lower band of its spectrum and on its shape in the physical domain. Figure 10 shows the IAS of the shaft for this reference simulation. In this figure, the peaks of the BPFOs are spotted with black squares and their harmonics are highlighted with black dots to distinguish them from the modulation and eigenmodes peaks.

Figure 9 : (a) Acceleration and (b) IAS complete Spectrum

The characteristic frequencies of both bearings are marked in Figure 10b. In this plot, the BPFOs of both bearings are pointed out with squares and their harmonics with dots. Concerning the axially loaded bearing, only the fundamental and the first harmonic are visible in the spectrum. The excitation due to the radially only loaded bearing is significantly stronger and composed of more harmonics. An oscillation of low amplitude and high frequency is noticeable in Figure 10a which is related to the system response. Finally, modulation side peaks are also visible around the peaks associated to the BPFO of the radially loaded bearing and its harmonics. The distance between the BPFO peaks and the modulations ones is &67 8 5 -67 8 '. It is important to notice that the modulations observed here correspond to the actual modulation phenomenon and are an output of the model. They are not imposed as an input and are not completely predictable a priori. Thus, this model supports strong couplings between bearing and shaft dynamics. In a more complex system including casing and gear, this would be a relevant tool to determine the location of frequencies of interest for bearing monitoring which are not always the BPFO and its harmonics. Again, all these IAS variations are of very small magnitude but they are deemed to be measurable by means of encoders. In order to highlight the system behavior changes under different operating conditions, all following results will be compared to this reference run.

II.2 -Nominal rotation speed influence

A change of nominal rotation speed has two major impacts on the bearing behavior. Firstly, it affects the characteristic frequencies of the bearings which are proportional to the rotation speed. Secondly, it modifies the centrifugal forces imposed on the rolling elements and thus modifies their quasi-static equilibrium. In order to observe the influence of the nominal rotation speed, a simulation has been performed at 20 000 rpm just by changing the resistive torque applied on the shaft. Every other parameter is kept equal to the reference run. Figure 11 shows the IAS variations obtained in this simulation compared with reference results. The mean rotation speed is detrend in order to compare IAS variations. Figure 11b compares the results obtained for nominal run and with higher nominal speed. The fundamental frequencies of the bearings are spotted with "+" and "x" markers respectively. The different harmonics are associated to square markers.

A reduction of the fundamental frequencies' magnitude and an increase of harmonics magnitude can be observed at higher rotational speed. It can be explained by the increase of centrifugal forces which maintain the balls deep in the outer race and increase the contact deflection on the OR. This creates an additional equivalent clearance and decreases the size of the angular interval of contact on the IR. Figure 12 shows the evolution of the contact normal force on IR. The lower curves represent the normal load on ball number 1. The angular interval of contact on IR decreases from 0.6 to 0.5 revolution of the IR for the radially loaded bearing. This phenomenon is even more visible for the axially loaded bearing whose balls were in contact with the IR during the whole rotation in the reference simulation, which is no longer true for higher speeds. The upper curves represent the evolution of the sum on every ball of normal contact load on IR for both bearings. This sum is closely related to the resistive torque in the bearing. It becomes more complex and more distant from the shape of a sine at high speed. This results in a decrease of fundamental magnitude and an increase of relative significance of harmonics in the IAS.

The natural response of the system is also more visible at high speed and peaks of eigenmodes appear. The time-frequency of eigenmodes does not depend on the rotational speed. This mean that their angular frequency is proportional to the inverse of the rotational speed. The more visible are peaks associated to the elastic coupling torsion modes located at 0.39 ev/rev and the second order torsion mode of the shaft at 25.44 ev/tr (red dots). This can be explained by the decrease of their angular frequency which is closer to the low harmonics or bearings BPFOs at high speed. The more impulsive shape of the resistive torque is also more favorable to the excitation of every frequencies. The modulations are more present in this simulation than in the reference one due to the increased number of available harmonics of the bearing BPFO. Finally, it is noticeable that the BPFO of the exclusively radially loaded bearing does not change while the BPFO of the axially loaded bearing decrease by 0.4% due to the modification of contact angles.

II.3 -Axial load influence

A variation of axial load has consequences on balls equilibrium, the most notable of which are contact angles modifications, contact stiffness changes (due to the non-linear contact law) and modifications of load distribution on rolling elements. Since only one bearing is axially loaded the behavior of the second one does not change and will not be discussed. A simulation with an axial load of 1 000 N has been performed and Figure 13 compares the IAS variations obtained with the reference results. Again, the BPFO and harmonics peaks are marked with red squares for the 1 000 N spectrum.

At high axial load, the excitation due to the bearing (BPFO and first harmonic) are not visible in the spectrum. Furthermore, the different modulations visible at 200 N also vanish at high load. This is explained by the fact that every ball is loaded through the whole rotation and the resistive torque due to rolling resistance phenomenon hardly varies. Since the BPFO is an output of the model, it is still possible to track its variations. In this case, the frequency has increased by 0.4% compared to the reference run. Other simulations, not presented here, were performed for intermediate axial loads and showed that, when axial load increases, the magnitude of the excitation due to the axially loaded bearing gradually decreases while the BPFO increases.

II.4 -Spall influence

A spall alters the dynamic response of bearings. In addition to the vibrations widely studied in the literature [START_REF] Sawalhi | Vibration response of spalled rolling element bearings: Observations, simulations and signal processing techniques to track the spall size[END_REF]- [START_REF] Antoni | Contrôle et diagnostic à partir des signaux acoustiques et vibratoires[END_REF], the passage of a ball above a spall induces a modification of the load distribution on the rolling elements and the resistive torque. A simulation has therefore been performed introducing a spall on the radially-only loaded bearing OR. It has been located in the maximum load angular area which means that it is in the direction of the radial load. This area has been chosen as it is the fastest one to degrade in actual conditions. The depth of the spall is set to 100 μm while its angular width is 30 degrees and corresponds to 25% of the race width. The angular length of the spall is set to 7 degrees which corresponds to about 20% of the angular space separating two adjacent balls. Figure 14 compares the obtained results to the reference ones. It appears that the presence of damage increases the IAS amplitude. It is possible to detect the entrance and the exit of a ball in the damaged area in the IAS signal in angular domain. Two vertical doted lines have been drawn in Figure 14(a) to highlight these events. The excitation due to the spalled bearing is close to an impulse and strongly excites the shaft eigen modes. Thus, the oscillations after the entry or the exit of a ball in the damaged area are generated by the response of the system. The excitation of the axially loaded bearing is now negligible compared to the one of the damaged bearing and modulations vanish. Once again, the nominal rotation speed has an impact on the IAS variation. A second simulation was performed introducing a spall at 20 000 rpm. As noticed before, the eigen frequencies of the system are constant in time domain and thus vary as the inverse of the nominal speed in angular domain. Even if the excitation due to the spall is similar at 6 000 and 20 000 rpm, the system response changes. Figure 15 compares the IAS variations in angular and frequency domains for both nominal speeds. As an example, the main torsion mode is located around 85 ev/rev at 6 000 rpm and 26 ev/rev at 20 000 rpm. For each run, the harmonics of the BPFO close to this mode clearly stand out. The excitation related to the spall also slightly differs because of the change of centrifugal forces. In Figure 16, the evolution of normal contact force on the IR for both runs is displayed. It is visible that the loss of charge at the entry of the spall is total at 20 000 rpm due to the centrifugal effects and the higher contact deflection on OR, which is not the case at 6 000 rpm. During the passage of a ball above the spall, the inner ring slowly moves to a new equilibrium position and the loads on every ball increase. This displacement is smaller at high speed because the time needed for a ball to pass above the spall is smaller. It is also important to note that the load on the ball right before the entry of a ball in the spall is smaller at low speed and the loss of normal force is therefore less significant. This phenomenon is purely related to dynamic effects and can widely differ for different systems or rotation speeds. It cannot be considered as a general observation.

Other simulations have showed that introducing the spall on the axially loaded bearing OR leads to smaller IAS variation amplitude increases because of the axial load. Due to the axial load, the load distribution is more homogenous and a decrease of contact load on one ball has a smaller impact on the resistive torque and consequently on the IAS. Moreover, depending on the contact angle, the portion of the contact ellipse that passes over the spall varies and affects the load decrease on the ball. However, the spall signature never disappears, even when the load is purely axial as long as the contact ellipse intersects with the spalled area. As an example, results obtained in the presence of a spall on the outer ring of the axially loaded bearing is given in Figure 17. For clarity purpose, a zoom on low frequency is performed in Figure 17b. The BPFO and harmonics of axially loaded bearing is now of same amplitude as the BPFO of the radially loaded one leading to numerous modulations that complexify the spectrum.

III -Conclusion

The deep-groove ball bearing model introduced in this paper is based on the Gomez' work and has been improved in order to take axial and centrifugal loads into account, to make it more suitable for high speed systems. The quasi-static equilibrium is then established accounting for nonzero contact angles on both races. A quasi-static validation was performed using commercial tools provided by Safran Helicopter Engines. Further validation based on dynamic experiment need to be carried out in order to ensure the validity of the assumption made in this model. An improvement of the model could also be carried out by enhancing the calculation of tangential contact forces. Finally, an adaptation of this work to three or four contact point bearings could be achieved in the future.

The model has then been investigated to simulate bearing behaviors under different operating conditions. For this purpose, an academic model of a rotating shaft has been elaborated in order to mix excitation from bearings. Several simulations where performed with different axial load and nominal speed. Results showed that bearings operating at high speed behave differently due to stronger centrifugal forces on balls. Furthermore, it appears that an axial load increase on a bearing leads to the homogenization of the load distribution on balls, to an increase of contact angles and to a mitigation or even a vanishing of the excitations related to this bearing. Eventually, the correct estimation of contact angle is important to introduce for a spall that does not cover the whole width of the race. Simulations involving spalls on bearing races showed that the presence of damages leads to an increase of the amplitude of IAS variations and a stronger excitation of the eigenmodes. However, this increase depends on the axial load which must therefore be determined precisely.

The minimum size of a spall that would be detectable using such a modeling is difficult to determine. Mathematically, in the simulation, any fault would affect the IAS no matter its size as long as it is located in the zone of load and its area is not zero. Some indicator based on the amplitude of the IAS variation due to bearing could be used, extracted from past experimental studies. Furthermore, as seen in this paper, the influence of a given spall on IAS highly depends on the axial and radial forces as well as the nominal rotating speed. A complete study would be required to determine the minimum detectable size of spall for different working conditions, different shape and location of spall. This remains for further works.

The simulations presented here were performed for stationary conditions. In case of large variation of rotation speed or load, the results can be even more complex, even for a system as simple as the one introduced here. Furthermore, real transmissions are composed of several bearing, and gears and the casing usually also influences the global dynamics. The influence of each part of the transmission on the vibrations and IAS are not predictable a priori and strong couplings are necessary in the model to estimate these governing couplings. To that goal, future works will focus on gear model in order to propose a complete modeling of a realistic transmission.

  Each ball j is located around the principal axis of rotation through an angle . The assumption A1 leads to a simple calculation of = + & OA'BP Q where is the angular position of the cage around the principal rotation axis which is obtained by integrating the angular velocity of the cage. The rotation speed of the cage is linked to that of the IR by the relation: = , .

Figure 2 :

 2 Figure 2 : Rolling resistance model Figure 3 : Force balance on the ball j in axial plane

Figure 5 :

 5 Figure 5 : Example of pressure distribution for a contact on a spall

Figure 6 :

 6 Figure 6 : Flowchart of the bearing model

Figure 7 :

 7 Figure 7 : Comparison of results provided by RBS2 (red) and by the proposed model NT (blue): (a) normal contact forces and (b) contact angles

Figure

  Figure 8 : Simulated mechanical architecture

Figure 10 :

 10 Figure 10 : IAS in (a) angular and (b) frequency domains (6000 RPM, Fa = 200N, Fr = 600N, no spall)

Figure 11 :Figure 12 :

 1112 Figure 11 : IAS variations in (a) the angular domain and (b) the frequencies domain. Influence of nominal rotation speed (Fa = 200N, Fr = 600N, no spall)

Figure 13 :

 13 Figure 13 : IAS in (a) angular domain and (b) frequency domain. Axial load influence (6000 RPM, Fr = 600N, no spall)

Figure 14 :

 14 Figure 14 : IAS variations in (a) angular domain and (b) frequency domain. Influence of a spall (6000 RPM, Fa = 200N, Fr = 600N)

Figure 15 :Figure 16 :

 1516 Figure 15 : IAS variations in presence of fault in (a) angular domain and (b) frequency domain, influence of the nominal rotation speed. (Fa = 200N, Fr = 600N, spalled OR)

Table 1 :

 1 8 : Simulated mechanical architecture Geometry and materials properties

	Parameter			value
	Shaft intern diameter	30 mm
	Shaft outer diameter	40 mm
	Shaft total length		400 mm
	IR diameter		67.3 mm
	OR diameter		92.7 mm
	Number of balls		left bearing: 11	right bearing: 10
	Specific mass		7800 kg/m3
	Poisson's coefficient		0.3
	Young's Modulus		200 GPa
		Reference run High speed run	High axial load run Damaged OR run
	Axial load	200 N	200 N	1 000 N	200 N
	Radial load	600 N	600 N	600 N	600 N
	nominal speed	6 000 RPM	20 000 RPM	6 000 RPM	6 000 / 20 000 RPM
	presence of spall No	No	No	Yes

Table 2 :

 2 Operating condition for each simulation

Acknowledgements

This work was supported by the French National Agency of Research and Technology (ANRT), CIFRE [grant number 2019/1353] and Company SAFRAN Helicopter Engines. The work is also involved in the framework of the CIRTrans. All of them are gratefully acknowledged for their support.