Three-Dimensional DEM Modelling of Ball Bearing with Lubrication Regime Prediction

This paper deals with an efficient 3D modelling of a radial ball bearing to predict the operating lubrication regime under mechanical loading and mounting conditions by using the Discrete Element Method (DEM). Due to the relevance of such an approach, especially for multicontact systems, the lubrication regime associated with specific operating conditions can be predicted accurately. By means of an elastohydrodynamic lubrication formulation depending on parameters related to the size of contact area, mechanical properties of materials, roughness and fluid viscosity, the lubricant film thickness is predicted and used to take into consideration the fluid film damping effect and friction coefficient variation. The lubrication regime can be identified according to Stribeck curve with the assumption of a piezo-viscous-elastic behaviour of the lubricant. The numerical simulations performed with MULTICOR-3D software on an operating ball bearing shown that the lubrication regime at the rolling element-raceway contact can be easily monitored and quantitatively identified. To assess the efficiency of the discrete modelling, a parametric study is carried out in order to exhibit how the operating conditions affect the lubrication regimes and the fluid film spread in the loaded zone. The adequacy between the choice of lubricant and the bearing tribofinition is sought to optimize the component lifetime.


Introduction
Rolling bearings are widely used components in both industrial and domestic rotating machinery.There is a large number of bearings for all possible applications, among them radial and oblique ball bearings, roller bearings, ball and roller thrust bearings, and many others.However, statistical studies show that in industrial applications more than 50% of malfunctions of rotating machineries are caused by shaft misalignment [1].Therefore this recurrent defect unfortunately plays a central role in the degradation of bearings and malfunctions in rotating machinery.In the modern area of the smart factory (Factory 4.0) [2,3], when attempts are being made to obtain relevant real time information on rotating machines, using connected devices called Internet of Things (IoT) is sought, and also with increasing demands on industry to reduce energy consumption and gas emissions, the effort to increase the efficiency of machine components is no doubt more intense than ever.In order to ensure the industrial systems' availability and the safety of goods and persons, the monitoring and diagnosis of bearing malfunctions have to be considered with prime importance, thereby ensuring good maintenance.Moreover, the challenges in terms of productivity and economic cost are obviously non-negligible.Therefore, defect tracking of rolling bearings components has led to extensive research over many years.Several experimental methods and various types of signal analysis have been developed for the detection and diagnosis of rolling bearing defects, such as vibration measurements [4], acoustic emission [5], defect signatures in the stator current of motors [6], thermal and lubrication analysis [7].
However, another very important and key aspect on which the performance of rolling bearings depends concerns lubrication.Because of its preponderant role, lubrication in fact plays a vital role in most machine components in ensuring that they work properly, are efficient and have a satisfactory service life.With clearance and alignment correction, lubrication is the main maintenance operation, without dismounting.Hamrock and Dowson [8,9] performed pioneering studies on lubricated contacts.The authors implemented a series of EHL (Elasto-Hydrodynamic Lubrication) numerical analysis for elliptical contacts and considered a wide range of load, velocity and material properties variations.Several curve-fit formulas derived form their numerical results characterize different fluid film thickness regimes, and which are still the most popular and frequently used today.More recently, Masjedi and Khonsari [10] have improved the EHL model, initially introduced by Hamrock and Dowson, proposing formulas to take into account the effect of surface roughness to predict central fluid film thickness, minimum film thickness and the asperity load ratio.It is, therefore, strongly advised that mechanical systems with non-conformal contacts should be mounted with sufficient clearance because the lubricant film should be thick enough as far as possible to avoid contact between the mating surfaces.This ensures the limitation of wear mechanisms that could occur at the contact interface between two surfaces.In this case, the friction coefficient is rather low.On the other hand, when the rotating machines are subjected to more severe operating conditions, typically under very high loads or at low speeds, then the load carried by the contacting asperities should be taken into consideration.According to Stribeck analysis [11], the relation between the friction coefficient µ and the fluid parameter Λ f (Sommerfeld number or Hersey number) is often used as an indicator to distinguish between the different lubrication regimes.Moreover, the increasing developments related to tribological studies have considerably enhanced our understanding of the mechanisms of lubrication that take place in highly loaded non-conformal contacts, in the case of rolling bearings, gears and cams [12,13], as well as for the identification of the fluid film thickness [14,15] and the characterization of the friction behaviour [16,17].
Wear and fatigue damage behaviour in bearings is hard to describe theoretically and numerically.Indeed, both wear and fatigue mechanisms are quite complex and essentially dependent on operating conditions, lubricant quality associated with a good mounting and a suitable bearing design.To overcome these limitations, a numerical twin of bearings, based on DEM using MULTICOR-2D software, has been developed to enable users to access outputs relative to damage initiation and to improve coupling models for monitoring [18,19].Moreover, since 2008, several studies carried out in our research laboratory (LTI EA-3899) have clearly established the relevance of electrical measurements located on dynamic interfaces and have generated a promising approach for the monitoring of bearings and preventing material wear in wind plant applications [20,21].By using a refined electrical resistivity circuit, some recent work has also shown the relevance of electrical measurements in characterizing the lubrication regime at the contact interface [22] and investigating the influence of the spin coating process on the oil film lifetime [23].Combining an improved electrical signal, where the lubricated contact and friction play a key role, with the existing monitoring methods [24] should lead to a high-performance rotating machine tracking tool.From an industrial point of view, obtaining real-time information including a qualitative definition of the rolling contact quality based on electrical measurements would allow the lifetime of rotating machines to be optimized.From a numerical point of view, generating mechanical responses of a bearing using DEM allows us to go further in the knowledge of monitoring signals, instead of investigating mechanisms responsible for wear.For example, wear debris can contaminate lubricated bearings and impede contacts resulting in lubrication failure and subsequent overheating.Our recent article [18] exhibits the potential of DEM simulations to investigate fatigue mechanism by accessing stress in rings.
A focus on the potential offered by the DEM for simulating the dynamic behavior of bearings with a 2D description is well described in our previous work [25,26], assuming a simple contact model without tribology enrichments.In this paper, a 3D DEM model is enriched with the phenomenological contribution of lubrication, in order to develop a monitoring tool based on an electrical measurement [21].In this sense, the lubricant consideration must remain sufficiently adapted to the scale of the apparent contact, where the competition between the metal-metal contact with asperities (resistive) and the contact through the fluid (capacitive) plays an essential role in electromechanical coupling.To achieve this proposal, from a mechanical standpoint, a three dimensional dynamic model for ball bearings based on MULTICOR-3D software is developed.The aim is to create a lubrication model by using a discrete formulation.The proposed model implements a major part of the functionalities according to the state of the art, such as the slice model, the mutual influences of neighbouring elements with clearance, elastic contact models, lubricant film estimation, damping forces, centrifugal forces, etc.In particular, the DEM model with strong assumptions introduces the basics of lubrication regimes in order to maintain efficiency in solving, as long as it could be simulated in an acceptable computational time.Technically, the smooth DEM methodology selects stiffnesses and damping coefficients (critical damping) and a constant coefficient of friction for simulations related to granular media [27,28].Following this tutorial dedicated to bearing application, stiffnesses are calibrated for bearings and typical load-deflections relationship are validated [20,29].Dampings [30] and friction [31] acquire a more physical meaning when a lubricant is considered [8,32,33].The friction coefficients, usually considered as an input parameter are now implicitly set at each contact.When the roughness is known, the parameter that indicates a good lubrication regime is the minimum fluid film thickness.This parameter must has values that do not fall below the roughness because a film thickness that is too small leads to risks of metal-metal contact, causing irreversible damage.Future numerical developments and experimental investigations will attempt to define electromechanical coupling models implemented in DEM simulations, based on the capacitive and resistive behaviour of lubricating contacts.When lubrication is introduced under a restrictive assumption on the fluid rheology, the resulting damping forces should be included.
The paper is divided into four sections, including two sections dedicated to the introduction and conclusion, respectively.The second section deals with a three-dimensional discrete modelling performed on a ball bearing, with a powerful computational tool to manage several contacts under dynamic loadings.The simulated radial ball bearing of 6208 series is designed fully and accurately in accordance with its geometrical conformity.Moreover, a mechanical formulation, based on contact stiffness for elliptic areas and viscous damping effects is developed.With the assumption of piezo-viscous-elastic fluid behaviour, an EHL model is implemented in the numerical software allowing the determination of the viscous damping due to the fluid film.Numerical simulations are carried out with MULTICOR-3D software to validate the proposed discrete model of radial ball bearing and are presented and discussed in Appendix A. The third section treats the qualitative prediction of the lubrication regimes in an operating radial ball bearing.To exhibit the relevance of the DEM approach, several operating conditions are considered.The effects of radial load, diametral clearance and shaft angular speed are considered and discussed, from a lubrication regime point of view.In the conclusion, we summarize the contributions of this work and we end with a few thoughts on future work.

Three Dimensional Modelling of Radial Ball Bearing Using DEM
To overcome the 2D modelling limitations from the out-of-plane effects (gyroscopic effect, axial loading, free contact angle, etc.) are ignored [20], we proposed a three-dimensional description of a ball bearing [29] using DEM (Figure 1b) with the aim of investigating the lubrication regime in relation to the kinematic and mechanical operating conditions.The radial ball bearing of 6208 series (SKF 6208 TN9/DBGA) is considered in this study (Figure 1a).
The dimensions and geometrical conformity of the bearing considered are listed in Table 1.The number of balls in contact with the inner and outer raceways is Z = 9.Nine additional smaller balls of radius R z in contact with two toruses are considered to play the role of a cage (Figure 1c).Even if the cage is different in shape in comparison to the standard form, the design does not at all affect the contact forces between rolling elements and cage [25].In addition, the proposed description of the cage mainly has two roles.The first is that it ensures a conform internal radial load distribution, while the second is that it allows intermittent rolling elements/cage contacts to be taken into account when the ball bearing operates.However, details should be supplied about the cage in the sense that the only other possible interactions involving the cage elements may be occur between rolling elements.The ball bearing design performed by DEM allowed us to numerically achieve bearing frequencies, namely FTF, BPFI, BPFO and BSF with a good accuracy and gives the opportunity to take into account all bearing components without any restriction.

Contact Stiffness and Damping Effect
The contact stiffness model used in the discrete modelling is based on a smoothed formulation established by Cundall et al. [27].The contact forces are obtained by means of Kelvin-Voigt Spring-Dashpot Model and a Coulomb's law if sliding occurs at the contact (Figure 2).The main assumption made with DEM is that the particles in contact are assumed perfectly rigid and therefore, there is no elastic deformation of particles.The elasticity is only considered at the contact between particles.The contact forces F n,t in normal and tangential directions according to the local frame at the contact plane, are described with an explicit mechanical model depending on elastic force displacement law, Coulomb friction and viscous damping coefficients (Equation ( 1)).
where δ n,t are the normal and tangential relative displacements, K n,t the normal and tangential stiffnesses, C n,t the normal and tangential viscous damping coefficients, v n,t the normal and tangential relative speeds according to the local frame and µ the Coulomb friction coefficient in the case of a sliding contact.The normal viscous damping coefficient C n to ensure a satisfactory mechanical steady state is a function of two damping types, namely hysteretic and elastohydrodynamic, related to the contact shape, roughness, fluid properties, temperature, operating conditions, etc.In considering a reasonable assumption, the expressions of both hysteretic and fluid damping coefficients are discussed later.The tangential viscous damping coefficient C t is more difficult to evaluate.For convenience, the ratio C t C n is considered to range in the interval [0, 1].The normal contact stiffness K n is related to mechanical characteristics and curvatures of the surfaces in contact according to non-linear model ; n .Considering a rolling element in contact with inner-raceway (or outer-raceway) the normal stiffness is given as function of the curvature sum radius and the approximate elliptic integrals as proposed in [34]: with R curve the curvature sum radius, R x the effective radius in x direction, R y the effective radius in F the approximate elliptic integral of first kind F = 1 + q lnα, q = π 2 − 1 and E the approximate elliptic integral of second kind E = 1 + q α (Figure 3).The equivalent radius R x is calculated by taking into account the cases of convex contact (ball-inner-raceway) and concave contact (ball-outer-raceway), as given in Equation ( 2).For the sake of clarity, we note that is the shear modulus of steel in the case of identical materials in contact, E the Young's modulus (E = 210 GPa) and ν the Poisson's ratio (ν = 0.3).When interactions between rolling element and cage ball occur the contact is assumed to be in boundary lubrication regime without friction dissipation and the damping coefficient is considered critical.The normal stiffness is only taken into account, and therefore the following equation is considered [35]: where R eq = r i r j r i +r j is the equivalent radius of two balls in contact of radius r i and r j , respectively.
The tangential stiffness is related to the normal stiffness so that the ratio K t K n ranges over the interval . The tangential stiffness K t can be calculated by means of the normal interaction F n , the curvature radii and the mechanical properties of both rolling element and cage ball [35].
The steady regime of the ball bearing in dynamic mode is achieved with viscous damping and local friction added as dissipative forces to the mechanical model (1).The primal viscous effect is the hysteretic damping coefficient C hyst in the normal direction n and related to the plastic strain of surface roughness due to the normal interactions (4).
where α e = 0.08 s•m −1 denotes the restitution coefficient of steel [36].In the case of a boundary lubrication regime in ball bearing, typically at heavy radial load and low rotation speed, this means that metal-metal contact predominates, involving a high coefficient of friction (µ 0.1).Bearing, cam mechanisms and gears are considered with non-conformal contacts endowed with EHL regime.Such mechanical systems are often designed for the lubricant film to be thick enough, more than 1 µm, to avoid contact between the mating surfaces.According to the relation given in [30], the damping coefficient C f luid due to the lubrication film in the normal direction of contact is written as follows (5): with η = 0.04 Pa•s the dynamic viscosity of the fluid at operating temperature and atmospheric pressure [37], a the semi-major axis of the Hertz's elliptical contact [9], h min the minimum fluid thickness of the ball-inner-raceway (or -outer-raceway) contact at an azimuth angle ψ.

Fluid Film Thickness and Fluid Parameters
In the field of elastohydrodynamic lubrication, grease rheology influences the fluid film thickness under operating conditions [38].Given that the viscosity of lubricant increases with respect to pressure, this means that the lubricant exhibits piezo-viscous behaviour.Temperature also plays a role in grease viscosity changes.Taking into account the piezo effects of the fluid, we have, in the discrete modelling, considered the piezo-viscous-elastic regime to evaluate the minimum fluid film thickness h min [34].For hard surfaces in contact, h min is expressed according to the following elastohydrodynamic Formula (6): E, U r , W are three dimensionless parameters: E = ξE , with ξ the viscosity-pressure coefficient of the lubricant at operating temperature and atmospheric pressure (ξ the effective elastic modulus of steel ; U r = ηU r E R x , with U r the rolling velocity ; W = , with Q ψ the ball-inner-raceway (or -outer-raceway) normal load (Figure 4).It should be noted that the deformation of the surface asperities is not taken into account in the fluid film thickness model.Masjedi and Khonsari [10] proposed a curve-fit expression for the film thickness and asperity load ratio in point-contact EHL of rough surfaces, in order to take into account the deformation of surface asperities.In the case of a deep groove ball bearing, the authors obtained a very low asperity load ratio for medium angular shaft speed (ω sha f t > 200 rad•s −1 ).However, for a rotational speed of one order of magnitude lower, the asperity load ratio at both outer and inner races reaches a value of the order of 10% to 12%.In the context of this work, the numerical simulations are at least performed with angular speed of two orders of magnitude.Therefore, it is reasonable to assume that the roughness deformation slightly affects the fluid film thickness.
Commonly, a simply way to identify which lubrication regime a system is running in is by the use of a dimensionless parameter, Λ f , proposed by Tallian and called the fluid parameter [39].Λ f is the ratio of the minimum fluid film thickness to the surface roughness.This parameter is also related to the friction coefficient µ as shown through a schematic representation of the Stribeck curve [11] for non-conformal contact (Figure 5).
where ω is the shaft angular speed (i.e., inner ring), σ a1 and σ a2 denotes the arithmetic mean roughnesses of the measured surfaces for the inner-raceway (or outer-raceway) and rolling elements, respectively.The finishing process known as "Superfinishing" enables surface condition values of a dozen nanometres [40,41].We assume that the races have σ a1 ∼ 3 × 10 −7 m and rolling element roughness can be neglected (σ a2 << σ a1 ).Ball-on-disk tests, known as tribometer usually involve similar orders of magnitude [15,42].At fully flooded condition, the normal damping coefficient C n is strongly dependent on the EHL regimes.Since the rheology of the fluid is taken into account through the assumption of a piezo-viscous-elastic model, we postulate that regimes simulated under the considered operating conditions do not increase K n with an additional fluid film stiffness by considering K n (Λ f ) [30].As long as the mechanics of the surfaces are responsible for the stiffness of the component, within the piezo-viscous-elastic limits, the film parameter allows us to select the "predominant" viscous damping C n (Λ f ) according to the following Table 2.
When the applied load is mainly carried by asperities of the surfaces in contact, this means that only the hysteretic effect plays a role in normal damping, C n = C hyst .In mixed lubrication regime the normal damping C n is considered to be a function of both hysteretic and fluid damping.By assuming that the two effects act in parallel, then C n = C hyst C f luid C hyst +C f luid .In the case of fully developed lubrication film, the hysteretic effect is not possible and only the fluid effect remains, so C n = C f luid .Figure 6 gives an illustration of viscous damping status according to the fluid parameter Λ f .This approach might appear simplistic in view of the complexity of the mechanisms put into in play at a lubricated contact.Nevertheless, this postulate is also correlated with the electromechnical analysis of a lubricated contact, where electrical contacts between asperities are resistive and the electrical response of the fluid film is capacitive [21].The electrical circuit of an apparent contact behaves as a parallel RC circuit according to the lubricated regime.The dissipative properties of a contact are updated according to the lubrication regime.To do this, it was necessary to identify the fluid parameter Λ f as a discriminant.

Space and Time Discretization of the DEM Ball Bearing Model
To perform numerical simulations using DEM, we assumed that the friction coefficient µ(Λ f ) ranges over the interval [0.01, 0.1] as for the case of rolling/sliding concentrated contacts [43].Typically, µ = 0.1 for Λ f < 1, which means that the applied load is completely carried by the asperities of the contacting surfaces.When Λ f increases up to 1, µ decreases considerably passing from 0.1 to 0.01.In this range, the mixed regime is reached so that the applied load is carried by both lubricant film and asperities of the contacting surfaces.If Λ f becomes greater than 3 the fully fluid regime is established and the friction coefficient is kept constant, µ = 0.01.However, it should be mentioned that the friction coefficient is usually assumed constant in numerical simulations by DEM.Taking the fact that the lubricated contacts in a ball bearing must be considered individually, namely with their own friction coefficient and viscous damping, therefore, this can be achieved through an identification on the Strikeck curve by means of the fluid parameter Λ f .This can be done through a linear fitting made on the Stribeck curve (Figure 5).From a numerical standpoint, at a given time step of DEM simulation, the kinematic and mechanical conditions are known at each contact (rolling speed, relative speed, normal load), so that the lubrication regime may individually be determined through the fluid parameter expression (7).Hence, the friction coefficient µ(Λ f ) and the viscous damping C n (Λ f ) are calculated to carry out mechanical resolution for the next time step.
We note that regarding the direction t, the effect of the tangential damping C t , in the context of radial ball bearing, is assumed to be weaker than the normal damping and could be then neglected.Therefore, following the tangential direction the dissipation is basically related to the friction coefficient µ.
The contact stiffness model based on the elliptic integrals and damping model related to lubrication regimes described above, are implemented in the numerical tool using DEM.Therefore, Equation ( 1) is rewritten under the following form (8): In Equation ( 8), the coefficient of friction µ(Λ f ) is then considered as function of the fluid parameter Λ f .It should, however, be stated that the correlation of the operating conditions to the friction coefficient remains difficult to apprehend numerically.M. Björling et al. [44] studied numerically and experimentally the variation of the friction coefficient at different Slide-to-Roll ratios (SRRs) and drive speeds.Guegan et al. [42] also investigated friction and film thickness in elastohydrodynamic contacts through different Stribeck curves, where the authors achieved a dependency on the SRRs and drive speeds.In fact, and based on these results, it seems crucial to be aware of the dependence of friction coefficient on the operating conditions.
The DEM modelling is by solving Newton's second law for each rolling element.Equation ( 9) takes into account dynamical effects of rolling elements, such as centrifugal forces and gyroscopic effects, when the bearing is functioning.
where the subscript i denotes a given rolling element (or cage element).m i and I i are the mass and the quadratic moment of inertia.üi and θi are the linear and angular accelerations.F ext i and M ext i are the external force and moment acting on element i. F j− →i and M j− →i are the reaction force and moment interaction of the element j on the element i.
The mechanical resolution of ball bearing, composed of 2 × Z balls (rolling and cage elements), requires an explicit time integration based on the velocity-Verlet scheme (10).
where u, u are the displacement and velocity vectors, respectively, t the current time and ∆t the integration time step.The numerical simulations are conducted with ∆t of 1 µs.Choosing a large time step induces numerical instabilities and moreover it does not allow a steady mechanical state to be achieved.The effectiveness of the DEM lies in the fact that a discrete formulation allows each contact (ball-inner-raceway, ball-outer-raceway or ball-cage) to be easily tracked instantaneously when the ball bearing under operating conditions is simulated.Before studying the lubrication regime when the ball bearing operates under dynamical conditions, we first carried out a quasistatically loaded ball bearing to investigate the validity of the numerical response.The loading protocol of the ball bearing is of prime importance in performing realistic simulations in accordance with standard theoretical models for bearings as described in [34,45].Parameters related to the normal loading and radial deflection in ball bearing together with numerical results are introduced and discussed in Appendix A.

Numerical Prediction of Lubrication Regimes in Operating Radial Ball Bearing
As the relevance of the DEM is clearly established for accurately predicting the mechanical state in a radial ball bearing (Appendix A), we will now focus on the lubrication regimes under operating conditions.Ball bearings require fluid lubrication in order to perform satisfactorily for long periods time.Elastohydrodynamic lubrication (EHL) theory is concerned with the formation of a thin fluid film at the contact area between rolling elements and raceways [9,46].Under favourable conditions of mounting, speed rotation, loading and fluid viscosity, the lubricant film is assumed sufficiently thicker to separate the rolling surfaces.Obviously, the rolling bearings operating under a fully fluid regime have significant reduction of wear, and therefore, the expected bearing lifetime will be much longer than in the case of bearings operating in boundary lubrication regime.Generally, the minimum thickness of the fluid film is of the order of the roughness; nevertheless, in the case of optimal conditions at high speed and adequate viscosity, the fluid film exceeds the size of the surface asperities so that a fluid film thickness of several micrometers can be reached, 3 to 10 × σ a [47,48].
In this section, numerical simulations carried out with MULTICOR-3D software will be discussed.The lubrication regimes of a radial ball bearing under operating conditions are investigated by taking into account both imposed mechanical loading and fixed kinematic conditions.Sensitivity analyses are performed on the main driven parameters, namely the radial load F r , the diametral clearance P d and the angular speed ω of the inner ring, which more or less may affect the lubrication regime.

Effect of Radial Load and Diametral Clearance
The lubrication regime related to the film parameter Λ f (ψ) is identified on a given rolling element at each azimuth angle ψ(t) with respect to time t.A rolling element at an initial azimuth angle ψ 0 is considered and followed during the simulation.The tracked rolling element exactly reflects the lubrication regime of the radial ball bearing, from the starting phase, when the rotation is imposed on the inner ring, until the slowdown phase after going through the steady phase.The three phases are then studied to investigate the lubrication regime.
As an example, numerical simulation, aimed to check the efficiency of the discrete formulation, is initially conducted with the following conditions: ω = 500 rad•s −1 , P d = 0 (ε = 0.5) and F r = 3000 N.For the sake of conciseness, we have only been interested in lubrication regimes at the contact interface between rolling elements and the outer-raceway.Obviously, the same numerical analysis can be done for the rolling elements in contact with the inner-raceway.In the steady state, the theoretical expression for angular speed in pure rolling, is given by ω ball = ω×R i The SRR is calculated and plotted for a single rolling element in order to track the status of the rolling element/inner-receway contact, inside and outside the loaded zone of the ball bearing, at start-up, steady and slowdown phases (Figure 7b).The SRR is very low because the rolling element has a pure rolling regime when the steady state is reached, except for start-up and slowdown regimes.Obviously, the SRR of the rest of the rolling elements is just the same that for the plotted one.
The transition from one lubrication regime to another at the ball-outer-raceway contact is well depicted in Figure 8a,b.The plotted curves exhibit the fluid film thickness changing, starting from the start-up phase to the slowdown one, passing through the steady phase.The fluid film thickness h out min increases from zero to steady values (curve in royalblue), corresponding to the start-up phase, while for the slowdown phase (curve in dark-yellow), h out min decreases from steady values to zero.The steady regime plotted in red describes the transition of the rolling element from the loaded zone (low value of h out min ) to the unloaded zone (high value of h out min ).When the steady regime is reached the h out min (curve in red) is still the same for each contact occurring between rolling element and outer-raceway, inside and outside the loaded zone.For these operating conditions for both start-up and slowdown phases the lubrication regime in the ball bearing changes according to the fluid parameter Λ f , so the value ranging of Λ f is between 0 to 3.However, the ball bearing only operates in mixed regime when the steady phase is reached, with Λ f varying between 2.15 and 3.
On the other hand, in Figure 8a,b at the steady phase we observe the influence of the loaded zone on the fluid film thickness and the fluid parameter, h out min and Λ out f , respectively.The fluid film thickness decreases for ψ ranging over the interval − π 2 , π 2 in accordance with the load parameter ε = 0.5.The minimum value of h out min is reached at ψ = 0, (h out min < 1 µm), in accordance with the applied radial load F r .Consequently, the numerical simulation using DEM may be considered as an efficient predictive tool for monitoring the transition between lubrication regimes in ball bearings and in operating conditions.
Secondly, we are now interested in the radial load effect on the lubrication regime in a ball bearing.The operating conditions and the diametral clearance are kept the same, with ω = 500 rad•s −1 and P d = 0 (ε = 0.5).We have only considered the steady phase and lubricated ball-outer-raceway contact for each simulated radial load F r .We assumed that both outer and inner contacts achieve similar lubrication regimes in the loaded zone of ball bearing for a fixed angular speed.Numerical simulations have carried out for six values of F r ranging in the interval [750 N, 7500 N].In Figure 9a is plotted the fluid parameter for all the simulated radial loads with respect to fixed operating conditions.It turns out that F r does not play a significant role because the lubrication of the ball bearing is still in mixed lubrication regime the radial load.This is well observed in Figure 9b on a logarithmic scale, where the fluid parameter Λ out f decreases slightly with increasing maximum normal load Q max as function of the radial load F r .In addition, the fit made on the curve of the fluid parameter Λ out f , as a function of the maximum normal load Q max , leads to a linear regression with a negative slope a 1 of −0.0727.This value is expected since it is related to the EHL model through the dimensionless parameter W raised to the power −0.073.Also, the second fit made on the curve depicted the variation of F r with respect to Q max leads to a slope a 2 very close to 1, which is in good agreement with Equation (A3).In view of the curves plotted in Figure 10a of the fluid parameter Λ out f , the diametral clearance P d does not notably affect the lubrication regime.In addition, Λ out f clearly describes a plateau at azimuth angle ψ = 0 (Figure 10b).However, the fluid parameter distribution differs considerably in shape when the diametral clearance varies, passing from positive values to negative values.Indeed, passing from a negative value of P d to a positive one, it directly affects the fluid parameter with respect to azimuth angle ψ of rolling element.The loaded zone of the ball bearing with a diametral preload (P d < 0, ε > 0.5) is larger than the case with a diametral clearance (P d > 0, ε < 0.5) (Figure 10a).Inversely, the maximum normal load Q max increases when a diametral clearance is considered (Figure 10b).In such a mechanical configuration, we have an overloading in the ball bearing with a reduction of the loaded zone size.Therefore, the ball bearing operates under severe conditions that could potentially induce early wear mechanisms because of the concentration of normal load Q max in the loaded zone.Moreover, the manufacturers advise users to avoid the case with a load parameter ε < 0.5 because of the risk of premature degradation of ball bearings.It is admitted by the manufacturers of bearings, that the lifetime of the components is optimized, provided that the mounting is done with a zero clearance (P d = 0, ε = 0.5) and a sufficient lubricant amount.One can state that the radial load F r , the diametral clearance P d > 0 and the diametral preload P d < 0 do not have a significant effect on the lubrication regime.As shown by the previous numerical simulations, when the operating steady phase is reached in the ball bearing, no change is expected in the established EHL regime.

Effect of Angular Speed
We are now interested in the influence of the shaft angular speed on the lubrication regime of the radial ball bearing.The angular speed ω (rad•s −1 ) ranges over the interval [100, 1700] and varies in steps of 200.A radial load F r of 3000 N and a zero clearance (P d = 0, ε = 0.5) are considered.The lubrication regimes are investigated for both ball-inner-raceway and ball-outer-raceway contacts.We anticipate that the lubrication regime will be more sensitive to the angular speed of the shaft this time as well because the lubricant gains in lift.The viscosity and the rolling speed have more significant effects on the minimum fluid film thickness h min as indicated by Equation (6).We note that the rolling velocity U r is for pure rolling; the relative sliding velocity is assumed negligible.However, if the sliding velocity becomes significant then the dimensionless rolling velocity U r is expressed as a function of both rolling and sliding velocities [37].
In addition, the fluid film thickness at the contact with the outer-raceway can be estimated at each revolution of the inner ring.Indeed, because of centrifugal effects the contacts are persistent between the rolling elements and the outer ring (Figure 11a).Furthermore, a disturbance of the fluid film thickness is observed at ω = 100 rad•s −1 in particular when the rolling element is outside the loaded zone.This effect is mainly due to changes in rolling velocity of the ball.However, when the angular speed of the inner ring becomes higher and higher, allowing the rolling element to reach a more steady velocity, the disturbance totally disappears.On the other hand, because of the load parameter ε = 0.5, corresponding to zero clearance (P d = 0), the fluid film thickness at the contact with the inner-raceway can only be predicted within the loaded zone varying in the interval − π 2 , π 2 (Figure 11b).Hence, we have loss of contact with the rolling element in the unloaded zone.Obviously, the latter is kept under a lubrication regime due to the switching with the loaded zone even if the fluid parameter Λ inn f is not predicted numerically.The numerical simulations conducted on radial ball bearings allow us to distinguish between the lubrication regimes for both ball-outer-raceway and ball-inner-raceway contacts.As is well depicted in Figure 11a,b, the film parameter Λ inn f related to the fluid film thickness at the ball-inner-raceway contact is lower because the equivalent radius of convex curvatures is lower than that of concave curvatures . Indeed, the comparison of the inner and the outer dimensionless loads W inn,out , upon which depends Equation (6), indicates a higher value for the ball-inner-raceway contact (Figure 12a).These results in higher contact forces act on the fluid film thickness and lead to a maximum pressure at the convex contact with the inner-raceway.Hence, it is sufficient to estimate the fluid film thickness h inn min at the contact with the inner-raceway to predict the lubrication regime.Nevertheless, for a ball bearing operating at high speed, the centrifugal force of the rolling element is added to the normal force at the contact with the outer-raceway.The fluid film thickness h out min may be only estimated at the contact with the outer-raceway.This centrifugal effect is clearly observed in Figure 12b at the ball-outer-raceway contact where the maximum normal load Q max increases in comparison to the ball-inner-raceway contact, where Q max remains close to the analytic solution in static mode, with Q th max = 1457 N.Moreover, the fit made on the curve of the fluid parameter Λ inn f , as function of the angular speed of shaft ω (Figure 11a), achieves a linear regression with a positive slope a of 0.6797.This value is very close to the exponent of the dimensionless term of the EHL model, related to the rolling velocity U 0.68 r .This confirms once again that the DEM approach is well done, and above all allows fast monitoring of the lubrication regime of the ball bearing under operating conditions.In Table 3 are summarized the lubrication regimes with respect to imposed angular speed.We remark when the ball bearing operates at a specific angular speed ω the rolling elements may simultaneously be under two lubrication regimes.For example, in the case of ω in the interval [900, 1700], we have at the same time both mixed and fully fluid regimes at the inner-raceway and outer-raceway contacts, respectively (Table 3).Unlike the cases when radial load F r and diametral clearance P d are considered to study their influence on the lubrication regime, the imposed angular speed plays a more significant role in lubrication regime transition.The EHL model is more sensitive to the rolling speed since the dimensionless rolling velocity U 0.68 r prevails over the dimensionless load W −0.073 .

Conclusions and Perspectives
Through this study we wanted to highlight the relevance of the DEM approach in describing the dynamic behaviour of multicontact systems such as rolling bearings.Under a quasi-static radial load, the proposed discrete description of a radial ball bearing, using a contact stiffness model based on elliptic integrals, achieves very similar radial deflection and normal load distributions in comparison to analytic formulation [45].The implementation of the discrete model in the MULTICOR-3D software allowed us to investigate the lubrication regimes from a qualitative point of view under operating conditions.Moreover, in the context of the EHL study, the DEM has proven to be a useful numerical tool since it enables continuous prediction of the lubrication regimes and dynamic management of the contact interface.By driving various parameters to be investigated, such as the radial load, diametral clearance and shaft angular speed, it has clearly been shown that the lubrication regime is basically dependent on the angular speed of the inner ring.In other words, the transition between the different lubrication regimes is mainly related to the rolling speed at the lubricated contact in the loaded zone of the ball bearing.As shown in figures consecutive to the parametric study, once the ball bearing reaches the steady state, both the radial load and diametral clearance have a minor effect on the lubrication regime.Starting from the mechanical behaviour at the contact interface and lubrication regime identification, the proposed discrete modelling may be extended to further studies, for example, dealing with a ball bearing operating under abnormal load or noisy kinematic conditions.In fact, such undesirable defects appear often when ball bearings are subjected to extreme operating conditions, typically under heavy loadings, in the presence of severe shocks, or in case of inadequate mounting.Furthermore, the innovative monitoring method recently developed and based on electrical measurements [20] may be also performed for ball bearing diagnosis taking into consideration the capacitive behaviour or the dielectric properties of the lubricant.In addition, thermal effects due to the operating conditions of the rolling bearing involve significant modifications of the lubricant properties, in terms of viscosity and contamination by debris.These thermal aspects must also be considered in the discrete modelling to investigate the lubricant aging effects on the lubrication regime.In Figure A2a,b are plotted both radial deflection δ ψ and normal load Q ψ for radial load F r starting at 1000 N and increasing in steps of 1000 N up to 3000 N. The numerical simulations achieve a mechanical state in the radial ball bearing in good agreement with analytic formulation.The more radial load F r there is, the larger the radial deflection δ ψ and normal load Q ψ .To prove the efficiency of the DEM approach, a simple assessment can be made on δ ψ and Q ψ at ψ = ± π 2 with zero clearance (P d = 0, ε = 0.5).In accordance with Equations (A1) and (A2), the numerical simulations lead to

Figure 2 .
Figure 2. Contact modelling in smoothed formulation of DEM.

Figure 3 .
Figure 3. Curvatures in contact in two orthogonal cross sections of a ball bearing: (a) x-z plane ; (b) y-z plane.

Figure 4 .
Figure 4. Radial shift δ r due to the imposed radial load F r .

Figure 5 .
Figure 5. Schematic representation of Stribeck curve: friction coefficient and fluid film thickness as function of fluid parameter.

Figure 6 .
Figure 6.The three status of viscous damping: hysteretic, mixed and fluid.

.Figure 7 .
Figure 7. (a) Imposed shaft angular speed ω and simulated angular speed ω ball of rolling element and (b) SRR for a single rolling element.

Figure 8 .
Figure 8.(a) Film fluid thickness h out min and (b) fluid parameter Λ out f with respect to azimuth angle ψ of the ball-outer-raceway contact.

Figure 9 .
Figure 9. (a) Fluid parameter Λ out f as function of azimuth angle ψ; (b) Fluid parameter Λ out f and maximum normal load Q max in the loaded zone at ψ = 0.

Figure 10 .
Figure 10.(a) Fluid parameter Λ out f as a function of azimuth angle ψ; (b) Fluid parameter Λ out f and load parameter ε as function of normal load Q max at ψ = 0.

Figure 11 .
Figure 11.Fluid parameter as function of azimuth angle ψ: (a) Λ out f for ball-outer-raceway contact; (b) Λ inn f for ball-inner-raceway contact.

Figure 12 .
Figure 12.(a) Fluid parameter and maximum contact pressure at ψ = 0 as function of angular speed ω; (b) Minimum fluid film thickness and maximum normal load at ψ = 0 as function of angular speed ω

Figure A2 .
Figure A2.(a) Radial deflection δ ψ and (b) normal load Q ψ distributions under several radial load F r .

Table 1 .
Geometrical characteristics of radial ball bearing of 6208 series.

Table 2 .
Identification of lubrication regimes for non-conformal contacts.