Issue 
Mechanics & Industry
Volume 21, Number 5, 2020
Scientific challenges and industrial applications in mechanical engineering



Article Number  506  
Number of page(s)  12  
DOI  https://doi.org/10.1051/meca/2020022  
Published online  10 August 2020 
Regular Article
Digital twin by DEM for ball bearing operating under EHD conditions
^{1}
Campus Universitaire  LTIEA3899,
48 rue d’Ostende,
02315
SaintQuentin,
France
^{2}
Université de Djibouti,
avenue Djanaleh
BP 1904, Djibouti
^{*} email: mohamed.guessasma@upicardie.fr
Received:
23
October
2019
Accepted:
11
March
2020
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.
Key words: Ball bearings / DEM / lubrication / EHD / electrical capacitance / friction / Stribeck curve
© S.O. Farah et al., published by EDP Sciences 2020
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 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 nonconformal 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–3]. Nevertheless, the electrical technique has previously been used with measurements of electrical resistance [4–6] and capacitance [7–9]. The electrical resistance can 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 deepgroove ball bearing using electrical capacitance measurements. The authors used insitu 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 parallelplate 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 threedimensional 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 piezoviscouselastic behavior of the lubricant film. The electrical capacitance of the lubricant is then introduced by means of the parallelplate 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.
2 Digitalt win of radial ball bearing
A threedimensional 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 ballsof 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.
Fig. 1 (a) Digital twin of radial ball bearing 6208, (b) SKF 6208 radial ball bearing. 
Geometric characteristics of 6208 radial ball bearing.
Fig. 2 Exploded view of digital twin (6208 radial ball bearing). 
3 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 KelvinVoigt SpringDashpot 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).
Fig. 3 Contact stiffness model. 
3.1 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). $$\text{{}{F}_{n,t}={K}_{n,t}\text{\hspace{0.17em}}{\delta}_{n,t}+{C}_{n,t}\text{\hspace{0.17em}}{v}_{n,t}\text{\hspace{0.17em}with\hspace{0.17em}}{F}_{t}=\text{min}\left({F}_{t},\text{\hspace{0.17em}}\mu {F}_{n}\right)\times sgn\left({v}_{t}\right)$$(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 $\frac{{C}_{t}}{{C}_{n}}$ 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 nonlinear model; ${F}_{n}\propto {\delta}_{n}^{\frac{3}{2}}$. For a given rolling element in contact with the innerraceway (or outerraceway), 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]: $$\{\begin{array}{l}{K}_{n}=\frac{2\pi \kappa G}{1\nu}\sqrt{\frac{2\overline{E}{R}_{\text{curve}}}{9\overline{F}{}^{3}}}{\delta}_{n}^{\frac{1}{2}}\hfill \\ {R}_{\text{curve}}=\frac{{R}_{x}{R}_{y}}{{R}_{x}+{R}_{y}}\hfill \\ {R}_{x}=\frac{{R}_{1x}{R}_{2x}}{{R}_{2x}\pm {R}_{1x}},{R}_{y}=\frac{{R}_{1y}{R}_{2y}}{{R}_{2y}{R}_{1y}}\hfill \end{array}$$(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 $\left(\kappa ={\alpha}^{\frac{2}{\pi}},\alpha =\frac{{R}_{y}}{{R}_{x}}\right)$, $\overline{F}$ the approximate elliptic integral of first kind $\left(\overline{F}=1+q\text{ln}\alpha ,q=\frac{\pi}{2}1\right)$ and $\overline{E}$ the approximate elliptic integral of second kind $\left(\overline{E}=1+\frac{q}{\alpha}\right)$. The equivalent radius R_{x} is calculatedby taking into account the cases of convex contact (rolling element/innerraceway) and concave contact (rolling element/outerraceway), as given in equation (2). For the sake of clarity, we note that R_{b} =R_{1x} = R_{1y}, R_{2x} = R_{i} and R_{2y} = R_{r} (Fig. 4). $G=\frac{E}{2\left(1+\nu \right)}$ 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]: $${K}_{n}=\frac{4{E}_{\text{eq}}\sqrt{{R}_{\text{eq}}}}{3\left(1\nu \right)}{\delta}_{n}^{\frac{1}{2}}$$(3)
with E_{eq} the equivalent Young’s modulus $\left(\frac{1}{{E}_{\text{eq}}}=\frac{1{\nu}^{2}}{E}+\frac{1{\nu}_{c}^{2}}{{E}_{c}}\right)$, E_{c} the Young’s modulus (E_{c} = 3.3 GPa) and ν_{c} the Poisson ratio (ν_{c} = 0.41) of cage balls made with polyamide 66 (PA66), ${R}_{\text{eq}}=\frac{{r}_{i}{r}_{j}}{{r}_{i}+{r}_{j}}$ 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 (Figure 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 $\stackrel{\u20d7}{n}$. $${C}_{\text{hyst}}=\frac{3{\alpha}_{e}{K}_{n}{\delta}_{n}}{2}$$(4)
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. $${C}_{\text{fluid}}=\frac{3\pi \eta a}{\sqrt{2}}{\left(\frac{{R}_{b}}{{h}_{\text{min}}}\right)}^{\frac{3}{2}}$$(5)
where η is the dynamic viscosity of the fluid at operating temperature and atmospheric pressure [24], a the semimajor axis of the Hertz elliptical contact [25], h_{min} the minimum lubricant film thickness of the ballinnerraceway (or outerraceway) 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 velocityVerlet scheme [17].
Fig. 4 Curvatures in contact in two orthogonal cross sections of a ball bearing: (a) x−z plane, (b) y−z plane. 
3.2 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 piezoviscous 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 piezoviscouselastic 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 pointcontact EHD, the expression of the minimum fluid film thickness (6) is as follows: $$\begin{array}{lll}\frac{{h}_{\text{min}}}{{R}_{x}}\hfill & =\hfill & 1.637{\overline{U}}_{r}^{0,711{\kappa}^{0,023}}{\overline{E}}^{0.65{\kappa}^{0.045}}{\overline{W}}^{0.09{\kappa}^{0.15}}\hfill \\ \hfill & \hfill & \times \left(10.974{e}^{0.676\kappa}\right)\hfill \end{array}$$(6)
where $\overline{E},{\overline{U}}_{r},\overline{W}$ are three dimensionless parameters: $\overline{E}=\xi {E}^{\prime}$, with ξ the pressureviscosity coefficient of the lubricant at operating temperature and atmospheric pressure, ${E}^{\prime}=\frac{E}{1{\nu}^{2}}$ the effective elastic modulus of steel; ${\overline{U}}_{r}=\frac{\eta {U}_{r}}{{E}^{\prime}{R}_{x}}$, with U_{r} the entrainment speed; $\overline{W}=\frac{{Q}_{\psi}}{{E}^{\prime}{R}_{x}^{2}}$, with Q_{ψ} the ballinnerraceway (or outerraceway) 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]: $${\Lambda}_{f}=\frac{{h}_{\text{min}}({Q}_{\psi},\omega )}{\sqrt{{\sigma}_{1}^{2}+{\sigma}_{2}^{2}}}$$(7)
where σ_{1} and σ_{2} denote the standard mean roughnesses of the measured surfaces for the innerraceway (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 $\left({\sigma}_{2}\ll {\sigma}_{1}\right)$. In addition, the schematic representation of the Stribeck curve (Fig. 5) for nonconformal 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, mixedfilm or fullyfluid film) in relation to fluid parameter Λ_{f}.
Fig. 5 Schematic representation of Stribeck curve: friction coefficient and fluid film thickness as function of fluid parameter. 
3.3 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 parallelplate formula the lubricant film thickness is inversely proportionalto the capacitance (8). $${C}_{\text{Hertz}}={\epsilon}_{0}{\epsilon}_{r}\frac{{A}_{\text{Hertz}}}{{h}_{c}}$$(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 $\left({\epsilon}_{0}=8.854\text{pF}.{\text{m}}^{1}\right)$. 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 semiaxes 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): $$\begin{array}{lll}\frac{{h}_{c}}{{R}_{x}}\hfill & =\hfill & 3.672{\overline{U}}_{r}^{0,663{\kappa}^{0,025}}{\overline{E}}^{0.502{\kappa}^{0.064}}{\overline{W}}^{0.045{\kappa}^{0.18}}\hfill \\ \hfill & \hfill & \times \left(10.573{e}^{0.74\kappa}\right).\hfill \end{array}$$(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 ClausiusMossotti model [32]: $$\frac{{\epsilon}_{r}1}{{\epsilon}_{r}+2}=\frac{{N}_{A}\beta \rho}{3M}$$(10)
where N_{A} is the Avogadro number, β the molecular polarizability, M the molecular weight and ρ the density. Thefactor $\frac{{N}_{A}\beta}{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): $$\rho ={\rho}_{0}\frac{0.59\cdot {10}^{9}+1.34\cdot {P}_{\text{Hertz}}}{0.59\cdot {10}^{9}+{P}_{\text{Hertz}}}$$(11)
with ρ_{0} the ambientpressure density and P_{Hertz} the Hertz contact pressure. In Figure 6a are plotted the dielectric constant evolution of synthetic oil film (Tab. Table 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 Fr = 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 ofcapacitances 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. $${C}_{\text{bearing}}={\sum}_{i=1}^{{N}_{l}}\left(\frac{{C}_{\text{Hertz}}^{\text{inner},i}\times {C}_{\text{Hertz}}^{\text{outer},i}}{{C}_{\text{Hertz}}^{\text{inner},i}+{C}_{\text{Hertz}}^{\text{outer},i}}\right)$$(12)
with N_{l} the number of rolling elements in the loaded zone of the ball bearing. The contact capacitances ${C}_{\text{Hertz}}^{\text{inner},i}$ and ${C}_{\text{Hertz}}^{\text{outer},i}$ between the ith rolling element and innerouter raceways, respectively, are determined according to the parallelplate 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].
Fig. 6 (a) Dielectric constant as function of Hertz contact pressure (6208 bearing – PAO VG48 at 40 °C). (b) Hertz contact pressure as function of azimuth angle ψ. 
Properties of the lubricants.
Fig. 7 Flowchart of equivalent circuit with parallel capacitance of ball bearing. 
4 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.
4.1 Effects of radial load and diametral clearance on electrical capacitance
According to the parallelplate 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 tothe 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 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 fluidfilm 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 $\frac{{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 $\frac{{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 $\frac{{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).
Fig. 8 (a) Capacitance distribution in the loaded zone with $\frac{{P}_{d}}{D}=0$ and U_{r} = 0.41 m/s. (b) Capacitance at the maximum radial load position $\left(\psi ={0}^{\xb0}\right)$ as function of F_{r}. 
Simulated clearances P_{d} and corresponding load parameter ε.
Fig. 9 Capacitance distribution in the loaded zone with F_{r} = 5000 N and U_{r} = 0.41 m/s. 
4.2 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 $\left(\epsilon =0.5\right)$. 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 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 conductivityin 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 bearingas 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 xaxis 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.
Fig. 10 Capacitance as function of entrainment speed U_{r}: (a) mineral oil, (b) synthetic oil PAO VG48. 
Fig. 11 Capacitancedielectric constant ratio as function of fluid film thickness of the elementinner raceway contact. 
Fig. 12 (a) Film parameter Λ_{f} as function of entrainment speed U_{r}. (b) Correlation of the fluid parameter Λ_{f} and ball bearing capacitance. 
5 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 parallelplate 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.
Nomenclature
α_{e}: Restitution coefficient of steel (s.m^{−1})
δ_{n,t}: Normal/tangential relative displacement (μm)
η: Dynamic viscosity of the fluid (m.Pa.s)
ω: Shaft/innerring angular speed (rad.s^{−1})
$\overline{E}$: Approximate elliptic integral of second kind
$\overline{F}$: Approximate elliptic integral of first kind
$\overline{E}$: Dimensionneless material parameter
${\overline{U}}_{r}$: Dimensionneless entrainment speed
${\overline{W}}_{r}$: Dimensionneless load parameter
ρ_{0}: Ambientpressure density of the lubricant
σ_{1,2}: Arithmetic mean roughnesses (μm)
ε_{0}: Vacuum permitivity (pF.m^{−1})
ξ: Viscositypressure coefficient of the lubricant (GPa^{−1})
a: Semimajor axis of the Hertz’s elliptical contact (mm)
b: Semiminor axis of the Hertz’s elliptical contact (μm)
C_{n,t}: Normal/tangential viscous damping coefficient (N.s.m^{−1})
D: Rolling element diameter (mm)
E: Young’s modulus of steel (GPa)
E′: Effective elastic modulus of steel (GPa)
G: Shear modulus of steel (GPa)
h_{c}: Central lubricant film thickness (nm)
h_{min}: Minimum film fluid thickness (nm)
K_{n,t}: Contact normal/tangential stiffness (N.m^{−1})
P_{d}: Diametral clearance (μm)
Q_{max}: Maximum normal load at ψ = 0 (N)
R_{b}: Rolling element radius (mm)
${R}_{c}^{i}$: Curvature radius of inner raceway (mm)
${R}_{c}^{o}$: Curvature radius of outer raceway (mm)
R_{curve}: Curvature sum radius (mm)
U_{r}: Entrainment speed (m.s^{−1})
v_{t}: Sliding velocity (m.s^{−1})
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, $\epsilon =\frac{1}{2}\left(1\frac{{P}_{d}}{2{\delta}_{r}}\right)$, 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]: $${Q}_{\psi}={Q}_{\text{max}}{\left[1\frac{(1\text{cos}\psi )}{2\epsilon}\right]}^{3/2}$$(A.1)
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\left({Q}_{\text{max}}=\frac{{F}_{r}}{Z{J}_{r}\left(\epsilon \right)}\right)$.
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 $\left(\frac{{P}_{d}}{D}=0.0045\right)$, zero clearance $\left(\frac{{P}_{d}}{D}=0\right)$ and positive clearance $\left(\frac{{P}_{d}}{D}=0.0045\right)$, respectively. The azimuth angle ψ of the rolling element is varied by driving the inner ring at a slow angular speed $\left(\omega =1\text{rad}\times {\text{s}}^{1}\right)$ 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].
Fig. A.1 Radial shift δ_{r} due to the imposed radial load F_{r}. 
Fig. A.2 Variation of $\frac{{Q}_{\psi}}{{F}_{r}}$ ratio as function of azimuth angle ψ under radial load F_{r} = 5000 N for three values of load parameter ε. 
Fig. A.3 Schematic representation of the axial free play. 
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: $$C=\underset{B}{\underbrace{\left(\frac{{R}_{c}^{i}+{R}_{c}^{o}}{D}1\right)}}\times D$$(A.2)
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=\left(\frac{2\times {R}_{c}^{i,o}}{D}1\right)$. From Figure A.3, one can express the free contact angle α^{0} made by thevertical line and the line passing through the points of contact between the rolling element and both raceways as follows: $${\alpha}^{0}={\text{cos}}^{1}\left(1\frac{{P}_{d}}{2C}\right).$$(A.3)
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: $${P}_{e}=2C\times \text{sin}{\alpha}^{0}.$$(A.4)
Numerical simulations with digital twin of ball bearing are conducted with several values of parameter B summarized in Table A.1. Raceway curvatures ${R}_{c}^{i,o}$ and conformities $\frac{{R}_{c}^{i,o}}{D}$ are deduced from the expression of the total curvature B.
Values of the total curvature B with their associated raceways curvatures ${R}_{c}^{i,o}$ and conformities $\frac{{R}_{c}^{i,o}}{D}$.
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.
Fig. A.4 Free contact angle and free endplay for several values of total curvature B. 
References
 A.D. Roberts, D. Tabor, Fluid film lubrication of rubber – an interferometric study, Wear, 11, 163–166 (1968) [Google Scholar]
 W.M. Nozhat, Measurement of liquidfilm thickness by laser interferometry, Appl. Opt. 36, 7864–7869 (1997) [CrossRef] [PubMed] [Google Scholar]
 S.Y. Jang, Computation of fluid film pressures by measuring the elastohydrodynamic lubrication film thickness with nanoscale resolution, In Experimental Mechanics in Nano and Biotechnology, volume 326 of Key Engineering Materials, pages 413–416. Trans Tech Publications, Switzerland, 2006, p. 12 [Google Scholar]
 J.F. Archard, M.T. Kirk, Lubrication at point contacts, Proc. R. Soc. London A 261, 532–550 (1962) [Google Scholar]
 T.E. Tallian, Y.P. Chiu, D.F. Huttenlocher, J.A. Kamenshine, L.B. Sibley, N.E. Sindlinger, Lubricant films in rolling contact of rough surfaces, ASLE Trans., 7, 109–126 (1964) [CrossRef] [Google Scholar]
 H.A. Spikes, Mixed lubrication – an overview, Lubrication Sci. 9, 221–253 (1997) [CrossRef] [Google Scholar]
 M.R. Özgü, J.C. Chen, N. Eberhardt, A capacitance method for measurement of film thickness in twophase flow, Rev. Sci. Instrum. 44, 1714–1716 (1973) [Google Scholar]
 J. Klausner, L.Z. Zeng, D.M. Bernhard, Development of a film thickness probe using capacitance for asymmetrical twophase flow with heat addition, Rev. Sci. Instrum. 63, 3147–3152 (1992) [Google Scholar]
 G.E. Thorncroft, J. Klausner, A capacitance sensor for twophase liquid film thickness measurements in a square duct, J. Fluids Eng. Trans. ASME 119, 03 (1997) [CrossRef] [Google Scholar]
 S. Murer, F. Bogard, L. Rasolofondraibe, B. Pottier, P. Marconnet, Determination of loads transmitted by rolling elements in a roller bearing using capacitive probes: Finite element validation, Mech. Syst. Signal Process. 5455, 02 (2015) [Google Scholar]
 C. Machado, M. Guessasma, V. Bourny, Electromechanical prediction of the regime of lubrication in ball bearings using Discrete Element Method, Tribol. Int. 127, 69–83 (2018) [Google Scholar]
 K.A. Jablonka, R. Glovnea, J. Bongaerts, Quantitative measurements of film thickness in a radially loaded deepgroove ball bearing, Tribol. Int. 119, 239–249 (2018) [Google Scholar]
 M.Y. Toumi, S. Murer, F. Bogard, F. Bolaers, Numerical simulation and experimental comparison of flaw evolution on a bearing raceway: Case of thrust ball bearing, J. Comput. Des. Eng. 5, 427–434 (2018) [Google Scholar]
 P.K. Gupta, Advanced dynamics of rolling elements, Springer, New York, 1984 [CrossRef] [Google Scholar]
 L. Xu, Y. Yang, Y. Li, C. Li, S. Wang, Modeling and analysis of planar multibody systems containing deep groove ball bearing with clearance, Mech. Mach. Theory 56, 69–88 (2012) [Google Scholar]
 C. Machado, M. Guessasma, E. Bellenger, An improved 2D modeling of bearing based on DEM for predicting mechanical stresses in dynamic, Mech. Mach. Theory, 113, 53–66 (2017) [Google Scholar]
 M. Guessasma, C. Machado, Threedimensional DEM modelling of ball bearing with lubrication regime prediction, Lubricants 6, 2 (2018) [Google Scholar]
 K. Bourbatache, M. Guessasma, E. Bellenger, V. Bourny, J. Fortin, DEM ball bearing model and defect diagnosis by electrical measurement, Mech. Syst. Signal Process. 41, 98–112 (2013) [Google Scholar]
 P.A. Cundall, O.D.L. Strack, A discrete numerical model for granular assemblies, Geotechnique 29, 47–65 (1979) [Google Scholar]
 B.J. Hamrock, W.J. Anderson, RollingElement Bearings. For sale by the National Technical Information Service, 1983. [Google Scholar]
 R.D. Mindlin, H. Deresiewicz, Elastic spheres in contact under varying oblique force, ASME J. Appl. Mech. 20, 327–344 (1953) [Google Scholar]
 K.H. Hunt, F.R.E. Crossley, Coefficient of restitution interpreted as damping in vibroimpact, J. Appl. Mech. 42, 440–445 (1975) [Google Scholar]
 T.L.H Walford, B.J. Stone, The sources of damping in rolling element bearings under oscillating conditions, The Institution of Mechanical Engineers, London, 1983, 197 [Google Scholar]
 A. Harnoy, Bearing Design in Machinery: Engineering Tribology and Lubrication. Marcel Dekker, New York, 2002 [Google Scholar]
 B.J. Hamrock, D. Dowson, Ball Bearing Lubrication: Elastohydrodynamics of Elliptical Contacts, John Wiley & Sons Inc, Chichester, 1981 [Google Scholar]
 C. Mary, D. Philippon, N. Devaux, N. Fillot, D. Laurent, S. Bair, P. Vergne, Bridging high pressure rheology and filmforming capacity of polymerbase oil solutions in EHL, Tribol. Int. 93, 502–510 (2016) [Google Scholar]
 M. Masjedi, M.M. Khonsari, On the effect of surface roughness in pointcontact EHL: Formulas for film thickness and asperity load, Tribol. Int. 82, 228–244 (2015) [Google Scholar]
 T.E. Tallian, Paper 14: Rolling contact failure control through lubrication. Proc. Inst. Mech. Eng. 182, 205–236 (1967) [Google Scholar]
 R. Stribeck, Ball bearings for various loads, Trans. ASME 29, 420–463 (1907) [Google Scholar]
 A.R. Wilson, The relative thickness of grease and oil films in rolling bearings, Proc. Inst. Mech. Eng. 193, 185–192 (1979) [CrossRef] [Google Scholar]
 P. Brüser, Untersuchungen uber die elastohydrodynamische schmierfilmdicke, 1972 [Google Scholar]
 J.D. Jackson, Classical Electrodynamics, 2nd edn., John Wiley & Sons, New York, 1999 [Google Scholar]
 D. Dowson, G.R. Higginson, Elastohydrodynamic Lubrication: the Fundamentals of Rollers and Gear Lubrication, Pergamon, Oxford, 1966 [Google Scholar]
 K.A. Jablonka, R. Glovnea, J. Bongaerts, Evaluation of ehd films by electrical capacitance, J. Phys. D 45, 85301–385309 (2012) [CrossRef] [Google Scholar]
 T.J. Harvey, R.J.K. Wood, H.E.G. Powrie, C. Warrens, Charging ability of pure hydrocarbons and lubricating oils, Tribol. Trans. 47, 263–271 (2004) [CrossRef] [Google Scholar]
 P. Brüser, Untersuchungen über die elastohydrodynamische schmierfilmdicke, 1972 [Google Scholar]
 T.A. Harris, M.N. Kotzalas, Advanced Concepts of Bearing Technology, In: Rolling Bearing Analysis, Fifth Edtion Series, Taylor & Francis, Abingdon, 2006 [Google Scholar]
Cite this article as: S.O. Farah, M. Guessasma, E. Bellenger, Digital twin by DEM for ball bearing operating under EHD conditions, Mechanics & Industry 21, 506 (2020)
All Tables
Values of the total curvature B with their associated raceways curvatures ${R}_{c}^{i,o}$ and conformities $\frac{{R}_{c}^{i,o}}{D}$.
All Figures
Fig. 1 (a) Digital twin of radial ball bearing 6208, (b) SKF 6208 radial ball bearing. 

In the text 
Fig. 2 Exploded view of digital twin (6208 radial ball bearing). 

In the text 
Fig. 3 Contact stiffness model. 

In the text 
Fig. 4 Curvatures in contact in two orthogonal cross sections of a ball bearing: (a) x−z plane, (b) y−z plane. 

In the text 
Fig. 5 Schematic representation of Stribeck curve: friction coefficient and fluid film thickness as function of fluid parameter. 

In the text 
Fig. 6 (a) Dielectric constant as function of Hertz contact pressure (6208 bearing – PAO VG48 at 40 °C). (b) Hertz contact pressure as function of azimuth angle ψ. 

In the text 
Fig. 7 Flowchart of equivalent circuit with parallel capacitance of ball bearing. 

In the text 
Fig. 8 (a) Capacitance distribution in the loaded zone with $\frac{{P}_{d}}{D}=0$ and U_{r} = 0.41 m/s. (b) Capacitance at the maximum radial load position $\left(\psi ={0}^{\xb0}\right)$ as function of F_{r}. 

In the text 
Fig. 9 Capacitance distribution in the loaded zone with F_{r} = 5000 N and U_{r} = 0.41 m/s. 

In the text 
Fig. 10 Capacitance as function of entrainment speed U_{r}: (a) mineral oil, (b) synthetic oil PAO VG48. 

In the text 
Fig. 11 Capacitancedielectric constant ratio as function of fluid film thickness of the elementinner raceway contact. 

In the text 
Fig. 12 (a) Film parameter Λ_{f} as function of entrainment speed U_{r}. (b) Correlation of the fluid parameter Λ_{f} and ball bearing capacitance. 

In the text 
Fig. A.1 Radial shift δ_{r} due to the imposed radial load F_{r}. 

In the text 
Fig. A.2 Variation of $\frac{{Q}_{\psi}}{{F}_{r}}$ ratio as function of azimuth angle ψ under radial load F_{r} = 5000 N for three values of load parameter ε. 

In the text 
Fig. A.3 Schematic representation of the axial free play. 

In the text 
Fig. A.4 Free contact angle and free endplay for several values of total curvature B. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.