Digital twin by DEM for ball bearing operating under EHD conditions

In this paper we investigate the elastohydrodynamic (EHD) lubrication in ball bearings by means of a digital twin. Based on an original description involving the discrete element method (DEM), the digital twin integrates all the components of ball bearings and enables realistic behavior under mechanical loading and kinematic conditions. In order to check the standard indicators recommended by most ball bearing manufactures, a stiffness model for elliptical Hertzian contact and an improved EHD formulation for lubricated contact are implemented in the numerical tool. In addition, we have introduced into the discrete modeling an electrical capacitance model correlated to the fluid film thickness and the contact pressure. The numerical predictions of lubricant film capacitance provided by digital twin are in good accordance, both qualitatively and quantitatively, with the experimental data available in the literature. The coupling of the discrete method with the electrical approach enables efficient solutions to be provided in terms of lubrication regime in relation to the lubricant properties to optimize the bearing lifetime.


Introduction
Industrial machines require a good lubrication regime with an appropriate lubricant to avoid wear mechanisms that could occur at the contact interface between two surfaces. In rotating machines, several components such as ball bearings, cams and gears operate under EHD lubrication regimes with non-conformal contacts. The nonconformity of the surfaces in contact means that the contact areas are subjected to a high contact pressure which can be in the range of several gigapascals. With a view to being able to assess whether the machine elements are working with safety conditions or are close to reaching boundary lubrication and ultimate failure, the fluid film thickness is consequently of a prime importance. Such mechanical components are often designed for the lubricant film to be thick enough, more than 1 µm, to avoid contact between the mating surfaces ensuring the expected bearing lifetime. For that purpose several experimental methods are dedicated to measuring the lubricant film thickness. The most widely used one is based on optical interferometry [1][2][3]. Nevertheless, the electrical technique has previously been used with measurements of electrical resistance [4][5][6] and capacitance [7][8][9]. The electrical resistance can * e-mail: mohamed.guessasma@u-picardie.fr indeed be useful when the contact between two surfaces is metallic which limits its use as a conducting contact interface. The electrical capacitance, on the other hand, is specifically linked to the contact interface thickness that accommodates the lubricant film. In the field of ball bearing diagnosis, the electrical techniques have been used to measure the load distribution by capacitance sensors [10] as well as to highlight the angular speed effect on the lubrication regime by means of the electrical resistance [11]. More recently, Jablonka et al. [12] proposed a quantitative evaluation of lubricant film thickness in a radially loaded deep-groove ball bearing using electrical capacitance measurements. The authors used in-situ monitoring tool to identify the lubrication conditions as function to the operating parameters and lubricant properties.
Regarding numerical modeling in the field of ball bearing there are basically two approaches to perform numerical computations, namely the finite element method (FEM) and the multibody method. Numerical tools allow analysis on specific ball bearing parameters, such as kinematics, dynamics, stresses, deformations, lubrication regime and lifetime. Because of the cyclic loadings, bearing rings are the most critical components since the flaking phenomena often appearing on rolling raceways is predictive of the onset of the damage mechanisms. In order to better understand the fatigue behavior of bearing rings, FEM was used to quantify the stress level and its distribution within the material according to the applied mechanical loading [13]. From a computational time point of view, however, FEM modeling remains costly due to the numerical management of the contact. Furthermore, FEM is often performed to describe only a portion of ball bearing. Unlike FEM, the multibody method offers the advantage of allowing a description taking into account the entire bearing, often with the exception of the cage, for modeling the loading distribution and lubrication effects [14,15]. Recently DEM has been used successfully for simulating the dynamic behavior of bearings with tribology enrichments [16] as well as for modeling lubrication regimes by means of an EHD model [11,17]. To achieve this, the digital twin of ball bearings implements a major part of the functionalities according to the state of the art, such as the slice model, the mutual influences of neighboring elements with clearance, elliptical Hertzian contact model, lubricant film estimation, damping forces, centrifugal forces, etc. In the following paper, the ball bearing is simulated with the assumption of EHD lubrication conditions. The coupling between the mechanical pressure and the dielectric behavior of the lubricant film thickness is carried out by means of the electrical capacitance. According to the model of a two parallel-plate capacitor separated by a dielectric material, the capacitance of the lubricant film thickness is extracted numerically. The influence of lubricants properties (viscosity and dielectric constant) on the capacitance is well highlighted. The numerical predictions of fluid film capacitance for the two simulated lubricants provide a promising result in comparison to the experimental data available in the bibliography [12].
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 modeling performed on a ball bearing, with a digital twin called MULTICOR3D. The simulated radial ball bearing of 6208 series is designed fully and accurately in accordance with its geometrical conformity. The third section focuses on the mechanical modeling of contact forces based on a smoothed formulation of DEM in the case of an elliptical Hertzian contact with the assumption of piezo-viscous-elastic behavior of the lubricant film. The electrical capacitance of the lubricant is then introduced by means of the parallel-plate capacitor formula depending on the contact pressure and fluid film thickness. The next section deals with the numerical simulations where the effects of operating conditions on the fluid film capacitance are highlighted and discussed. Some comparisons with experimental measurements extracted from the bibliography are also discussed in this section. In the conclusion, we summarize the contributions of this work and we end with a few thoughts on future work.

Digital twin of radial ball bearing
A three-dimensional digital twin of ball bearing (Fig. 1a) is developed by means of MULTICOR3D software using DEM [17]. The radial ball bearing of 6208 series (SKF  6208 TN9/DBGA) has been considered in the present study (Fig. 1b).
The dimensions and geometrical conformity of the 6208 radial ball bearing are provided in Table 1. The number of rolling elements which are likely to be in contact with the inner or outer raceways is Z = 9. Furthermore, nine additional smaller balls of radius R z in contact with two toruses are considered to reproduce the cage effects ( Fig. 2). 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 [17].
We note that this description of the cage mainly has two roles. The first is that it ensures a conforming 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 balls may be occur with rolling elements. By considering all bearing components, the proposed modeling gives the opportunity to develop a digital twin without any restriction.
We note that the digital twin design performed by DEM allowed us to numerically achieve standard bearing frequencies [18]. For validation purpose, normal load distribution for several diametral clearances and free contact angle variation as a function of the raceway conformities are shown in Appendix A.

Mechanical modeling of contact forces
The contact stiffness model that we used is based on a smoothed formulation of DEM initially proposed by Cundall and Strack [19]. The contact interactions are obtained by means of Kelvin-Voigt Spring-Dashpot Model and a Coulomb's law if sliding occurs at the contact (Fig. 3). The main assumption made with DEM is that the particles in contact are assumed perfectly rigid and therefore, there is no elastic deformation of all bearing components. The elasticity is only considered at the contact between, rolling element/raceways (inner or outer), cage ball/toruses and rolling element/cage ball (Fig. 2).

Hertzian contact stiffness model
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 (1).
with K n,t the normal and tangential stiffnesses, C n,t the normal and tangential viscous damping coefficients, δ n,t the normal and tangential relative displacements, 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 mechanical steady state is a function of two damping types, namely hysteretic and EHD, related to the contact shape, roughness, fluid viscosity, temperature, entrainment speed, etc. The tangential viscous damping coefficient C t is more difficult to evaluate. The ratio Ct Cn is usually considered to range within the interval [0, 1]. The normal contact stiffness K n is related to the mechanical characteristics and curvatures of the surfaces in contact according to non-linear model; F n ∝ δ 3 /2 n . For a given rolling element in contact with the inner-raceway (or outer-raceway), with the assumption of an elliptical Hertzian contact, the normal stiffness is written as function of the curvature sum radius and the approximate elliptic integrals as introduced by Hamrock [20]: R2y−R1y (2)  where R curve is the curvature sum radius, R x the effective radius in x direction, R y the effective radius in y direction, κ the ellipticity parameter κ = α 2 π , α = Ry Rx , 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 α . The equivalent radius R x is calculated by taking into account the cases of convex contact (rolling element/inner-raceway) and concave contact (rolling element/outer-raceway), as given in equation (2). For the sake of clarity, we note that (Fig. 4).
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 without friction and the damping coefficient is considered critical. The normal stiffness (3) only is taken into account, and therefore the following equation is considered [21]: the equivalent radius of rolling element and cage ball in contact of radius r i and r j , respectively. To achieve a steady regime of the ball bearing under operating conditions, viscous damping and local friction are added as dissipative forces to the mechanical model (1). When a contact occurs between rolling element and raceways (inner or outer) two kinds of normal viscous damping could be involved: hysteritic and(or) fluid. The viscous effect due to hysteretic damping, which is related to the plastic strain of surface roughness (4), is expressed by the coefficient C hyst in the normal direction n.
where α e denotes the restitution coefficient of steel [22]. A value of 0.08 s.m −1 is considered to compute the hysteritic damping C hyst .
Regarding the viscous damping due to the fluid film, Walford and Stone [23] established a damping coefficient C fluid related to the lubrication film in the normal direction of contact (5). Since the bearing contains nonconformal contacts endowed with EHD lubrication regime, the lubricant film must be thick enough to avoid contact between the mating surfaces.
where η is the dynamic viscosity of the fluid at operating temperature and atmospheric pressure [24], a the semi-major axis of the Hertz elliptical contact [25], h min the minimum lubricant film thickness of the ball-innerraceway (or -outer-raceway) contact. The DEM modeling is by solving Newton's second law for each rolling element taking into account dynamic effects, such as centrifugal forces and gyroscopic effects, when the bearing is functioning. The mechanical resolution requires an explicit time integration based on the velocity-Verlet scheme [17].

EHD lubrication model
In the field of EHD lubrication, oil (or grease) rheology influences the fluid film thickness under operating conditions [26]. Given that the viscosity of a lubricant increases with respect to pressure, this means that the lubricant exhibits piezo-viscous behavior. Furthermore, temperature also plays a role in lubricant viscosity changes with the consequence of diminishing its viscosity. Starting from these assumptions with a piezo effect of the lubricant, we have implemented an EHD model in the digital twin with piezo-viscous-elastic regime to predict the minimum fluid film thickness h min . As initially proposed by Hamrock and Anderson [20], the general form of the EHD model is a function of parameters related to the speed, load and material properties. Starting from the previous studies in the field of EHD lubrication, Masjedi and Khonsari [27] improved relatively recently the EHD model by introducing the ellipticity effect. Since changing bearing conformity affects the contact ellipticity, we have implemented the new EHD formulation in the digital twin in which the parameters of the model are a function of the ellipticity. In the case of smooth point-contact EHD, the expression of the minimum fluid film thickness (6) is as follows: where E, U r , W are three dimensionless parameters: E = ξE , with ξ the pressure-viscosity coefficient of the lubricant at operating temperature and atmospheric pressure, , with Q ψ the ball-inner-raceway (or -outer-raceway) normal load at an azimuth angle ψ (Fig. A.1).
It is commonly agreed that the identification of the lubrication regime in which a system is running is by the use of the dimensionless fluid parameter Λ f such defined by (7). Λ f is the ratio of the minimum fluid film thickness to the surface roughness [28]: where σ 1 and σ 2 denote the standard mean roughnesses of the measured surfaces for the inner-raceway (or outerraceway) and rolling elements, respectively. However, it is admitted that only the raceway roughness is considered, and the rolling element roughness can then be neglected (σ 2 σ 1 ). In addition, the schematic representation of the Stribeck curve (Fig. 5) for non-conformal contact shows the dependence of the friction coefficient µ on the dimensionless fluid parameter Λ f . The friction coefficient µ is therefore dependent on the lubrication regime [29]. In addition, the Stribeck curve mainly defines three lubrication regimes (boundary, mixed-film or fully-fluid film) in relation to fluid parameter Λ f .

Electrical capacitance model for lubricated Hertzian contact
In the context of the identification of a lubrication film by means of electrical measurements, the fluid film thickness can easily be extracted from the capacitance of the contact [30]. According to the parallel-plate formula the lubricant film thickness is inversely proportional to the capacitance (8). with C Hertz the Hertzian contact capacitance between two parallel conducting plates of area A Hertz separated by a dielectric material where ε r is the dielectric constant of the fluid and h c its central thickness. ε 0 is the vacuum permitivity ε 0 = 8.854 pF.m −1 . The theoretical form of the contact area between the rolling element and raceway for a perfect elliptical Hertzian contact is A Hertz = πab, with a and b the major and minor semi-axes of the ellipse, respectively. In the case of a smooth surface [27] the central film fluid thickness h c integrating the ellepticity effect is as follows (9): From experimental stand point, the qualitative conversion of the measured capacitance to lubricant film thickness involves first the estimation of the capacitance related to the nearby area outside the Hertzian contact and some background capacitance. The latter two are then subtracted from all measurement data. The remaining value corresponds to the contact capacitance due to the fluid film thickness [31]. From a numerical point of view, the contact capacitance is easily computed from expression (8) after predicting the central fluid film thickness by (9).
The dielectric constant ε r of the lubricants varies with Hertz contact pressure. For nonpolar lubricants ε r can be calculated from Clausius-Mossotti model [32]: where N A is the Avogadro number, β the molecular polarizability, M the molecular weight and ρ the density. The factor N A β 3M can be deduced by considering permittivity and density at atmospheric pressure. According to the formula established by Dowson and Higginson [33], the variation of the lubricant density as a function of the contact pressure is given by following equation (11): ρ = ρ 0 0.59 · 10 9 + 1.34 · P Hertz 0.59 · 10 9 + P Hertz (11) with ρ 0 the ambient-pressure density and P Hertz the Hertz contact pressure. In Figure 6a are plotted the dielectric constant evolution of synthetic oil film (Tab. 2) with respect to Hertz contact pressure at inner and outer raceways. This can easily be understood by the fact that the higher a pressure, the greater the dielectric constant of the lubricant film. As depicted by the curves in Figure 6b, the maximum Hertz contact pressure at the inner raceway is greater than reached at the outer raceway due to the difference of the ellepticity of both inner and outer raceways. The contact pressures are simulated under dynamic operating conditions, with radial force F r = 5 kN and shaft angular speed ω = 477.5 rpm. From electrical point of view, a single rolling element in contact with raceways can be considered as an equivalent circuit composed of capacitances between the inner and outer raceway contacts, which are in series. It should be pointed out that the cage components (PA66) are insulating which allows the ball bearing to be considered as a set of capacitances in parallel (Fig. 7). The total capacitance of the ball bearing is calculated by summing the equivalent contact capacitances for all rolling elements located in the loaded zone.
Hertz × C outer,i Hertz C inner,i Hertz + C outer,i Hertz (12) with N l the number of rolling elements in the loaded zone of the ball bearing. The contact capacitances C inner,i Hertz and C outer,i Hertz between the ith rolling element and inner-outer raceways, respectively, are determined according to the parallel-plate formula (8). It must be emphasized that, unlike the case of experimental measurements, the capacitance model implemented in the digital twin does not take into consideration the contribution of capacitance between inner and outer ring and some background capacitance [31,34].

Numerical simulations with the digital twin of ball bearing
This section is dedicated to numerical simulations performed on the digital twin of ball bearing. The correlation between contact capacitance and fluid film thickness is investigated when the ball bearing operates under 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 shaft on which the inner ring is mounted. We recall that for benchmarking purpose some numerical simulations considering radial loading, free contact angle and free endplay in ball bearing have been carried out and given in Appendix A.

Effects of radial load and diametral clearance on electrical capacitance
According to the parallel-plate capacitor formula the capacitance of the film fluid thickness is extracted by considering a single rolling element after one revolution.
Simulations are then carried out in order to exhibit the influence of the radial load F r and diametral clearance P d on the capacitance of the fluid film thickness. Comparisons with respect to experimental measurements [12] are performed in order to highlight the performances of the digital twin of ball bearing. Two kinds of lubricants are considered based on the experimental study given in [12], namely mineral oil and polyalphaolefin synthetic oil (PAO). In Table 2 are summarized the physical properties of both lubricants. Figure 8a depicts the sensitivity of the capacitance to the radial load F r ranging over the interval [1−7 kN] in the case of mineral oil. The capacitance is proportional to F r , which means that the higher the radial load F r , the greater the change in fluid film capacitance. From a numerical point of view, when the rolling element loses contact with the rings of the ball bearing, the capacitances in series related to the inner and outer raceways are disabled (Fig. 7). In addition, as can be observed in Figure 8a, the predicted capacitance for F r = 5 kN is   in quite good agreement with the experimental measurements obtained with the same entrainment speed in the case of 6306 ETN9 deep groove bearing with a polymer cage [12].
For comparison purpose we have added to our numerical results the capacitance related to the background and the capacitance due to the space between inner and outer raceways. The equivalent capacitance due to these latter effects is about 30 pF according to the experimental results given in [12]. Therefore, the additional capacitance is assumed to be of the same order of magnitude as for the simulated ball bearing 6208 series (SKF 6208 TN9/DBGA). Good agreement with experimental Table 3. Simulated clearances P d and corresponding load parameter ε. measurements can also be observed on the curve fit applied to the numerical predictions for a zero azimuth angle ψ (Fig. 8b) where the radial load reaches its maximum. This confirms that the digital twin of ball bearing is able to predict fluid film capacitance fairly close to experimental data. The next simulation deals with the effect of the load parameter ε on the capacitance distribution in the ball bearing by considering several values of P d D ratio. The capacitance distribution in the loaded zone is plotted for three configurations, typically with negative clearance, zero clearance and positive clearance. The values of P d D ratio that have been considered are summarized in Table 3.
The effect of the diametral clearance P d is well highlighted in Figure 9. Because there are more contacts between elements and raceways in the case of P d < 0, this increases the number of capacitances operating in the electrical model (Fig. 7). As can be observed in Figure 9, for P d D ratio of −0.0035, corresponding to a load parameter ε of 1.32, the capacitances of the electrical model are all activated. Conversely, in the case of positive clearance more than 50% of the capacitances are disabled because of the reduced size of loaded zone (Fig. 9).

Capacitance variation with respect to entrainment speed and fluid film thickness
In this part dedicated to the study of capacitance variation as a function of the entrainment speed U r and the fluid film thickness h min , we have considered two radial loads F r = 1 kN and F r = 5 kN. The diametral clearance P d is kept constant and equal to zero (ε = 0.5). First of all, we have plotted the capacitance variation with respect to U r for the two lubricants listed in Table 2. We note a sharp variation in capacitance for low rolling speed under a given radial load F r . Clearly, there is a strong correlation between the capacitance and the contact area (F r ) on the one hand and fluid film thickness (U r ) on the other. In addition, an increase in U r carries a slight decrease in the capacitance ( Fig. 10a and b). We note that the capacitance values for the synthetic lubricant are higher then those for the mineral lubricant. The difference between the electrical behaviors of the two lubricants is mainly driven by the chemical composition of each oil. According to the experimental measurements performed on PAO6 (synthetic oil [35]), as expected, is characterized by a very low electrical conductivity, while mineral oil is a rather better electrical conductor.
Once again, we have carried out a comparison of the numerical simulations with the measurements resulting from the experimental study performed with 6306 ETN9 deep groove bearing [12]. According to the previous capacitance distributions (Fig. 8a) plotted as function of the angular position of one rolling element, the capacitance value at an azimuth angle ψ = 0, corresponding to the maximum normal load Q ψ , provides more valuable indications about the lubrication regime. Figure 10a for mineral oil and Figure 10b for PAO VG48 exhibit quite a good agreement between numerical predictions and the experimental results. In the case of radial load F r = 1 kN the capacitance predictions as function of entrainment speed are quite close to experimental one for both lubricants. However, for F r = 5 kN, there are some discrepancies with experimental measurements in Fig. 11. Capacitance-dielectric constant ratio as function of fluid film thickness of the element-inner raceway contact. particular for PAO VG48 lubricant (Fig. 10b). This is related in part to the fluid film thickness difference at the contact between the rolling elements and raceways. In fact, the asperities effect is more significant at a low entrainment speed U r combined with a high radial load F r , which explains why the measured capacitance is higher than the numerical capacitance. In [27] the authors have suggested a parameter (asperity load ratio) in order to quantify the load carried by the asperities. Because of a better electrical metal conductivity in comparison to lubricant, the higher the asperity load ratio, the higher the contact capacitance. Furthermore, the effects of the outside zones close to the Hertzian contact [36] have to be taken into consideration for a better agreement of numerical predictions of contact capacitance with experimental results. As mentioned in [12], the background and ring capacitances are subtracted from all experimental measurements.
To avoid the effect of the lubricant properties, Figure 11 depicts the overall capacitance of the ball bearing as a function of the fluid film thickness at the contact between rolling element and inner raceway. Only the lubricant film thickness at the inner raceway has been considered because of a higher contact pressure. Assuming that the two radial loads provide curves with similar trends, only the case of F r = 5 kN is simulated. As for the previous capacitance curves (Fig. 10a and b) a sharp decrease in capacitance is observed for low fluid film thickness, which corresponds to low entrainment speed U r . It stands to reason that from a numerical point of view, it is easy to correlate the capacitance to the fluid film thickness by means of equation (6). In this sense, the digital twin gives benefit over an experimental measurement because it easily allows the capacitance change to be followed according to the lubrication regime when the bearing is operating. Figure 12a and b depict the dependence of fluid parameter Λ f on entrainment speed U r and the correlation between Λ f and the overall capacitance of the ball bearing. In the interval of entrainment speeds displayed on the x-axis and under a radial load F r = 5 kN, the lubrication regime in the ball bearing skips from boundary to mixed according to the Stribeck curve (Fig. 5). For a given entrainment speed U r , we clearly observe a difference in lubrication regime when the ball bearing operates with mineral oil and when it operates with synthetic lubricant PAO (Fig. 12a). Because the two tested lubricants have different viscosities, the ball bearing rapidly reaches the mixed regime with mineral oil in comparison to the PAO. This helps in predicting suitable operating conditions with an appropriate lubricant in order to increase the ball bearing lifetime. In Figure 12b the intervals of boundary and mixed regimes are plotted with respect to the ball bearing capacitance. For the considered ball bearing and radial loading, the change from boundary to mixed regime occurs when capacitance reaches 60 pF.

Conclusion
Throughout this work we have attempted to highlight the relevance of the digital twin when a ball bearing operates under EHD lubrication regimes. A mechanical modeling based on DEM has been carried out with the aim of providing the realistic behavior of a ball bearing in terms of kinematic parameters, contact stiffness and lubricant film effects. By means of a parallel-plate capacitor formula the capacitance of the fluid film thickness is numerically extracted with respect to the EHD lubrication regimes. Furthermore, qualitative and quantitative comparisons of the numerical predictions with experimental data highlight the interest of the digital twin for industrial applications. The digital twin could indeed be a useful predictive tool for recommended lubrication regimes in bearings with a view to enhancing the performance of industrial machines. Nevertheless, the proposed modeling could be improved in order to take into account the effect of asperities in the EHD formulation.

Appendix A: Normal load and axial free endplay
According to the standard theory, the mechanical state of a ball bearing under a radial load F r is provided by a dimensionless load parameter, ε = 1 2 1 − P d 2δr , related to the diametral clearance P d and the radial shift δ r of the inner raceway (Fig. A.1).
The normal load Q ψ due to radial load F r at an azimuth angle ψ is function of the load parameter ε such that [37]: with Q max the maximum normal load at ψ = 0. Stribeck [29] proposed that Q max is function of the radial integral J r and the number of rolling elements Z Q max = Fr ZJr(ε) . Figure A.2 exhibits the normal load Q ψ plotted with respect to rolling element position at an azimuth angle ψ. Three values of load parameter ε are considered, 1.7, 0.5 and 0.29, which corresponds to negative clearance P d D = −0.0045 , zero clearance P d D = 0 and positive clearance P d D = 0.0045 , respectively. The azimuth angle ψ of the rolling element is varied by driving the inner ring at a slow angular speed ω = 1 rad × s −1 to avoid any centrifugal effect. One can conclude that the digital twin of ball bearing achieves results in terms of radial load distribution very close to the analytic solution given in [37].
The next comparison with the analytic solution was made on the free contact angle and free endplay. In fact a  radial ball bearing is commonly designed with a diametral clearance which causes both radial and axial degrees of freedom. Under a free axial thrust the rolling element and raceway contact become oblique [37]. Figure A.3 shows the effect of axial play on radial ball bearing when the inner ring moves to the right. The distance between the centers O i and O o of the inner and outer raceway curvatures respectively, is: with B the total curvature of the bearing. We note that the numerical simulations have been carried out with identical conformities of inner and outer races which yields total curvature: B = can express the free contact angle α 0 made by the vertical line and the line passing through the points of contact between the rolling element and both raceways as follows: As mentioned before, because of diametral clearance P d the bearing is free to move axially without mechanical loading. One can define the free endplay as the maximum relative axial displacement of the outer raceway with respect to the inner raceway:   As shown in Figure A.4, the simulations achieve results in good agreement with analytical model [37]. In conclusion, the benchmarking performed between the theoretical formulation and numerical simulations shows that DEM is well adapted to ball bearing modeling and provides accurate results.