Rolling of an Elastic Ellipsoid Upon an Elastic-Plastic Flat

The paper presents a numerical analysis of the rolling contact between an elastic ellipsoid and an elastic-plastic ﬂat. Numerical simulations have been performed with the help of a contact solver called Plast-Kid®, with an algorithm based on an integral formula-tion or semi-analytical method. The application of both the conjugate gradient method and the discrete convolution and fast Fourier transform technique allows keeping the computing time reasonable when performing transient 3D simulations while solving the contact problem and calculating the subsurface stress and strain states. The effects of the ellipticity ratio k—ranging from 1 to 16—and of the normal load—from 4.2 GPa to 8 GPa —are investigated. The reference simulation corresponds to the rolling of a ceramic ball on a steel plate made of an AISI 52100 bearing steel under a load of 5.7 GPa. The results that are presented are, ﬁrst, the permanent deformation of the surface and, sec-ond, the contact pressure distribution, the von Mises stress ﬁeld, the hydrostatic pressure, and the equivalent plastic strain state within the elastic-plastic body. A comparison with an experimental surface deformation proﬁle is also given to validate the theoretical background and the numerical procedure.


Introduction
Many structural materials exhibit a strain limit under load beyond which full recovery of the initial geometry is not possible when the load is removed. A bearing steel loaded in compression behaves in a similar manner. Thus, when a ball is normally loaded on a bearing raceway, an indentation may remain in the raceway and the ball may exhibit a flat spot after unloading if the yield stress is exceeded ͓1͔. These permanent deformations may be at the origin of unexpected vibrations and/or entail the bearing life by the accumulation of plastic strain or by stress concentration due to the localized change of the surface geometry.
The basic static load rating of a rolling bearing is defined as that load applied to a nonrotating bearing that will result in permanent deformation of 0.0001ϫ D at the weaker of the inner or outer raceway contact occurring at the position of the maximum loaded rolling element, D being the ball diameter. Empirical formulas are given by Harris ͓1͔. This paper presents a numerical analysis of the permanent print produced by the rolling of an elastic body upon an elastic-plastic ͑EP͒ flat when plastic flow occurs.
The EP contact solver is based on a semianalytical method ͑SAM͒. The original concept and first attempt to solve 3D problems by a SAM are due to Jacq et al. ͓2͔. Compared to the finite element method ͑FEM͒, the SAM shows much shorter computation times, typically by several orders of magnitude. In this method, analytical formulas are derived using Green's functions, commonly called influence coefficients in the discrete form. Quantities are then obtained by numerical computing using accelerating techniques, leading to extremely short computation times. Among many numerical methods one of the most efficient procedures consists of coupling the conjugate gradient method ͑CGM͒ to solve the contact problem ͑Polonsky and Keer ͓3͔͒ and the discrete convolution fast Fourier transform ͑DC-FFT͒ introduced by Liu et al. ͓4͔ in computing the integrals.
The SAM has already been successfully applied in solving EP contacts, vertically loaded or subjected to a rolling load ͓2͔, with or without considering the frictional heating at the contact interface ͓5͔, considering also the effect of a Coulomb friction coefficient ͓6͔, and more recently in a displacement driven formulation for the tugging between asperities ͓7͔. The method has been successfully used to determine the microyield stress profile of a nitrided steel from nanoindentation experiments ͓8͔, to investigate the rolling of a load on a smooth, dented, or rough surface ͓9,10͔, and to simulate fretting wear ͓11,12͔ and the running-in or wear of initially smooth or rough surfaces ͓13͔. The EP contact has been also studied with a very similar method by Wang and Keer ͓14͔, who investigated the effect of the hardening model, and by Popescu et al. ͓15,16͔, who considered an initial stress state. The way to consider the rolling motion of the contacting bodies consists of solving the EP contact at each time step while upgrading the geometries as well as the hardening state along the moving direction. This paper investigates the effect of an overload on the permanent deformation of the surface and subsequent subsurface stress and strain states. The effect of the ellipticity ratio k = c / a is first studied, with values ranging from 1 to 16 and a semiminor axis a along the rolling direction. In a second step and for an ellipticity ratio of 8, the effect of the normal load-from 4.2 GPa to 8 GPa-is also investigated. Results are of prime importance for rolling bearing manufacturers and users.
2 Elastic-Plastic Contact Model 2.1 Theory. The current model is based on the work of Liu et al. ͓4͔ and of Nélias and co-workers ͓5-13͔, who developed an EP contact solver using an integral formulation with intermediate analytical solutions, also called SAM. The three-dimensional EP contact model was initially built to study the vertical or rolling loading of a smooth body against a dented ͓2͔ or rough surface ͓10͔. This model includes two important modules, the contact solver that determines the contact area and the pressure distribution, and the plasticity loop that is used to calculate the plastic region ͑where the plastic strain is not nil͒. Note that the presence of a plastic region will modify the geometry of the contacting surfaces; therefore, it is required to iterate between the contact solver and the plasticity loop. The CGM is used in the contact solver to accelerate the computation. Convolution products are

Vincent Boucly
LaMCoS, INSA-Lyon, CNRS UMR5259, Villeurbanne F69621, France performed using the DC-FFT, following the method proposed by Liu et al. in Ref. ͓3͔. In order to simulate the rolling/sliding contact, a load-driven formulation was first used by applying a normal load ͑vertical loading͒ prior to the tangential displacement of the load ͑rolling load͒.
The set of equations to be solved simultaneously for the contact problem is recalled below: u ij being the composite displacement of both bodies at point ͑i , j͒ of the contact area I c , h ij the surface separation, ␣ the rigid body approach, p ij the pressure, a x and a y the grid spacing in x and y directions, respectively, and W the total load. In addition to the main hypothesis of the contact mechanics, which is that the contact area is small in comparison to the dimensions of bodies justifying the assumption of half-spaces, displacements and strains are small, alloying the superimposition of inelastic and elastic contributions. Plasticity effects are introduced through Betti's reciprocal theorem, leading to a formulation where the surface displacements u and the subsurface stresses ij are expressed as the summation of the contribution of the elastic field, the residual state, and the thermally induced strains ͑see Eqs. ͑2͒ and ͑3͒͒. u͑A͒ = u e ͑A͒ + u r ͑A͒ + u t ͑A͒ A being a point on the surface ͑2͒ ij ͑B͒ = ij e ͑B͒ + ij r ͑B͒ + ij t ͑B͒ B being a point in the volume

͑3͒
Plasticity is an irreversible phenomenon that requires an incremental description. In a general incremental formulation of plasticity, a plastic strain increment depends on the stress, the stress increment, and on the hardening parameters. A return mapping algorithm with an elastic predictor/plastic corrector scheme was implemented in the plasticity loop for the integration of the EP constitutive equations ͑see Fotiu and Nemat-Nasser ͓17͔͒. The geometry of the contacting bodies can also evolve during the loading history due to the occurrence of a permanent deformation of the surface when the yield stress is exceeded, which will indeed modify the contact pressure distribution. Therefore, the relation between plastic strain and contact pressure must also be incremental. A two-stage incremental formulation of the EP contact is then used.
To complete this incremental formulation, the loading history must be defined. Two types of load increments are considered. The first one is a vertical loading or unloading without a rolling movement. The only change is an increase or a decrease of the applied external load W. The second type of load increment corresponds to a rolling movement of the load. The applied external load does not change. Considering a mark attached to the contact, the plastic strain, hardening state, and contact pressure must be shifted from one iteration to the next one. The contact radius for the reference simulation is a = 310 m, and the mesh grid is here equal to 25 m in both rolling and transverse directions.

Reference Simulation and
An isotropic hardening model along with the von Mises criterion has been chosen for the representation of the yield surface. The hardening law of the EP body is described by Swift's law ͑Eq. ͑4͒͒, which presents the numerical advantage of being continuously derivable. Parameters describing the hardening of the AISI 52100 bearing steel are B = 945 MPa, C = 20, and n = 0.121, with p as the equivalent plastic strain, eq = B͑C + 10 6 ϫ p ͒ n ͑4͒ The first step consists in determining the subsurface stress and strain states found under load or after unloading, as well as the shape of the permanent print found on the surface of the EP flat. Figure 1 shows the longitudinal profile ͑permanent͒ found after one loading cycle and after unloading with and without the rolling of the ball under a contact pressure of 5.7 GPa ͑Hertz pressure, i.e., elastic͒. The curve with square symbols corresponds to the axisymmetrical print found after unloading and without rolling. The curve with small dots presents the permanent print after a vertical loading, rolling of 4 mm from left to right, maintaining the normal load constant, and unloading ͑vertical͒. Results are presented in a dimensionless form with the abscissa and ordinate normalized by the contact radius a.
It should be noted that the maximum depth is higher for the rolling load compared to the purely vertical loading/unloading. One may observe a decrease of the print depth after a transient period when the load is rolling, as previously observed by Jacq et al. ͓2͔. This depth reaches an asymptotic value when rolling a few times the contact radius, which corresponds to a steady-state regime. Finally, a plastically deformed shoulder higher at the unloading position than at the initially loaded area can be observed.
A comparison of the surface displacement found at the center of the rolling track after unloading for the three first cycles and a rolling distance of 2 mm is shown in Fig. 2 for the longitudinal profile and in Fig. 3 for the transverse profile in the steady-state region. It can be observed that the permanent print reaches quickly a stabilized value after very few cycles, i.e., two or three rolling loading cycles, which is due to the isotropic hardening model used in this study. The evolution between the first and second cycles is explained by the change of the surface conformity due to the plastic print that occurs mostly after one loading cycle.
An experimental profile in the stationary region is also given in Fig. 3, and a very good agreement is found. Note that starting from Fig. 3, all results presented in the transversal direction will Fig. 1 Comparison of a vertical loading/unloading with and without rolling "k =1−P Hertz = 5.7 GPa, ceramic ball "elastic…/flat 52100 "elastic-plastic… surface, after a single cycle, rolling distance of 4 mm… be presented in the stationary region. Figures 4 and 5 show a cross-section view of the von Mises stress field under load and after unloading, respectively, in a plane perpendicular to the rolling motion and in the steady-state region.
The stress values, which are normalized by the elastic contact pressure P Hertz , are given after two loading cycles. The loaded case ͑Fig. 4͒ exhibits a maximum stress that is lower than the elastic solution due to plasticity. More interesting is the distribution of the residual stress ͑see Fig. 5͒. Three areas with high residual stress coexist, in fact, three bands along the rolling direction, one at the Hertzian depth and two others at the surface near the border of the contact but slightly inside the permanent print plotted in Fig. 3. The residual stress distribution in the plane of symmetry ͑y =0͒ is also given in Fig. 6.
It is not possible to know whether a region is in a tensile or a compressive state when using the von Mises equivalent stress. It is why the hydrostatic stress, also called the mean stress, is now used to describe the residual stress. A positive value indicates that a region is in a tensile state, while a negative value corresponds to compression. It becomes clear from Fig. 7 that the near surface area is in a tensile state after unloading, while subsurface plasticity induces compressive stresses at the Hertzian depth. This result is of prime interest for fatigue applications.
A cross-section view of the residual hydrostatic pressure is presented in Fig. 7. Conversely to Fig. 5, which corresponds to the von Mises stress, a tensile zone lying on the surface is now observed, with two local maxima on both sides of the rolling track.
This peculiar feature is of prime interest for a rolling contact fatigue since these tensile areas may favor the initiation and propagation of surface or near surface cracks due to surface inhomogeneities-such as embedded inclusions-or geometrical defects-such as grinding furrows or debris denting.
The distribution of the equivalent plastic strain in the plane y = 0 is given in Fig. 8 for a rolling of 4 mm from the left to the right and after two loading cycles. One may notice a maximum value of approximately 0.9% at the Hertzian depth.

Effect of the Ellipticity Ratio
Some trends when the ellipticity ratio k increases up to 16 are now presented to give more realistic results for rolling bearing applications. As in ball bearings the elliptical contact is elongated in the direction perpendicular to the rolling movement. In all simulations, the grid size is kept equal to 25 m along the rolling direction but is multiplied by the ellipticity ratio in the transversal direction ͑i.e., ⌬y = k ϫ 25 m͒. The distance of rolling is 4 mm, Fig. 2 Surface displacement at the center of the rolling track after unloading "k =1−P Hertz = 5.7 GPa, ceramic ball "elastic…/flat 52100 "elastic-plastic… surface, after one, two, and three cycles, rolling distance of 2 mm… Fig. 3 Transverse surface profile after unloading "k =1−P Hertz = 5.7 GPa, ceramic ball "elastic…/flat 52100 "elastic-plastic… surface, after one, two, and three cycles… i.e., 12.9 times the semiminor axis of the contact, which is required to reach the steady-state regime for high ellipticity ratios.
As in the previous section, all results will be given in a dimensionless form, except for the equivalent plastic strain which is given in percent, i.e. the coordinates are normalized by the contact radius of the reference test ͑a = 310 m͒, and the stresses by the Hertzian pressure ͑5.7 GPa͒. Again the cross-section views correspond to a transversal cut at an abscissa in the steady-state region, i.e., far from the loading and unloading points.
The geometry of the elastic body leading to an elliptical contact of various ellipticity ratios is summarized in Table 1. It could be observed that the radius ratio R y / R x = 77.88 produces a semimajor axis of the contact ellipse 16 times greater than the semiminor axis. In the meantime, the normal load should be multiplied by the ellipticity ratio to keep the Hertzian contact pressure constant. Note that the contact dimensions and normal load indicated in Table 1 correspond to a contact pressure of 5.7 GPa as an elastic solution ͑Hertz͒, which indeed differs from the EP solution.
Figures 9 and 10 present a longitudinal and a transversal view of the permanent print found after the first and the second cycles for different ellipticity ratios. It is interesting to note that despite the same maximum contact pressure of 5.7 GPa, the print depth decreases when the ellipticity ratio increases ͑Fig. 9͒ whereas its width increases ͑Fig. 10͒. Another interesting feature is the formation of a material pileup in front of the unloading point and on each side of the rolling track. For k = 1, the maximum height is found on each side of the rolling track, the shoulder ahead of the track at the unloading position ͑Fig. 9͒ being slightly lower than on the lateral sides. Conversely, at higher ellipticity ratios, one may observe that the height of the shoulder found in front of the rolling track increases significantly ͑Fig. 10͒, whereas it decreases besides the rolling track. It should be also noted that the longitudinal print found for k = 16 is close to the line contact solution.
The effect of the ellipticity ratio on the steady-state contact pressure distribution is given in Fig. 11. A comparison with the elastic solution shows a more pronounced decrease of the maximum EP contact pressure when increasing the ellipticity ratio, until an asymptotic solution corresponding to the line contact solution. More surprising is the associated profile for the equivalent plastic strain ͑see Fig. 12͒, which indicates a decrease of the maximum of the equivalent plastic strain found at the Hertzian depth when the ellipticity ratio increases. The depth at which this maximum is located also increases with the k ratio, which is in agreement with the elastic theory, i.e., 0.48a for a circular point contact and 0.8a for a line contact. Finally, the flattening of the contact pressure distribution when increasing the ellipticity ratio could be attributed to the extension of the plastic zone in depth, as seen in Fig. 12, despite a lower plasticity level ͑denoted by the maximum value of the equivalent plastic strain͒. Figures 13-15 describe the stress profiles found below a point located at the center of the rolling track and in the steady-state region. The von Mises stress distribution is first given under load in Fig. 13, then after unloading in Fig. 14. One may observe in Fig. 6 Distribution of the residual stress in the plane y =0 "k =1−P Hertz = 5.7 GPa, ceramic ball "elastic…/flat 52100 "elastic-plastic… surface, after the second loading cycle… Fig. 7 Cross-section view of the hydrostatic pressure after unloading "residual stress… in the steady-state region "k =1 − P Hertz = 5.7 GPa, ceramic ball "elastic…/flat 52100 "elastic-plastic… surface, after the second loading cycle… Fig. 8 Isovalues of the equivalent plastic strain "in %… in the plane y =0 "k =1−P Hertz = 5.7 GPa, ceramic ball "elastic…/flat 52100 "elastic-plastic… surface, after the second loading cycle… Table 1 Effect of the ellipticity ratio k = c / a on the radius of the elastic body and, subsequently, on the normal load to keep the Hertzian contact pressure equal to 5.7 GPa  Fig. 13 a decrease with the ellipticity ratio of the maximum stress found under loading, along with an increase of the depth at which the maximum value is found, in agreement with the trends for the equivalent plastic strain ͑Fig. 12͒. The residual von Mises stress plotted in Fig. 14 shows a high level of residual stress, up to 10% of the Hertz pressure for k = 1, i.e., 570 MPa. It can also be observed that there is almost no difference between the stress state after the first and the second cycles. A convenient way to state whether a region is in a tensile or a compressive state is to plot the hydrostatic pressure, as in Fig. 15. It becomes clear that the stress state within the depth is a succession of tensile, compressive, and tensile regions, i.e., compressive at the Hertzian depth and tensile elsewhere, including at the surface. It should be noted that, coming back to Fig. 7, the highest level of tensile stress found at the surface is located on each side of the rolling track and not in the plane of symmetry.

Effect of the Normal Load
The effect of the maximum contact pressure-ranging from 4.2 GPa to 8 GPa-is now investigated. The corresponding normal load is given in Table 2. In the present simulations, the ellip- Fig. 9 Surface displacement at the center of the rolling track after unloading for various ellipticity ratios and after one and two cycles "P Hertz = 5.7 GPa, ceramic ellipsoid "elastic…/flat 52100 "elastic-plastic… surface…. Loading, rolling, and unloading are represented by the thick arrows Fig. 10 Transverse surface profile after unloading for various ellipticity ratios and after one and two cycles. Profiles taken in the steady-state region "P Hertz = 5.7 GPa, ceramic ellipsoid "elastic…/flat 52100 "elastic-plastic… surface… ticity ratio k has been chosen to be equal to 8, which is a close value to what is found in ball bearings. In this section, results are presented with engineering values, i.e., distances in micrometers and stress in mega pascals, since both the contact ellipse and the maximum pressure are varying.
All results presented in this section correspond to a profile along a vertical line from a surface point located at the center of the rolling track and in the steady-state region. The equivalent plastic strain profile is first given ͑see Fig. 16͒, then the residual von Mises stress ͑after unloading͒ Fig. 17, and finally the hydrostatic pressure found after unloading ͑see Fig. 18͒. It is clear from Fig. 16 that the maximum strain found at 4.2 GPa, i.e., 0.05% is far below the conventional yield strain defined at 0.2%. Higher contact pressures result in a significant level of plasticity corresponding to 0.4%, 0.8%, and 2.1% of plastic strain at 5.7 GPa, 6.5 GPa, and 8 GPa, respectively. The depth at which these maxima are found increases with the load level, as predicted by the conventional elastic theory. Note that the plastic volume increases significantly with the load. Here, the von Mises stress profile under loading is close to the elastic solution since the plastic strain remains lower than 2.1% at the highest ͑for 8 GPa͒. Fig. 11 Contact pressure distribution found in the steady-state region for various ellipticity ratios and after the second cycle "P Hertz = 5.7 GPa, ceramic ellipsoid "elastic…/flat 52100 "elastic-plastic… surface… Fig. 12 Profile of the equivalent plastic strain found at the center of the track and in the steady-state region for various ellipticity ratios and after one and two cycles "P Hertz = 5.7 GPa, ceramic ellipsoid "elastic…/flat 52100 "elastic-plastic… surface… More interesting is the level of residual stress, either the von Mises stress ͑Fig. 17͒ or the hydrostatic pressure ͑Fig. 18͒, which exhibits a high level of stress, i.e., 40 MPa, 253 MPa, 445 MPa, and 935 MPa for the equivalent von Mises stress at 4.2 GPa, 5.7 GPa, 6.5 GPa, and 8 GPa, respectively. Figure 18 indicates that the combination of a high ellipticity ratio ͑k =8͒ with a high contact pressure ͑6.5 GPa or higher͒ results in the disappearance of the tensile zone found at the surface for the circular point contact. It leads to only two zones, one in compression from the surface to a depth approximately two times the Hertzian depth, followed by a tensile zone at a higher depth. This mean stress varies from −390 MPa to 180 MPa at 8 GPa.
Finally, Fig. 19 presents the permanent print found at the center of the track along the rolling direction for different normal loads. Note that coordinates x and z are given in a dimensionless form, i.e., divided by the semiminor axis a found at 4.2 GPa ͑i.e., 310 m͒. Note that the maximum depth is found for the maximum load ͑8 GPa͒ and corresponds to 6.93 m near the starting point and to 2.96 m in the stationary region.

Conclusion
Mechanical components of the high level of reliability are usually designed to operate with a low level of stress, ideally below Fig. 13 Profile of the von Mises stress found under loading at the center of the track and in the steady-state region for various ellipticity ratios and after one and two cycles "P Hertz = 5.7 GPa, ceramic ellipsoid "elastic…/flat 52100 "elastic-plastic… surface… Fig. 14 Profile of the residual von Mises stress found after unloading at the center of the track and in the steady-state region for various ellipticity ratios and after one and two cycles "P Hertz = 5.7 GPa, ceramic ellipsoid "elastic…/flat 52100 "elastic-plastic… surface… the yield stress. A contact overload may, however, be encountered accidentally or may be the result of a transient operation regime. In other components operating under severe or extreme conditions, the yield stress is often exceeded during the normal running conditions.
Prior to studying the effect of the material hardening on the fatigue life, it is interesting to describe the effects of such an overload on the contact behavior and residual state. The EP response of a half-space normally loaded by a rolling elastic sphere  . 15 Hydrostatic pressure found after unloading at the center of the track and in the steadystate region for various ellipticity ratios and after one and two cycles "P Hertz = 5.7 GPa, ceramic ellipsoid "elastic…/flat 52100 "elastic-plastic… surface… Fig. 16 Profile of the equivalent plastic strain found at the center of the track and in the steady-state region for various normal loads and after one and two cycles "k = 8, ceramic ellipsoid "elastic…/flat 52100 "elastic-plastic… surface… or ellipsoid has been investigated, and results are presented. The effect of the ellipticity ratio and contact pressure on the surface permanent deformation and subsurface stress and strain states has been analyzed. For a circular point contact, it has been shown that the permanent surface print is deeper for a rolling load than for a point experiencing vertical loading and unloading. The steady-state regime is reached after a very few number of cycles when using an isotropic hardening model, typically 2 or 3. Moreover, a tensile zone lying on the surface of the EP body on both sides of the rolling track is found.
For an elongated contact with a semimajor axis perpendicular to the rolling direction, a shallower but wider surface print and a lower level of residual stress and strain is found when increasing the ellipticity ratio while maintaining the maximum Hertz pressure constant. This point should be carefully considered when extrapolating experimental fatigue life data obtained on any analytical ball on a flat test apparatus to real rolling bearings.
Less surprising is the effect of the contact pressure on the surface and subsurface stress and strain for an elastic-plastic contact with an ellipticity ratio of 8. Residual stress and strain are found to increase with the normal load, up to 2.1% for the equivalent plastic strain and 935 MPa for the equivalent von Mises stress for the AISI 52100 bearing steel with a frictionless rolling load of 8 GPa. Fig. 17 Profile of the von Mises stress found after unloading at the center of the track and in the steady-state region for various normal loads and after one and two cycles "k = 8, ceramic ellipsoid "elastic…/flat 52100 "elastic-plastic… surface… Fig. 18 Profile of the hydrostatic pressure found after unloading at the center of the track and in the steady-state region for various normal loads and after one and two cycles "k = 8, ceramic ellipsoid "elastic…/flat 52100 "elastic-plastic… surface… Fig. 19 Displacement at the center of the rolling track after unloading for various normal loads and after one and two cycles "k = 8-ceramic ellipsoid "elastic…/flat 52100 "elastic-plastic… surface…