1 Introduction

Rolling element bearings are designed to enable smooth rotation during service. Severe loadings, starved lubrication condition and dynamic loadings cause detrimental effects on bearing life [1]. Stress histories in a bearing element during rolling cycles are termed rolling contact fatigue (RCF). The two regimes of fatigue high cyclic fatigue and low cyclic fatigue are distinguished due to the number of running cycles required for failure and the rate of irreversible plastic flow which accumulates with corresponding stress cycles. RCF is, however, different from the conventional low cyclic and high cyclic fatigue due to complex multiaxial nature and non-proportional loading [2]. Typically in a rolling bearing element, carbide volume fraction, retained austenite, microtexture, microstructure and deep zone residual stresses with non-proportional cyclic stress histories make bearing material heterogeneous in nature [3]. A generic representation of a typical bearing life exhibits three phases of microstructural changes and residual stress evolution during RCF loading starting from initial plastic shakedown (stage I) to steady-state response (stage II) which is followed by the instability (stage III) leading to bearing failure [4]. Initial plastic shakedown refers to the initial accumulation of plasticity during RCF cycles; afterwards, there is little to no further plasticity generation and bearing operates mostly in the elastic regime.

According to Hertzian contact theory, the maximum pressure (Pmax) required to cause plastic yielding of a subsurface material point with a static contact load is given by Pmax = 1.6 σy, where σy is the yield strength of the virgin material [5]. Rolling bearing elements are typically operated in <2.5 GPa pressure so the loading response of the material remains globally elastic. However, at microstructural level, the local stress concentration due to carbides and non-metallic inclusions causes subsurface plasticity, also named microplasticity [6]. However, the range of plasticity accumulation (∆εp) with progressive RCF cycle beyond the shakedown limit is minimal, i.e. in the range of ∆εp<0.003 and is spread over a wide range, i.e. ~1010 cycles [7]. During this steady-state stage, the microstructure of the stressed zone under Hertzian contact evolves gradually, which causes material phase transformation, changes in residual stresses and change in yield strength unceasingly. The ability of the material to work harden under such conditions is quite significant and the spatial and temporal evolution of this hardness change is very critical for damage accumulation as hardness change defines the material resistance to plastic flow. The continuously evolving cyclic stress–strain fields under the influence of residual stresses makes complexities in modelling rolling contact fatigue. Despite the wide usage of high strength bearing steels, limited information on material hardening is available for standard AISI 52100 bearing steel only under monotonic tension and compression loading and static torsional loading [7] which cannot address the true micromechanical response of bearing material operated under triaxial loading in RCF. Due to the absence of relevant material hardening data, it is quite challenging to undertake realistic studies on RCF-induced damage evaluation and resulting residual stress estimation.

The residual stresses during RCF arise due to non-uniform volume expansion (decomposition of retained austenite) and non-uniform plasticity due to complex subsurface stress state. Deep zone residual stresses formed during cyclic loading ranging from 3b to 9b contact depth play a vital role in bearing life [8], where b represents the semi-minor axis. Residual stresses are especially valuable when employed with higher strength metals. It has been reported [9,10,11] that the fatigue life of a bearing element can be improved with the help of compressive residual stresses (CRS). These compressive residual stresses significantly affect fatigue life by crack closure and hindering crack growth [11, 12]. Moreover, it also alters the alternating applied stress range which governs the fatigue life of bearing material. These residual stresses can be decreased or totally lost by stress relaxation which is associated with the cyclic yielding and higher strain amplitude. However, at very small strain amplitude (during steady-state stage), the relaxation ceases at stable value [13, 14]. Voskamp [9], after endurance testing of deep groove ball bearings (DGBB), reported onset and progressive change in residual stresses and structural alterations with rolling cycles. The gradual structural alterations have been termed Dark Etching regions (DERs) and White Etching bands (WEBs) [2, 15, 16] owing to decomposition of retained austenite and parent martensite. Later experimental studies [10, 17, 18] on ball bearings owing to rolling contact fatigue showed significant contribution of deep zone residual stresses towards bearing failure. Moreover, researchers [19, 20] have reported that a close proximity exists between residual stress distribution with subsurface microstructural alterations; however, the assessment of residual stresses during subsurface alterations along with plasticity accumulation has not been evaluated.

A realistic finite element model to RCF requires prior knowledge of cyclic response of material along with a reliable material model. As discussed earlier, during rolling contact fatigue, material undergoes cyclic plasticity which is accumulated at subsurface level over hundreds of million cycles along with the kinematic translation and radial expansion of yield surface. It is reported [21, 22] that the plastic strain accumulation during RCF cycles results in contact track widening and change in contact geometry. Previous research studies [5, 22,23,24,25,26] employed the power-law hardening model and elastic linear kinematic hardening plastic (ELKP) model where the hardening parameters were derived from uniaxial compression test or torsional test of standard bearing steel. These models described the plastic distortion, cyclic plastic strains and residual stresses in early life (RCF cycle < 106) of a bearing in good agreement with the experimental data. However, it under-predicted the residual stresses in the later life (RCF cycles > 106) of a bearing element, as reported by Bush [27] and Voskamp [9]. Besides, these conventional models failed to incorporate the stochastic nature of the bearing material due to cyclic yielding during RCF which needed to be defined with a more complex and sophisticated hardening model.

The non-linear isotropic and kinematic hardening (NIKH) model was presented by Chaboche and Lamaitre to study the cyclic response of material under multiaxial fatigue [28]. The NIKH model is capable of simulating complex material behaviour during RCF loading, i.e. plastic shakedown, cyclic hardening, and ratcheting in the steady-state regime of RCF cycles [29]. It has been reported [30] that the ratcheting based response of the bearing material microstructure can be modelled with the help of non-linear kinematic hardening response. Nevertheless, the initial response of the bearing in the early stage is best represented by the combined effect of isotropic expansion and kinematic expansion. Current study also employs the combined isotropic and non-linear kinematic hardening model with the combination of nanohardness measurement. The reason for opting such material behaviour is to successfully model the initial plastic shakedown along with the ratcheting behaviour of bearing material. The material constants required for NIKH material model can be derived with the help of expanding cavity model presented by Gao et al. [31]. The expanding cavity model describes the indentation deformation of elastic power-law hardening and elastic linear hardening materials and reported that, for given indenter geometry, indentation hardness depends on Young’s Modulus, yield stress and strain hardening index ‘n’ of the indented surface. Earlier studies [32, 33] have employed micro-indentation experiments on RCF damage evaluation and reported a gradual increase in material hardness with the linear accumulation of plastic strains. However, they were unable to evaluate the very localised damage of the RCF damaged zone due to the limitation of micro-indenter that only gives an average hardness over the grain [34]. Current research work overcomes that problem with the use of a Berkovich nanoindenter. The use of three-faced pyramidal Berkovich indenter having the same depth-to-area ratio as a Vickers indenter made it convenient to measure accurately the change in hardness of the RCF-affected zone. The precise temporal and spatial evaluation of this hardness change due to RCF-induced damage resulted in cyclically evolved yield stress which can be used to estimate the cyclic hardening parameters for NIKH material model. The resulting model can accurately capture the evolved elastoplastic response of the material in the later life of RCF loading and resulting residual stresses can be assessed.

The following sections describe RCF experiments with a ball-on-ball rolling contact in a rotary tribometer which was followed by surface/subsurface investigations. Next, development of a 3D finite element model for RCF is presented. The material model relies on the constitutive response of the bearing material after RCF testing. The constitutive response is based on the cyclic hardening parameters obtained from nanoindentation in conjunction with expanding cavity model. Simulation results showed a significant difference in estimation of deep zone residual stresses employing newly developed semi-empirical material model when compared with the conventional uniaxial power-law hardening in the later life of rolling cycles.

2 RCF Experimentation

A high-speed Microprocessor Rotary Tribometer (TE-92) with 4-ball configuration was employed to run the RCF test on through-hardened AISI 52100 bearing steel balls of 12.7 mm diameter. In this configuration, three lower balls are kept in a cup which acts as outer race of the bearing. The upper ball is mounted on a shaft with the help of a collet and rotates on top of lower three balls which forms a contact track similar to the inner raceway of a bearing [35, 36], following Fig. 1. Industrial grade ISO-VG-320 (kinematic viscosity of 300 cSt at 40 °C) was used as a lubricant to avoid asperity contacts so that pure subsurface damage can be evaluated with no contact track wear. The cleanliness of the added lubricant in the lower cup and calculated lambda ratios (described in section Appendix) ensured no wear phenomenon during the test. Tests were carried out with two different load levels at the spindle with 40 °C lubricant temperature which was kept constant throughout experiments. Tests were stopped either at required RCF cycles or whenever a ball fails due to spalling which was detected by a vibration cut-off sensor. Schematic for 4-ball configuration is shown in Fig. 1 and conducted experiments are listed in Table 1. After experimentation, samples were sliced and then polished with diamond suspension to make the sample suitable for nanoindentation. A 3% Nital solution (97% Ethanol, 3% Nitric acid) was also used later on to reveal the surface as an etchant. A uniaxial compression test was also carried out to obtain the flow curve of the bearing material with similar composition as of RCF-tested samples. The stress–strain curve is plotted in Fig. 2 and corresponding material properties are referred in Sect. 2.2 for subsurface hardness evaluation.

Fig. 1
figure 1

Four-ball test configuration shown with 12.7-mm-diameter bearing balls; contact track profile is shown in the left top corner. Post-RCF-tested sample sliced in the circumferential direction (parallel to contact track) for subsurface analysis

Full size image
Table 1 Four-ball experiments conducted for RCF test under high and low load conditions
Full size table
Fig. 2
figure 2

Flow curve obtained from the compression test of bearing steel at room temperature with similar composition as RCF-tested samples

Full size image

The 4-ball kinematics is described as the upper ball (with radius 6.35 mm) rotates on top of lower three balls having same radius, the upper ball in collet undergoes 2.25 times stress cycle and develops a contact track which continues to grow with respective stress cycles. This contact track was investigated with the help of a 3D optical interferometer for each sample during RCF test. The samples were cleaned with the ultrasonic cleaner to remove any lubrication prior to interferometry and measurements were taken after specified running cycles. A typical interferometer analysis is shown in Fig. 1. The consistency and repeatability of the contact track measurements indicate the reliability of the acquired data.

2.1 Hardness Measurement

Subsurface hardness was evaluated for the virgin material with a Vickers Indenter (HV 0.3). For through-hardened bearing steel, it was expected to be nearly constant over entire depth. After RCF experimentation, samples were sliced in a plane parallel to contact track, see Fig. 1, and polished with a final diamond suspension spray of 0.5 μm particle size to remove unwanted surface roughness. For nanoindentation, a 3-faced pyramidal Berkovich Indenter was used which measures the indentation hardness (HIT) as a function of load–displacement curve. The indentation hardness (also known as Instrumentation hardness) is defined by the applied load (F) divided by the projected area of indent (A) as

$${H}_{IT}=\frac{F}{A\left({h}_c\right)}.$$
(1)

The projected area of indent is a function of contact depth (hc) and is calculated from the knowledge of indenter geometry and stiffness of the contact. To determine the projected area function of a typical Berkovich indenter, a reference sample with known elastic modulus is indented and a relationship is obtained between contact depth and projected area. Further details of area function analysis are described in section Appendix.

The load–unload curve from nanoindentation is fitted to power-law fit and contact compliance and plastic depth (hc) is calculated using Oliver and Pharr analysis [37]. Once the plastic depth is computed, the instrumented hardness can be obtained in GPa from equation (1). From ISO 14577-1:2015, the ratio of surface area to projected area is kept constant for a perfect Vickers or modified Berkovich indenter, which leads to the correlation of indentation hardness with the conventional Vickers hardness as

$$HV=0.0945\times {H}_{IT}\ \left( in\ MPa\right).$$
(2)

The instrumentation hardness obtained from the nanoindentation test for each sample can be converted in Vickers hardness scale. Fig. 3a reveals the Nital-etched microstructure after 6 GPa test and Fig. 3b represents nanoindentation performed in rows from the surface towards the core of sample. The applied load is kept constant in each row of indent and corresponding instrumentation hardness is evaluated from equation (1) and (2). Fig. 3c shows the nanoscale measured hardness in HIT and respective estimated Vickers hardness (HV) for both samples. For comparison of the converted hardness scale, an HV0.3 hardness profile is also plotted which shows consistent hardness response from virgin ball bearing sample. It is evident that the baseline hardness for the nanoindentation profiles can be marked near 9.01 ± 0.1 GPa (~ 853 HV in converted scale) which is reasonably in range with the measured HV0.3 values. It has been observed that changing the nanoindentation load greatly affects the observed hardness where the baseline hardness shifts its position. In order to compensate the size effect of applied load, the nanoindentations were conducted with applied loads of 100, 50 and 10 mN and series of tests were carried out to check the repeatability and reliability of the measured hardness.

Fig. 3
figure 3

Subsurface investigation of axially sliced RCF-tested samples: (a) Nital-etched sample after 42.7 million RCF cycles at 6 GPa contact stress, (b) nanoindentation performed at 100, 50 and 10 mN; zoomed-in area represents Berkovich indent lies in the vicinity of carbides, (c) indentation hardness HIT, plotted for both RCF-tested samples along with reference virgin material hardness HV0.3

Full size image

Considerable microstructural alterations (Dark etching regions, DERs) can be observed under the contact track, marked as RCF-affected zone in Fig. 3a. Such microstructural alterations are attributed to the decay of parent structure due to the subsurface plastic deformation during high cyclic loadings which accumulates over millions of RCF cycles [9]. The associated bearing material response can be characterised by an increase in highly localised hardness with reference to the observed baseline hardness. No observable change in the microstructure was obtained for 4 GP tested sample; however, a slight change in subsurface hardness can be observed.

For elastic power-law hardening materials, Gao [31] presented an expanding cavity model which correlates the material hardness with elastic modulus and flow stress as

$$\frac{H}{\sigma_y}=\frac{2}{3}\ \left\{\left(1-\frac{1}{n}\right)+\left(\frac{3}{4}+\frac{1}{n}\right){\left(\frac{1}{3}\ \frac{E}{\sigma_y}\ Cot\ \alpha \right)}^n\right\},$$
(3)

where ‘H’ represents hardness in MPa, ‘n’ is strain hardening exponent, derived from monotonic compression test, Young’s Modulus E= 210 GPa and ‘α’ is effective cone angle (α=70.2996° for Berkovich Indenter).

The highly localised nanoindentation hardness estimated in Vickers hardness scale can be transformed into evolved flow stress of the material using equation (3). As discussed earlier, expanding cavity model works well with the strain hardening materials, so an appropriate value of strain hardening exponent was derived from uniaxial compression test (see Fig. 2), and corresponding flow stresses were evaluated under the contact track. Maximum evolved flow stress obtained at 155 μm below contact path was 2967 MPa (HIT = 10.47 GPa), which corresponded to a total hardness increase of ~417 MPa (2967 – 2550 = 417 MPa). This 417 MPa flow stress was evolved during RCF loading at 6 GPa contact stress after 42.75 million RCF cycles. This increase in hardness is in comparison with the previous studies [33, 38] where substantial hardening of subsurface region during RCF loading was reported. The precise measurement of change in hardness for RCF sample tested at 4 GPa contact stress also reported an increase of hardness (i.e. 9.91 GPa) at 95 μm depth with reference to the baseline hardness. This RCF-affected evolved flow stresses represent the evolved yield stress of the subsurface in the later life rolling cycles and can be used as an input parameter for the finite element model.

3 Finite Element Modelling

To simulate the RCF loading of a 4-ball test, a 3D finite element model, shown in Fig. 4, has been developed in a commercial solver Abaqus. The model comprises of a deformable sphere with rigid plates rotating over the contact track where the whole model was designed to the original scale. A 4-node linear hexahedral element type has been used to model the entire contact track with relatively finer mesh size, while the rest of the geometry was meshed with quadratic tetrahedral elements. Use of cyclic symmetric model, rigid surfaces and optimised mesh element size (5 um for hexahedral brick element) made the model possible to solve with higher accuracy and less computational time. The load was applied at an angle of 35.26° with the vertical axis (as demonstrated in Fig. 1) and was ramped during multiple cycles so that stresses can be stabilised and smoothly distributed across the contact track. For interacting bodies, a surface to surface contact pair was defined with a frictionless contact.

Fig. 4
figure 4

FE Model for simulating rolling contact fatigue in a 4-ball test

Full size image

The non-linear kinematic isotropic material model was employed to mimic the cyclic response of AISI 52100 bearing steel. This model is based on two types of material behaviour, i.e. an isotropic part which causes the yield surface of material to expand and a non-linear kinematic part which causes the yield surface to translate with a back stress. The model is based on von Mises or maximum distortion energy yield criterion. The generalised associative flow [28] is given as

$$f\left(\text{S},\text{q},{\varepsilon}_{pl}^{-}\ \right)=\left|\ S-q\ \right|-{\sigma}_y,$$
(4)

where ‘S’ is the stress vector, ‘q’ is the back-stress tensor, \(^{\prime }{\varepsilon}_{pl}^{-\prime }\) is the accumulated plastic strains and ‘σy’ is the instantaneous yield strength. For non-linear isotropic hardening, this instantaneous yield strength will be

$${\sigma}_y={\sigma}_y^o+Q\left(1-{e}^{-\upbeta\ {\varepsilon}_{pl}^{-}}\right).$$
(5)

Here ‘\({\sigma}_y^o\)’ is initial yield strength of the material, and ‘Q’ represents the maximum expansion of yield surface with rate ‘β’. For non-linear kinematic hardening,

$${q}^{\prime }=\frac{2}{3}\ C\ {\varepsilon}_{pl}^{.}-\gamma\ {\varepsilon}_{pl}^{.-},$$
(6)

where ‘q′’ denotes increment in back-stress tensor, ‘C’ is initial kinematic hardening modulus, ‘γ’ is rate of decreasing ‘C ‘and ‘\({\varepsilon}_{pl}^{.}\)’ is incremental equivalent plastic strains. The kinematic term added in the model as AISI 52100 steel represents ratcheting (continuous accumulation of plasticity) during cyclic loading. In order to stabilise the ratcheting response, isotropic term was added so that a realistic response towards saturation can be obtained.

Recall from Fig. 3d that the virgin material hardness can be obtained at a depth far off Hertzian contact which was found to be around 853 HV (9.01 ± 0.1 GPa), and the maximum change in instrumented hardness for 6 GPa and 4 GPa RCF test was found to be 1.4714 GPa and 0.789 GPa which yielded approximately 417 MPa and 225 MPa flow stress. The change in flow stress under a given loading condition is attributed to the maximum change in the yield surface, ((∆σy)max). This maximum change in yield surface can be integrated into the isotropic and kinematic hardening and expressed as

$${\left(\Delta \sigma y\right)}_{max}\Big)=Q+\sqrt{\frac{2}{3}}\frac{C}{\gamma }.$$
(7)

Equation (7) represents the maximum possible change in the yield surface under RCF loading and corresponds to the saturation of cyclic hardening. This total change in yield surface can be split into isotropic expansion and kinematic expansion equally as an initial assumption since no accurate cyclic data are available for AISI 52100 bearing steel as discussed earlier. This initial assumption would lead to a 50% expansion and 50% translation of yield surface. Next step is to estimate the cyclic hardening parameter ‘C’, which is the initial kinematic hardening modulus. From compression stress–strain curve, Fig. 2, it can be noted that once the stress in a material exceeds the yield limit, the material experiences elasticity to plasticity transition state and its elastic modulus changes (E) to the plastic modulus ‘M’. It is known that with an increase in the plastic strains, ‘M’ tends to reduce, regardless of its stress state (i.e. uniaxial or triaxial) and eventually comes to zero representing a saturation in the strain hardening. However, during transition stage from elasticity to plasticity, the rate of plastic strain accumulation is very small and plastic modulus (M) is not much different from the elastic modulus (E). Therefore, it is safe to say that the parameter ‘C’ which is the initial kinematic hardening modulus can be set to a value slightly smaller than E = 210 GPa. With a continuous increase in the plasticity, the kinematic hardening modulus keeps on decreasing and is being governed by the factor ‘γ’. From equation (7), if the rest of the parameters are known, ‘γ’ (rate of decreasing ‘C’) can be easily computed.

As the total change in yield strength has been divided into isotropic and kinematic part equally, therefore it is vital to define the expansion of yield surface (Isotropic hardening) which is controlled by the hardening parameter ‘β’. This rate parameter only affects the number of cycles required for isotropic hardening hence changing the value of ‘β’ does not affect the overall hardening of the material. However, to reduce the number of RCF cycles for lesser simulation run time, a constant value of ‘β’ has been fixed initially. As an initial assumption, 50% isotropic and 50% kinematic hardening contribution was agreed to the total extent of hardening and corresponding cyclic hardening parameters are listed in Table 2. Previous studies [30, 39] have shown that in combined hardening models, 90% of hardening is associated with kinematic hardening. However, in kinematic hardening mode, multiple back-stress tensors are required which creates complexities during parameters estimation. The initial assumption of combined hardening in the current study provides preliminary understanding due to the absence of any known cyclic hardening parameters. The cyclic hardening parameters for 4 GPa RCF test are also estimated and given in Table 2.

Table 2 Cyclic hardening parameters listed for Kine50% as an initial estimation
Full size table

4 Result & Discussion

The newly developed Finite Element Model (FEM) simulates the RCF stress cycles as the load moves from right to the left end of the domain, as shown in Fig. 5. Simulations were conducted employing a semi-empirical NIKH material model (cyclic hardening parameters listed in Table 2) and power-law hardening model (uniaxial compression flow curve, Fig. 2). The material model follows the J2-based plasticity criteria. As load moves from right to left in the 3D domain, marked as RD (rolling direction), the material undergoes yielding resulting in a subsurface plasticity zone at the centre of contact track. The initial RCF cycles have been simulated to ramp the load gradually so that a uniform stress distribution can be obtained under the contact track.

Fig. 5
figure 5

(a) Direction of loading in 3D domain and plasticity displayed in a cylindrical coordinate system denoted as 11: radial, 22: circumferential and 33 as tangential direction. (b) Octahedral shear stress represented in a circumferential plane sliced through the centre of contact track

Full size image

The plastic strains at subsurface lead to the plasticity-induced residual stresses in a 3-dimensional state. With repeated RCF cycles, the cyclic stress–strain loop stabilises and reaches a steady-state indicating plastic shakedown. Fig. 6a and b represents a stabilised hysteresis loop obtained for power-law hardening model and NIKH material model, respectively. Each hysteresis loop is plotted for shear stress and shear strain components at a representative depth of maximum stress state. The number of cycles required to achieve the steady-state response of the hysteresis loop depends upon the loading behaviour and the cyclic hardening parameters. For the semi-empirical numerical model used in this study, the stabilised hysteresis loop was obtained after 5 contact cycles, whereas, for power-law hardening model, it was obtained after 3 contact cycles. The reason for such early plastic shakedown at 6 GPa is attributed towards the initial assumption of 50% isotopic/kinematic hardening contribution.

Fig. 6
figure 6

Shear stress–strain cycle obtained for (a) power-law hardening and (b) NIKH material model; the steady-state hysteresis loop is represented by black solid lines

Full size image

Result obtained after a few RCF cycles were acquired for a high load of 857 N and are presented in Fig. 7a and b. Fig. 7a represents the von Mises stress plotted with equivalent plastic strains for multiple load cycles operated under NIKH material model. The von Mises stresses continue to increase with increasing load and then start to saturate after a few load cycles as defined by the hardening parameters of the constitutive model. It can be noted that after first rolling pass, the area under contact track (subsurface region) accumulates a considerable amount of plasticity and then it starts to saturate once the stabilised cyclic stress–strain loop is achieved. A stabilised cycle is obtained just after four rolling cycles, where no further plasticity is generated under same RCF loading. The von Mises response (Fig. 7a) during rolling passes shows the generation of residual stress along with plasticity which cumulates with each rolling pass and its build-up tend to cease. The accumulation of plasticity in each load cycle in the presence of mean stress is termed ratcheting. The ratcheting phenomenon causes residual stress to relax. However, beyond shakedown limit the stress relaxation is minimal due to low plasticity damage as discussed earlier in the previous section, until softening of the material happens after hundreds of millions of RCF cycles. However, due to the limited scope of this research, such material response has been excluded in modelling and will be added as future work. Fig. 7a reports the potential generation of residual stress (RS) below the contact path which was found to be distributed below the contact path in the circumferential and tangential direction; however, the maximum of the RS distribution was observed at centreline contact path. Fig. 7b depicts the maximum orthogonal shear stress experienced by an element at a depth z/b = 0.48 (represents normalised centreline depth) with plastic strains growth. It is important to know that the plastic strains accumulate in each cycle where reversal of orthogonal shear stress takes place until shakedown and there is little to no accumulation of plasticity in the following cycles.

Fig. 7
figure 7

Accumulation of plastic strains with (a) von Mises stress (b) Maximum orthogonal shear stress at representative depth below contact track

Full size image

After load removal, the residual stresses and contact track geometries were investigated. The contact track results obtained from the finite element model and 4-ball test are listed in Table 3 for comparison. Contact track semi-major axis observed from interferometer at both loading condition was found close to the NIKH material model in comparison with elastic (Hertzian theory) and elastoplastic model (power-law hardening). Table 3 also shows the maximum contact stress obtained from analysis which was found to be nearly the same for low resultant load of 265 N; however, a significant difference was noticed at a high resultant load of 857 N. Power-law hardening model and NIKH hardening model showed 5605 MPa and 5414 MPa, respectively. The difference in the obtained contact stresses is attributed to the contact track plastic deformation and resulting track widening. The obtained half-width radius ‘a’ from numerical and experimental findings is listed where the value found by NIKH model was observed closer to the experimental results.

Table 3 Contact track numerical and experimental analysis
Full size table

During RCF cycles, bearing material undergoes multiaxial loading, resulting in subsurface plastic damage and contact track widening. This damage is being accumulated over millions of cycles and can be related to early life and later life of the bearing material. Contact track profile measured after RCF test results from cyclic deformation due to RCF loading. A 2D profile acquired normal to the contact track represents the maximum change in contact track profile (width and depth). Fig. 8a represents the interferometer results obtained for 6 GPa contact stress after 0.37 million and 42.75 million RCF cycles and compared with the FEM results. It can be seen that the contact track profile obtained after 0.37 million cycles corroborates with the power-law hardening model representing a potential early life plastic shakedown of contact track. Moreover, the contact track profile obtained after 42.75 million cycles corresponds with the combined NIKH hardening model which endorses that the combined hardening model can better estimate the later life bearing material elastoplastic behaviour. In order to verify this argument, the newly developed finite element model was extended to the 4 GPa contact stress with corresponding hardening parameters listed in Table 2. The plastic deformation at contact track was obtained once the elastic-plastic shakedown is attained in the simulation. Fig. 8b represents the contact track profiles for experimental analysis as well as simulation from combined material model. It can be seen that the experimental curve fitted well in the range of the experimental observation after 37.5 million RCF cycles under 3980 MPa true contact pressure (Table 3). The close agreement and consistency of plastic deformation at contact track observed from experimental and simulation provides a linkage to the RCF subsurface damage and resulting continuously evolving residual stresses generation in the later stages of RCF life cycles.

Fig. 8
figure 8

Contact track profiles obtained from FEM models and RCF test conducted at (a) 6 GPa and (b) 4 GPa contact stress

Full size image

Once the cyclic hardening parameters are testified for combined hardening model, the next step is to evaluate the cyclic hardening parameters which can best describe the evolved elastoplastic behaviour of bearing steel. The initial assumption of 50% kinematic hardening contribution in the total hardening can be varied and resulting contact track plasticity can be analysed. Table 4 lists a change in kinematic hardening contribution to 10%, 50% and 70% of the total hardening and resulting values of Q, C and γ are updated. Fig. 9a represents the contact track plasticity obtained from the variation in cyclic hardening parameters and compared with the experimental data. A 10% contribution of the kinematic hardening results in a maximum contribution of isotropic hardening which resists plasticity by a linear expansion of yield surface. Resulting contact track profile is plotted as Kine10% and is found quite low as compared to the experimental curve, whereas a 50% and 70% kinematic hardening curve lies in the vicinity of experimental data. It should be noted that these simulations are carried out with a fixed value of β = 50, which defines the rate of isotropic expansion. Changing the β value for current simulations affects the cycles required for plastic shakedown. A lower β = 15 results in slightly higher plasticity accumulation (PEEQ) as observed in Fig. 9b, however, requires a considerably higher computational time in terms of stress cycles required for plastic shakedown. It was noted that a large value of β results in lower computation time with a maximum of 5% variation in contact track deformation.

Table 4 Evaluation of cyclic hardening parameters
Full size table
Fig. 9
figure 9

Evaluation of cyclic hardening parameters: (a) Contact track profile with variation in parameters obtained for 6 GPa contact stress. (b) Effect of changing β at subsurface plasticity accumulation

Full size image

Fig. 10 represents the residual stress profiles obtained from experimental and numerical observation of through-hardened bearing steel with similar composition. Residual stress component, σ22, (stress component in the circumferential direction, Fig. 5a) was acquired from 4 GPa contact stress simulations employing power-law hardening and NIKH material model (parameters listed in Table 4) and compared with the residual stresses data presented by A. P. Voskamp [9]. The extensive testing data obtained from the fatigue testing of deep groove ball bearing (DGBB) using X-ray diffraction technique are considered to be a benchmark for residual stress estimation under RCF load. From Fig. 10, power-law hardening model showed a reduced level of the residual stress generation during rolling cycles when compared with experimental curves obtained at 1, 10 and 100 million rolling cycles, whereas the residual stress results from combined NIKH material model lies closer to the experimental data. The underprediction of residual stresses is attributed towards the elastic-plastic response caused by the power-law hardening material model which cannot mimic the cyclic hardening during rolling cycles resulting in lesser yielding of the subsurface. It can be noted that the residual stress maxima was found dissimilar for both material models. It is important to note that the residual stress data [9] plotted in Fig. 10 are reported from standard DGBB of 6309 type. The non-conformal contact in ball bearing at ball-inner raceway makes an elliptical contact which is different from a point contact of a 4-ball configuration. The difference in the osculation of contact geometries between simulated and reported experimental makes this comparison challenging which has been addressed by the normalised centreline distance. The non-linear kinematic hardening model modified with the expanding cavity model employing nanoindentation yielded a much more plausible residual stress estimation. It should be noted the NIKH residual stress profile presented in Fig. 10 corresponds to a test with 37.5 million RCF cycles at 3980 MPa maximum contact stress so it could fall somewhere between 10 million and 100 million RS profiles obtained from DGBB test. It should be noted that the near-surface higher residual stresses in experimental data are associated with the machining processes involved to manufacture the bearing samples and were ignored during finite element analysis.

Fig. 10
figure 10

Residual stress estimation from finite element models and experimental data* obtained from RCF testing of DGBB 6309 [9] after 1, 10, 37.5 and 100 million cycles

Full size image

As discussed in Sect. 1, the gradual changes in residual stress during rolling cycles can be correlated to the changes in yield strength as well as to the changes in specific volume arising during microstructural alterations. The decomposition of retained austenite is associated with a volume expansion, whereas the decay of martensite is related to the volume contraction; hence, the total effect can be considered as the relative volume fraction of these phase changes. Although the microstructural phase changes progress during RCF, the bearing material undergoes considerable work hardening in the subsurface plastic zone as reported in 4-ball experimentation. The relative increase in hardness observed at 6 GPa and 4 GPa contact stress supports the consistency of this hardness change with progressive RCF cycles. The current study uses the change in subsurface local hardness as an indicator to estimate the gradual change in residual stress, which was not possible earlier due to the complexities of the rolling contact fatigue. Fig. 10 also shows that the residual stress profile predicted by the combined NIKH material model slightly deviates from the given experimental data. The main reason for such behaviour is the limitation of current material model which cannot take into account the decay of retained austenite resulting into internal compressive stresses and makes the residual stress profile broader at greater depths. Future work is recommended to integrate the presented FE model with the decay of retained austenite which can estimate the further development of residual stresses due to material phase transformations. Nevertheless, the close evaluation of the predicted residual stresses profiles can be seen as more profound towards greater depth at higher contact stress levels where compressive residual stresses tend to widen at greater depth. However, the normalised centreline contact path showed maxima for low and high resultant load at representative depth as shown in Fig. 11, where ‘z’ represents centreline distance from contact track and ‘b’ represents the semi-minor axis.

Fig. 11
figure 11

Later life residual stress estimation from modified NIKH material model with corresponding true contact pressure

Full size image

Recall from Sect. 1 that compressive residual stress generation under cyclic loading influences the fatigue life by altering the alternating stress amplitude especially in the later stage of RCF. To understand the change in alternating stress range driving RCF, the Mises stress profiles at centerlines contact track have been plotted for 4 GPa and 6 GPa contact stress in Fig. 12a and b. Maximum Mises stresses during loading state and after unloading have been shown for NIKH material model, power-law hardening and Hertzian elastic theory. The pure elastic contact between ball-on-ball point contact for 4 GPa contact stress has shown maxima of 2650 MPa which is slightly reduced at 2430 MPa employing power-law hardening model, representing initial plastic shakedown. However, the NIKH material model employed in current study has shown reduced Mises maxima of 2170 MPa. The 11% reduction in Mises stress at critical depth, i.e. z/b = 0.48, represents the evolved alternating stress (Mises) range after 37.5 million RCF cycles. Similarly, at 6 GPa contact stress, the maxima for applied Mises stress can be found as 3930 MPa, 2700 MPa and 2040 MPa for elastic, power-law and NIKH material model, respectively. The reduction in applied stress was found nearly 24% for NIKH material model. Moreover, it can be noticed that the evolved von Misses stress with the superposition of residual stress tend to shift its maxima from the critical depths to deeper or closer towards the surface. It was reported [9, 40] that the subsurface microstructural alterations, i.e. dark etching regions (DERs) and white etching bands (WEBs), are directly associated with the accumulation of plasticity which in turn related to the applied alternating stress range; the evolution of such Mises stress profile also helps to retard the development of visible microstructural alterations. Recent studies [34, 36] have shown that microstructural alterations are governed by the maximum distortion energy yield criteria during multiaxial RCF loadings where subsurface plasticity is governed by equivalent von Mises stresses. Inclusion of deep zone residual stresses with von Mises stress profile can better explain the formation mechanism of such microstructural alterations. The continuous development and progression of DERs and WEBs for ball-on-ball point contact under 4 GPa and 6 GPa contact stress are obtained from [36] and plotted in Fig. 12. The spread of dotted lines represents the development of minimum and maximum depth of microstructural alterations from a few million cycles to tens of million cycles. It can be seen that the wider spread of von Mises maxima with progressive RCF cycles corroborates well with formation of DERs which is comparable from a recent study [36]. With the subsurface plasticity accumulation, lowering of von Mises stress profile indicates lowering of yield stress which coincides with the depth range of WEBs formation. Fu et al. [40] has stated the ferrite bands within WEBs manifest the accumulated plastic strains (PEEQ) and its reported value is 0.01863 which is quite close to our current FEM results. It should be noticed that no WEBs were observed in current study due to low operating temperatures, i.e. 40 °C, as thermal diffusivity of Carbon can significantly affect the formation of WEBs as discussed in previous studies [36, 41]. Nevertheless, the evolved Mises criteria can define the evolution of microstructural alterations in well agreement with the experimental data [36] and can be used to develop formation models for DERs.

Fig. 12
figure 12

Simulated and numerical stress profiles acquired for (a) 4 GPa contact stress and (b) 6 GPa contact stress. Dotted lines represent microstructural alteration depth range under contact track, *obtained from 4-ball RCF testing [36], whereas the spread of these lines shows the continuous growth with progressive stress cycles

Full size image

It is evident from Fig. 12 that the later life residual stresses play a significant role in altering alternating Mises stress for each stress levels where the significant material strain hardening governs the evolution of residual stress. This evolution of residual stress and its effects on mean stress can be understood with a modified Goodman approach [42], where the residual stress affects the mean stress levels with the relation as \(\frac{{\text{S}}_{\text{max}}-{\text{S}}_{\text{min}}}{2}+\text{RS}\), where Smax and Smin represent the maximum and minimum amplitude of alternating stress in RCF. The presence of subsurface residual stresses along with mean stress is extremely beneficial and hence optimum microstructural design for significant work hardening of bearing material and resulting compressive residual stresses is a key for a prolonged service life of bearing elements. Current research work is focused on the semi-empirical approach to model RCF simulation and evaluating subsurface plasticity-induced residual stresses and highlights its importance to increase the fatigue life of bearing elements in operation and understanding microstructural alterations. The novelty of the current study lies in the methodology to carry out this scientific research where only four cyclic hardening parameters are required to accurately capture the evolved bearing material response. The combined hardening model (NIKH) employed in the current study incorporates the isotropic expansion and non-linear kinematic hardening of a bearing material operated under RCF. In order to obtain the cyclic hardening parameters, a three-faced pyramidal Berkovich nanoindenter has been used to evaluate the subsurface yield stresses in terms of instrumentation hardness (GPa) which evolves gradually with progressive RCF cycles. The precise spatial evaluation of the very localised hardness change and the integration of the expanding cavity model facilitates the analysis of evolved subsurface yield stresses which were then employed into FEM simulations. The constitutive material constants employed in the current study are obtained directly from the RCF-tested bearing samples which represent the true elastic-plastic nature of evolved bearing steel which makes the current model more feasible and viable to apply in industrial research and applications. Furthermore, the semi-empirical material model can be coupled with continuum damage mechanics to investigate the failure and degradation of the material. This is significant in terms of design for the durability of bearing elements.

5 Conclusions

A comprehensive finite element model has been developed for evolved bearing material elastoplastic response, prediction of plasticity-induced deep zone residual stresses and subsequent effects on microstructure. For RCF testing, standard AISI 52100 bearing steel balls have been tested in a rotary tribometer at 4 GPa and 6 GPa contact stress. Post-RCF subsurface analysis has been conducted to obtain the hardening response of the bearing material. This research concludes the subsequent key points:

  • The subsurface microscopic investigation has revealed substantial microstructural alterations at 6 GPa contact stress, whereas no considerable microstructural changes were observed at 4 GPa. Nanoindentation analysis of the RCF-tested samples revealed a significant change in subsurface hardness. For both contact stresses, substantial work hardening of bearing material was observed at representative depths.

  • The cyclic hardening parameters were estimated with instrumented hardness measurement of the evolved subsurface in conjunction with the expanding cavity model. The cyclic hardening parameters were estimated for combined non-linear isotropic kinematic hardening (NIKH) behaviour and employed in newly developed 3D finite element model to mimic the cyclic hardening of rolling contact fatigue.

  • The modelling results for contact track profiles and residual stresses revealed that the simplistic power-law hardening model can estimate the early plastic shakedown of RCF, whereas the combined NIKH hardening model corroborates well with the elastoplastic behaviour in the later stages of RCF stress cycles.

  • Microplasticity-induced residual stresses showed a significant effect on the alternating Mises stress range. The 11% and 24% reduction of Mises stress just after few tens of million cycles shows the evolution of subsurface material response under RCF. The inclusion of deep zone residual stress with the von Mises stress criteria explain the progression of microstructural alterations and can be extended to develop the prognostic models for DER/WEB formation.

  • The quantitative model for plasticity-induced residual stress generation employed in the current study minimises the complexities of RCF modelling and can be applied to any through-hardened bearing steel.