A Multiscale Overview of Modelling Rolling Cyclic Fatigue in Bearing Elements

During service, bearing components experience rolling cyclic fatigue (RCF), resulting in subsurface plasticity and decay of the parent microstructure. The accumulation of micro strains spans billions of rolling cycles, resulting in the continuous evolution of the bearing steel microstructure. The bearing steel composition, non-metallic inclusions, continuously evolving residual stresses, and substantial work hardening, followed by subsurface softening, create further complications in modelling bearing steel at different length scales. The current study presents a multiscale overview of modelling RCF in terms of plastic deformation and the corresponding microstructural alterations. This article investigates previous models to predict microstructural alterations and material hardening approaches widely adopted to mimic the cyclic hardening response of the evolved bearing steel microstructure. This review presents state-of-the-art, relevant reviews in terms of this subject and provides a robust academic critique to enhance the understanding of the elastoplastic response of bearing steel under non-proportional loadings, damage evolution, and the formation mechanics of microstructural alterations, leading to the increased fatigue life of bearing components. It is suggested that a multidisciplinary approach at various length scales is required to fully understand the micromechanical and metallurgical response of bearing steels widely used in industry. This review will make significant contributions to novel design methodologies and improved product design specifications to deliver the durability and reliability of bearing elements.


Background
Rolling bearing elements (RBEs) are fundamental parts of the automotive industry, turbines, and precision machinery. Bearings are an integral part of different mechanisms and assemblies. RBEs are used to transfer loads and motion while providing structural support to the machine components [1]. A bearing element consists of an outer part and an inner part which are connected to the machine parts and several rolling elements. Bearing elements are categorised depending upon the roller type (i.e., cylindrical roller, spherical roller, tapered roller) and can be used based on load and lubrication requirements. During service, rolling bearing elements are subjected to highly varying loads in high-speed turbines, precision machinery, automotive, and aerospace industries. The rolling bearing elements are operated under non-proportional loadings over hundreds of millions of stress cycles. The bearing material response under RCF is dependent on many factors, i.e., contact stress, frictional coefficient, residual stresses, carbide volume fraction, inclusions, and bearing material microstructure [2]. During operation, microstructure, microtexture, compressive residual stresses, and localised hardness in the subsurface area evolve with progressing stress cycles [3,4], making the RCF phenomenon much more complex and heterogeneous. Moreover, non-metallic inclusions, surface defects, slippage, and starved lubrication conditions create further scatter and uncertainty in the service life of a rolling element [5,6].
Rolling contact fatigue (RCF) in bearings is manifested as subsurface plastic damage accumulated under the rolling contact path. This subsurface damage is highly localised and is accumulated at the microstructural scale [7] (typically in the vicinity of carbides due to stress concentration). The ability of a material to strain-harden with the continuous accumulation of plasticity is very high where the material yield point translates (kinematic hardening) and expands (isotropic hardening). Therefore, the continuous evaluation of material hardening/softening becomes critical for researchers to acquire the cyclic hardening response of evolved bearing steel, which is quite difficult to predict with simplistic material hardening models. It is well known that this subsurface plastic deformation involves a forward flow of material parallel to rolling direction, which causes compressive residual stresses in the axial and circumferential direction [8,9]. The deep zone residual stresses arise due to the decomposition of parent microstructure and accumulation of non-uniform plasticity [10] during RCF. Earlier studies [3,11,12] reported that the life span of bearing elements during operation can be altered due to deep zone compressive residual stresses (CRS). The CRS help in the closure of crack propagation by impeding the stress concentration at the notch [12,13]. It can also affect the amplitude of applied stress amplitude. During microplasticity, the microstructure of the subsurface region evolves, and the material phase transformation takes place, resulting in formation of microstructural features, i.e., dark etching regions (DERs) and white etching bands (WEBs). Dark etching regions (DERs) are observed as black patches under an optical microscope after etching with Nital (3% nitric acid solution in ethanol) solution. DERs are believed to be a tempered form of martensite [14] and are formed much earlier than white etching bands (WEBs). WEBs, in contrast to DERs, are developed in the later life of RCF and appeared as bright contrast under an optical microscope. The WEBs are considered ferritic in nature and remain throughout the RCF cycles. These elongated white bands are reported as rectangular disc-shaped, which are inclined at an angle of approximately 30 • and 80 • with the contact path and are termed as lower angle bands (LABs) and higher angle bands (HABs), respectively [15][16][17]. The scanning electron microscopy (SEM)-along with electron backscattering diffraction (EBSD)-investigations revealed that WEBs are shear bands which are formed at a depth of principal and orthogonal shear maxima [18,19]. The redistribution of the carbon atoms during RCF was reported as the significant factor to drive this phase transformation phenomenon. Despite the advancements in the material characterisation technique, resulting in new aspects of the localised RCF damage, researchers are still unable to comprehend the formation mechanism of such microstructural features [20,21]. The stress component responsible for DERs formation, along with the unique directionalities of WEBs, have been up for debate. Moreover, it has been reported [22,23] that close proximity exists between CRS and depth distribution of microstructural features (i.e., DERs, WEBs); however, how the microstructural alterations are affected by the presence of compressive stresses has not been quantified.

Bearing Steel
Rolling contact fatigue (RCF) is quite distinct from conventional uniaxial fatigue (i.e., low and high cyclic fatigue), where a bulk of specimen material is tested to generate the stress 'S' versus cycles to failure 'N' (known as S-N curves) [20]. These two regimes of fatigue are defined by the number of cycles to failure. The low cycle fatigue (LCF) usually operates at higher loads, where the plasticity is accumulated homogenously, appears widespread, and is characterised by repeated plasticity (plastic deformation in each cycle). The high cycle fatigue (HCF) is characterised by elastic deformations, where the stress levels are below the macroscopic yielding, and plasticity is accumulated in the isolated regions under loading (i.e., near non-metallic inclusion due to stress concentrations). There is no fixed transition line where the two regimes can be separated, as this transition depends upon material ductility [24]. RCF, in contrast, occurs as a highly localized plastic damage in bearing elements, where the rated life is typically expressed in billions of rolling cycles. This plastic damage is accumulated under the rolling contact path and can be investigated in the axial and circumferential section of bearing components, as illustrated in Figure 1. damage in bearing elements, where the rated life is typically expressed in billions of rolling cycles. This plastic damage is accumulated under the rolling contact path and can be investigated in the axial and circumferential section of bearing components, as illustrated in Figure 1. The standard process to manufacture martensitic steel is the partial austenitization of steel, which consists of austenite and almost 4 vol% cementite (θ carbides). The austenite phase is followed by the quenching process, where almost 10-11% of the austenite phase remains in the martensitic matrix due to high carbon content and is termed as retained austenite (RA) [25]. During rapid cooling, the carbon cannot diffuse out of the solid solution and forms a supersaturated martensitic matrix in a body-centred tetragonal (BCT) structure. It is reported [26] that the quenched martensitic solid solution contains internal stress due to material phase transformations, along with resulting volume expansion/contraction, and is compensated by a large number of dislocations within the matrix. The resulting martensitic structure can attain hardness as high as 900 HV, but is highly brittle and contains internal residual stresses. The brittleness of the microstructure can be overcome by tempering the steel at 200 °C for 1-2 h to reduce the hardness and introduce ductility. During tempering, carbon leaves the solid solution and creates nanosized precipitates (ε and η), also called tempered carbides, or transition carbides. The size and morphology of these transition carbides depend upon the tempering condition. The proper balance of the primary carbides (θ) and transition carbides (ε and η) in the microstructure ensures the bearing material toughness and strength [27]. Transmission electron microscopy (TEM) of the through-hardened bearing steel has revealed that the θ carbides have (Fe, Cr)3C structure and are made up of 12 wt.% chromium [28]. The addition of chromium (Cr) to the bearing steel increases the wear resistance, strength, and surface hardness. The overall microstructure of standard martensitic steel contains needle-like martensite plates, 10-11% RA, along with undissolved carbides (θ, ε, and η), which can be observed in Figure 2. The standard process to manufacture martensitic steel is the partial austenitization of steel, which consists of austenite and almost 4 vol% cementite (θ carbides). The austenite phase is followed by the quenching process, where almost 10-11% of the austenite phase remains in the martensitic matrix due to high carbon content and is termed as retained austenite (RA) [25]. During rapid cooling, the carbon cannot diffuse out of the solid solution and forms a supersaturated martensitic matrix in a body-centred tetragonal (BCT) structure. It is reported [26] that the quenched martensitic solid solution contains internal stress due to material phase transformations, along with resulting volume expansion/contraction, and is compensated by a large number of dislocations within the matrix. The resulting martensitic structure can attain hardness as high as 900 HV, but is highly brittle and contains internal residual stresses. The brittleness of the microstructure can be overcome by tempering the steel at 200 • C for 1-2 h to reduce the hardness and introduce ductility. During tempering, carbon leaves the solid solution and creates nanosized precipitates (ε and η), also called tempered carbides, or transition carbides. The size and morphology of these transition carbides depend upon the tempering condition. The proper balance of the primary carbides (θ) and transition carbides (ε and η) in the microstructure ensures the bearing material toughness and strength [27]. Transmission electron microscopy (TEM) of the through-hardened bearing steel has revealed that the θ carbides have (Fe, Cr)3C structure and are made up of 12 wt.% chromium [28]. The addition of chromium (Cr) to the bearing steel increases the wear resistance, strength, and surface hardness. The overall microstructure of standard martensitic steel contains needle-like martensite plates, 10-11% RA, along with undissolved carbides (θ, ε, and η), which can be observed in Figure 2. The nominal stress of a bearing element during service is usually below 2.5 GPa; consequently, the applied stress is unable to cause material flow. However, at a microstructural scale, the stress amplifies due to stress intensity factor around carbides and nonmetallic inclusions causing an irreversible plastic flow in the subsurface region of bearing material, which is accumulated with each rolling cycle [30,31]. It is reported [32] that the  The nominal stress of a bearing element during service is usually below 2.5 GPa; consequently, the applied stress is unable to cause material flow. However, at a microstructural scale, the stress amplifies due to stress intensity factor around carbides and non-metallic inclusions causing an irreversible plastic flow in the subsurface region of bearing material, which is accumulated with each rolling cycle [30,31]. It is reported [32] that the RCF-induced cyclic hardening exponent is different from the monotonic strain hardening exponent, which greatly affects the fatigue life of rolling elements. During RCF cyclic loadings, there are four types of material plasticity zones: (i) the elastic region, (ii) the elastic shakedown, (iii) the plastic shakedown, and (iv) ratcheting, as demonstrated in Figure 3. Initially, when the material is loaded elastically, it recovers back to its initial stage and is described in the elastic region. When the load exceeds the yield limit of the material, the material experiences plasticity and resulting residual stress, which resides in and is added to the subsequent loading cycles. The generation of residual stresses is manifested due to two reasons: cyclic plasticity and the transformation of RA in the parent martensitic matrix [33]. If the following cycles do not create any plasticity due to yield point offset, then the material is said to be in the elastic shakedown region. However, if the subsequent loading cycles create combined stress well above the elastic shakedown limit, then the material can experience two possible plasticity modes, e.g., plastic shakedown and ratchetting. In plastic shakedown, the material undergoes closed-loop stress-strain hysteresis, whereas in ratcheting, the material undergoes continuous accumulation of plasticity in an open hysteresis loop [34]. Generally, the RCF life is categorised into three stages: stage I as elastic/plastic shakedown, stage II, where the bearing undergoes steady-state response (which spans to billions of RCF cycles), and stage III as instability or the failure stage. The microstructural alterations, as discussed previously, are dominant in stage II of RCF.  Bearing material under cyclic loadings manifests considerable strain hardening of the subsurface region due to plasticity, which continues for billions of cycles until failure [36]. The prolonged accumulation of subsurface plasticity results in the considerable hardening of bearing material [37,38] and a change in subsurface microstructure [3,39]. During the early shakedown stage of RCF (stage I), the continuously evolving cyclic plasticity results in the build-up of residual stress, while the material yield strength increases due to significant work hardening, as explained earlier. It is noted [40] that a higher load applied at the initial stage results in a significant work hardening, which causes an extended fatigue life of bearing elements. Experimental findings [41] from subsurface hardness measure- Bearing material under cyclic loadings manifests considerable strain hardening of the subsurface region due to plasticity, which continues for billions of cycles until failure [36]. The prolonged accumulation of subsurface plasticity results in the considerable hardening of bearing material [37,38] and a change in subsurface microstructure [3,39]. During the early shakedown stage of RCF (stage I), the continuously evolving cyclic plasticity results in the build-up of residual stress, while the material yield strength increases due to significant work hardening, as explained earlier. It is noted [40] that a higher load applied at the initial stage results in a significant work hardening, which causes an extended fatigue life of bearing elements. Experimental findings [41] from subsurface hardness measurements after progressive RCF cycles have revealed that the subsurface hardening of bearing elements can persist for a longer period, up to 240 million cycles.
Besides considerable work hardening, subsequent softening is also reported [35,[42][43][44] in the later stage of RCF loadings (i.e., stage II). As mentioned previously, bearing elements under cyclic loadings experience initial plastic shakedown in the early stage of the RCF loading cycles. The shakedown limit of an elastic-plastic material is dependent on the shear cyclic yield strength 'k' of bearing material [8,25]. Generally, the shakedown is described in the form of P max /k, where P max is the maximum Hertzian contact pressure. For frictionless line contact and point contact, yielding will occur if this ratio becomes equal to 4.0 and 4.7, respectively [45]. The shear cyclic yield strength 'k' can be estimated as k = where σ y is the tensile cyclic yield strength [45,46]. The cyclic strains which are developed after the shakedown stage are dependent on cyclic yield strength and are typically found in the range of 0 ≤ ε p ≤ 0.002. The inevitable development of microstructural alterations under contact track (driven by the diffusion of carbon in the parent structure [47,48]) and the local temperature peaks during cyclic loadings can trigger potential slip systems, resulting in material softening [20,36,49]. Moreover, for bearing elements under elastic shakedown conditions, the overall material response under loading is generally considered elastic. However, under plastic shakedown conditions, the stress-strain response of the bearing material forms a closed loop, as shown in Figure 3, where the local plastic deformations accumulate near carbides and non-metallic inclusions, causing microstructural alterations and cyclic softening, which results in further increase in the cyclic strain amplitude, leading to failure or instability (stage III) [2]. The bearing material hardening and softening response when subjected to RCF is demonstrated in Figure 4. softening, which results in further increase in the cyclic strain amplitude, leading to failure or instability (stage III) [2]. The bearing material hardening and softening response when subjected to RCF is demonstrated in Figure 4.

Material Modelling
When two rollers roll together with a sufficiently high load, the subsurface region experiences plasticity, which is accumulated under the contact track and the surface laye of the cylinder, displacing in the forward direction of motion [8]. This phenomenon wa initially described in 1963 [50], when it was stated that the forward flow of the surface

Material Modelling
When two rollers roll together with a sufficiently high load, the subsurface region experiences plasticity, which is accumulated under the contact track and the surface layer of the cylinder, displacing in the forward direction of motion [8]. This phenomenon was initially described in 1963 [50], when it was stated that the forward flow of the surface layer of a cylinder is a result of a complex stress state, as described in Figure 5a,b. During the rolling process, compressive residual stresses are evolved in the initial first few loading cycles, provided that the load exceeds the yield limit of the material. Bhargava et al. [51] were among the early researchers to present a numerical model for RCF using a 2D finite element analysis of elastic-plastic behaviour. The centreline hardness is shown as a function of depth. Hardening takes place at a depth of nearly 450 µm. A softened region can also be observed at 500-700 µm depth. Maximum hardening can be observed after 246 million RCF cycles [35].

Material Modelling
When two rollers roll together with a sufficiently high load, the subsurface region experiences plasticity, which is accumulated under the contact track and the surface layer of the cylinder, displacing in the forward direction of motion [8]. This phenomenon was initially described in 1963 [50], when it was stated that the forward flow of the surface layer of a cylinder is a result of a complex stress state, as described in Figure 5a,b. During the rolling process, compressive residual stresses are evolved in the initial first few loading cycles, provided that the load exceeds the yield limit of the material. Bhargava et al. [51] were among the early researchers to present a numerical model for RCF using a 2D finite element analysis of elastic-plastic behaviour.  The initial theoretical and numerical analyses [8,50,52,53] highlighted the linear elasticplastic behaviour of bearing material. Later on, Howell and Bower et al. [54,55] described the isotropic and nonlinear kinematic hardening of bearing material that denotes the linear expansion and translation of the yield surface in the stress space. A three-dimensional semianalytical code was presented by Jacq et al. [56], which was later extended by Chaise and Chen [57,58] in order to estimate the residual stress/strain after a few rolling cycles. They reported that the maximum contact stress and the equivalent plastic strains (PEEQ) were dependent on the ellipticity ratio of the interacting surfaces. The change in composition and heat treatment processes employed in the widely used bearing steels was also characterised. From cyclic push-pull experiments on the three variants of bearing steel, it was reported [46] that a standard martensitic hardened bearing steel represents pronounced hardening when subjected to cyclic loadings, when compared with the bainitic variant of steel. This is due to the fact that RA decays at a much lower stress level, leading to the onset of plastic deformation and consequently, cyclic hardening.
A bilinear elastic-linear kinematic plastic (ELKP) model was presented [45,59] with the help of multivalued cyclic stress-strain relationships. The material parameters utilised in the ELKP material model are elastic modulus 'E' , yield strength 'S y ', kinematic yield strength 'σ k ', and plastic modulus 'M', which are nominally derived from extensive experimental data of fatigue tests. It is believed that the idealisation of the ELKP material model from torsional testing data is feasible for finite element analysis of rolling contact fatigue and is utilised by different researchers [10,[59][60][61][62] for RCF modelling. Figure 6a shows a generic representation of the stress-strain curve for ELKP relationships whose simulation parameters are acquired from the torsional testing of standard AISI 52100 bearing steel. The resulting shear cyclic stress-strain relationship depicts a close hysteresis loop because of rolling contact loading and is represented in Figure 6b. The stabilised shear stress-strain loop (indicated in dark colour) shows the steady state response indicative of plastic shake- down. Hahn et al. [43] presented a 2D finite element model of the ELKP response, and the results were compared with the available experimental data on deep groove ball bearings (DGBB) [3,14]. The finite element simulation results, based on the torsional fatigue testing of standard bearing steel, revealed that the cyclic plasticity and modest residual stresses corroborate well with the early life of RCF (i.e., less than 10 6 cycles). However, for the later life of RCF (i.e., 10 8 ≤ N ≤ 10 10 cycles), the simulation results underpredicted the formation of residual stress, in which the bearing material experienced significant metallurgical changes under RCF [43]. It was stated that the later formation of residual stresses, accompanied by metallurgical change, depends upon the relative volume fraction of the decayed microstructure. However, no direct relationship of compressive residual stresses with the subsurface microstructural alterations was considered in detail. The underprediction of calculated residual stresses from finite element simulation leads to the argument that a more complex material model is required to compensate for the mechanical-metallurgical response associated with the RCF. Later, it was also reported [63] that the bearing material subjected to rolling cycles exhibits nonlinear kinematic hardening, which can be described by the Ramberg-Osgood model. A later study [10] presented J2 plasticity-based finite element simulations incorporating linear kinematic hardening and nonlinear kinematic hardening, consecutively. Nevertheless, the cyclic behaviour of a hardened bearing under RCF is complex and challenging to define with simplistic isotropic hardening [64] and kinematic hardening [65][66][67]. To incorporate the complex hardening of bearing steel, a combined nonlinear i tropic and kinematic hardening (NIKH) material model was employed to model subs face plasticity during cyclic loadings [68]. The NIKH material model was originally p sented by Chaboche [69] and is capable of simulating various material responses un multiaxial cyclic loading, i.e., the Bauschinger effect, with mean stress relaxation, pla shakedown, and ratcheting. Due to its complex nature, the NIKH model has been wid implemented in fatigue applications [70][71][72]. It was reported [68,73] that the NIHK ma rial model can accurately mimic the cyclic hardening and continuous accumulation plasticity (ratcheting) in the vicinity of carbides. However, the estimation of material c stants for such a complex model is very crucial and requires significant experimental d for calibration. As discussed earlier, such material data can be acquired from conventio cyclic stress-strain curve or torsional fatigue tests of bearing steel [74].
Pandkar et al. [68,73] presented a semi-empirical approach based on combined n linear isotropic/kinematic hardening to model the ratcheting response of case harden bearing steel and to estimate the continuous accumulation of plastic strains in the vicin of carbides. The cyclic hardening parameters for a case depth of M50-NiL bearing st were estimated from the subsurface micro-hardness measurement. The NIKH mate model employed in this study is based on von Mises yield criterion, with an associat flow rule and an additive decomposition of strain tensor. The generalised associative fl rule for the NIKH material model is given in Equation (1) as [68], f S, q, ε | S q | σ To incorporate the complex hardening of bearing steel, a combined nonlinear isotropic and kinematic hardening (NIKH) material model was employed to model subsurface plasticity during cyclic loadings [68]. The NIKH material model was originally presented by Chaboche [69] and is capable of simulating various material responses under multiaxial cyclic loading, i.e., the Bauschinger effect, with mean stress relaxation, plastic shakedown, and ratcheting. Due to its complex nature, the NIKH model has been widely implemented in fatigue applications [70][71][72]. It was reported [68,73] that the NIHK material model can accurately mimic the cyclic hardening and continuous accumulation of plasticity (ratcheting) in the vicinity of carbides. However, the estimation of material constants for such a complex model is very crucial and requires significant experimental data for calibration. As discussed earlier, such material data can be acquired from conventional cyclic stress-strain curve or torsional fatigue tests of bearing steel [74].
Pandkar et al. [68,73] presented a semi-empirical approach based on combined nonlinear isotropic/kinematic hardening to model the ratcheting response of case hardened bearing steel and to estimate the continuous accumulation of plastic strains in the vicinity of carbides. The cyclic hardening parameters for a case depth of M50-NiL bearing steel were estimated from the subsurface micro-hardness measurement. The NIKH material model employed in this study is based on von Mises yield criterion, with an associative flow rule and an additive decomposition of strain tensor. The generalised associative flow rule for the NIKH material model is given in Equation (1) as [68], where 'S' is the stress vector, 'q' is back-stress tensor, 'ε − pl ' is the accumulated plastic strain, and 'σ y ' is the instantaneous yield strength. The depth distribution of the change in hardness of the evolved microstructure of the case-hardened bearing steel was integrated into the NIKH material model. Figure 7 shows the ratcheting response of bearing steel subjected to uniaxial loading in a stress-controlled environment. It can be seen that equivalent or von Mises stresses indicate substantial work hardening, with continuous accumulation of PEEQ. Simulation results showed that the cyclic hardening via ratcheting is promoted at the scale of the carbide microstructure. The hard spheroidised carbide particles act as local stress increasers and are the primary driver for the localised subsurface hardening. A similar approach was extended in a more recent research work [75] in order to simulate the deep zone residual stresses in an evolved bearing steel microstructure. The highly localised subsurface hardness change was incorporated into the Tabor rule [76] to convert the nanoindentation hardness to equivalent flow stresses. The evolved flow stresses were then fed into the newly developed 3D finite model in order to mimic the cyclic hardening response of the evolved bearing steel microstructure.  A continuum damage mechanics (CDM) model was presented by Morris et al. [23], in which microstructural deterioration, material degradation, phase transformation, and resulting residual stress formation were presented. As discussed previously, the development of residual stresses during RCF is associated with the microplasticity and the decay of RA [1,22]. In a 2D rolling fatigue simulation based on CDM [12], Shen et al. presented the effects of RA and compressive residual stresses on the fatigue life of carburised AISI 8620 steel. It was reported that the presence of compressive residual stresses is beneficial for fatigue life by reducing the damage evolution rate. Moreover, the higher amounts of RA resulted in the increased RCF life of carburised steel. In addition, the initial presence of higher RA also leads to prolonged subsurface-initiated spalling life [77]. Morris et al. [23] presented a novel modelling approach to estimate residual stresses as a function of the decay of RA. A CDM model was integrated into a finite element analysis to capture A continuum damage mechanics (CDM) model was presented by Morris et al. [23], in which microstructural deterioration, material degradation, phase transformation, and resulting residual stress formation were presented. As discussed previously, the development of residual stresses during RCF is associated with the microplasticity and the decay of RA [1,22]. In a 2D rolling fatigue simulation based on CDM [12], Shen  8620 steel. It was reported that the presence of compressive residual stresses is beneficial for fatigue life by reducing the damage evolution rate. Moreover, the higher amounts of RA resulted in the increased RCF life of carburised steel. In addition, the initial presence of higher RA also leads to prolonged subsurface-initiated spalling life [77]. Morris et al. [23] presented a novel modelling approach to estimate residual stresses as a function of the decay of RA. A CDM model was integrated into a finite element analysis to capture the microstructural deterioration, phase transformation, and residual stress formation. The model was able to predict the residual stress generation under RCF; however, it neglects the plasticity-induced residual stresses during cyclic loadings.
To model the subsurface plasticity and corresponding stress state to initiate the spalling in rolling contact fatigue, a finite element model was integrated with continuum damage mechanics (CDM) [62]. It was reported [13,43,78] that the plastic deformation of the subsurface region during RCF eventually leads to the nucleation of subsurface cracks. A similar approach was utilized by Walvekar [79] in a three-dimensional finite element analysis to report 3D spall formation and the fatigue life scatter of bearing components. A randomly generated Voronoi tessellation technique [80] was implemented to replicate the intergranular paths of the microstructure. A similar meshing technique was implemented in substantial fatigue models [60,62,79,[81][82][83][84][85] to imitate the characteristic grain boundary interactions. Figure 8a represents a three-dimensional topology of the circular rolling contact, modelled using the Voronoi tessellation. The evolution of subsurface damage, represented in Figure 8b, indicates subsurface crack initiation and then propagation to form a spall, in conjunction with the continuous accumulation of plastic strains. Figure 8c shows the Weibull probability plot of RCF failure, as acquired from simulation and experimental results. It was testified that the fatigue damage is induced due to the PEEQ around broken intergranular joints leading to crack propagation and eventually, RCF failure. in substantial fatigue models [60,62,79,[81][82][83][84][85] to imitate the characteristic grain boundary interactions. Figure 8a represents a three-dimensional topology of the circular rolling contact, modelled using the Voronoi tessellation. The evolution of subsurface damage, represented in Figure 8b, indicates subsurface crack initiation and then propagation to form a spall, in conjunction with the continuous accumulation of plastic strains. Figure 8c shows the Weibull probability plot of RCF failure, as acquired from simulation and experimental results. It was testified that the fatigue damage is induced due to the PEEQ around broken intergranular joints leading to crack propagation and eventually, RCF failure. The cyclic hardening and fatigue life models are based on the constitutive properties of the bearing steel obtained from either the monotonic stress-strain response or the conventional torsional fatigue test and are therefore unable to present a true picture of the multiaxial loading state during RCF [20,29]. Moreover, the cyclic lives of torsional tested steel are seven to nine orders less than the actual life cycle of well-lubricated bearing steel operated under rolling contact conditions [45]. Future work is recommended to incorporate a complex cyclic hardening response with the gradual evolution of bearing steel microstructure. Bridging the gap between the mechanical and metallurgical aspects in modelling RCF would be highly beneficial in estimating residual stresses and the evolved subsurface response.

Microstructural Modelling
Preliminary investigations [86,87] suggested that microstructural alterations, i.e., DERs and WEBs, are developed due to tempering and localised temperature increase The cyclic hardening and fatigue life models are based on the constitutive properties of the bearing steel obtained from either the monotonic stress-strain response or the conventional torsional fatigue test and are therefore unable to present a true picture of the multiaxial loading state during RCF [20,29]. Moreover, the cyclic lives of torsional tested steel are seven to nine orders less than the actual life cycle of well-lubricated bearing steel operated under rolling contact conditions [45]. Future work is recommended to incorporate a complex cyclic hardening response with the gradual evolution of bearing steel microstructure. Bridging the gap between the mechanical and metallurgical aspects in modelling RCF would be highly beneficial in estimating residual stresses and the evolved subsurface response.

Microstructural Modelling
Preliminary investigations [86,87] suggested that microstructural alterations, i.e., DERs and WEBs, are developed due to tempering and localised temperature increase during cyclic loadings. However, no experimental evidence was found to support such an argument. Nevertheless, the threshold stress for such microstructural alterations represents this mechanism to be stress-induced, regardless of bearing running time. The microstructural phase transformation of the martensite (parent material) into ferrite in the affected zones of RCF can be explained by the carbon migration and diffusion theories. Mitamura et al. [88] conducted RCF studies at tempering conditions of bearing steel and reported that the activation energies of the fatigue process (from Arrhenius plots) correspond to those of carbon diffusion in body-centred-cubic (bcc) ferrite, as shown in Figure 9. The experimental and numerical results indicated that the microstructural alterations during RCF are controlled by the diffusion of carbon atoms in the ferrite matrix. Solid state diffusion is classified as vacancy diffusion and interstitial diffusion, where the latter is relevant to structural changes during rolling contact fatigue [89]. During interstitial diffusion, carbon atoms diffuse from one interstitial site to a neighbouring site, without permanently displacing the martensite matrix. during cyclic loadings. However, no experimental evidence was found to support such an argument. Nevertheless, the threshold stress for such microstructural alterations represents this mechanism to be stress-induced, regardless of bearing running time. The microstructural phase transformation of the martensite (parent material) into ferrite in the affected zones of RCF can be explained by the carbon migration and diffusion theories. Mitamura et al. [88] conducted RCF studies at tempering conditions of bearing steel and reported that the activation energies of the fatigue process (from Arrhenius plots) correspond to those of carbon diffusion in body-centred-cubic (bcc) ferrite, as shown in Figure  9. The experimental and numerical results indicated that the microstructural alterations during RCF are controlled by the diffusion of carbon atoms in the ferrite matrix. Solid state diffusion is classified as vacancy diffusion and interstitial diffusion, where the latter is relevant to structural changes during rolling contact fatigue [89]. During interstitial diffusion, carbon atoms diffuse from one interstitial site to a neighbouring site, without permanently displacing the martensite matrix. The development of WEBs was initially defined as shear bands by K. L. Johnson [90]. It was claimed that the white bands originate at a depth of maximum Hertz shear stress and initiate on planes of maximum shear stress. Polonsky and Keer [15], in a detailed critical analysis, reconstructed the theoretical and analytical models presented by previous researchers [90][91][92] to describe the formation mechanism of white etching bands. It was suggested that the diffusional outflow of carbon from WEBs to form lenticular carbides and ferrite plays a crucial role in WEBs formation. It was reported that the unique orientations of LABs (30°) and HABs (80°) arise from the maximisation of relative normal stress (deviatoric stress) and relative shear stress across the WEBs, instead of from a static state The development of WEBs was initially defined as shear bands by K. L. Johnson [90]. It was claimed that the white bands originate at a depth of maximum Hertz shear stress and initiate on planes of maximum shear stress. Polonsky and Keer [15], in a detailed critical analysis, reconstructed the theoretical and analytical models presented by previous researchers [90][91][92] to describe the formation mechanism of white etching bands. It was suggested that the diffusional outflow of carbon from WEBs to form lenticular carbides and ferrite plays a crucial role in WEBs formation. It was reported that the unique orientations of LABs (30 • ) and HABs (80 • ) arise from the maximisation of relative normal stress (deviatoric stress) and relative shear stress across the WEBs, instead of from a static state of stress, as suggested previously by Zwirlein et al. [92,93]. A similar theoretical model was presented by Abdullah et al. [94] in recent research work to explain the unique directionalities and formation mechanisms of white bands. Figure 10 represents the distribution of relative normal stress and relative shear stress at various orientations of WEBs. It can be seen that the LABs tend to adopt such orientation, where the relative normal stress maximum is last among all major extrema. To develop predictive models for microstructural alterations (DERs/WEBs), a realistic mechanism for the dissolution and redistribution of carbon atoms in a martensite/ferrite matrix is an important aspect. It was suggested by [95] that the development of white bands is a phenomenon of recrystallization, and it involves carbon diffusion resulting in the development of carbon depleted discs (ferrites). These ferrite regions/bands are immune to etching and appear as white contrast bands. Mirza et al. [96] linked thermomechanical and microstructural development with recrystallization and textural evolution during the rolling passes. The carbon from the martensitic matrix becomes trapped at defects (dislocation) or precipitates as cementite bands, creating lenticular carbides. The theory of carbon diffusion during the martensite-ferrite phase transition was supported by previous research work [39,87,97]. Following that, an austenite-ferrite phase transformation model was reported by Pernach et al. [98] which employed a finite difference scheme to solve Fick's diffusion equation. The model successfully described ferrite volume fraction, carbon segregation and kinetics of phase transformation, but it was based on the pre-existing carbon concentration gradients and hence, was unable to state the true picture of RCF. Zhao et al. [99] presented a sequentially coupled deformation-diffusion analysis to model oxygen diffusion along grain boundaries in a nickel-based superalloy. A crystal plasticity model was used to define the material constitutive behaviour, and the diffusion of oxygen was calculated in the presence of a stress tensor. This stress-assisted diffusion model was further extended by Warhadpande et al. [89] by coupling the stressassisted diffusion flux with J2 plasticity. It was described that the dissipated plastic strain energy during cyclic plasticity drives the carbon outflow from the parent matrix. Figure  11a shows the distribution of dissipated plastic energy per unit volume as the load moves from left to right in a 2D domain. Correspondingly, the distribution of carbon and unique To develop predictive models for microstructural alterations (DERs/WEBs), a realistic mechanism for the dissolution and redistribution of carbon atoms in a martensite/ferrite matrix is an important aspect. It was suggested by [95] that the development of white bands is a phenomenon of recrystallization, and it involves carbon diffusion resulting in the development of carbon depleted discs (ferrites). These ferrite regions/bands are immune to etching and appear as white contrast bands. Mirza et al. [96] linked thermomechanical and microstructural development with recrystallization and textural evolution during the rolling passes. The carbon from the martensitic matrix becomes trapped at defects (dislocation) or precipitates as cementite bands, creating lenticular carbides. The theory of carbon diffusion during the martensite-ferrite phase transition was supported by previous research work [39,87,97]. Following that, an austenite-ferrite phase transformation model was reported by Pernach et al. [98] which employed a finite difference scheme to solve Fick's diffusion equation. The model successfully described ferrite volume fraction, carbon segregation and kinetics of phase transformation, but it was based on the pre-existing carbon concentration gradients and hence, was unable to state the true picture of RCF. Zhao et al. [99] presented a sequentially coupled deformation-diffusion analysis to model oxygen diffusion along grain boundaries in a nickel-based superalloy. A crystal plasticity model was used to define the material constitutive behaviour, and the diffusion of oxygen was calculated in the presence of a stress tensor. This stress-assisted diffusion model was further extended by Warhadpande et al. [89] by coupling the stress-assisted diffusion flux with J2 plasticity. It was described that the dissipated plastic strain energy during cyclic plasticity drives the carbon outflow from the parent matrix. Figure 11a shows the distribution of dissipated plastic energy per unit volume as the load moves from left to right in a 2D domain. Correspondingly, the distribution of carbon and unique orientations of WEBs (i.e., LABs and HABs) can be seen in Figure 11b in terms of carbon concentration. It is suggested that during cyclic rolling loadings, carbon migrates in specific orientations in the material domain, corresponding to the directions of WEB. The unique orientation of carbon outflow from the parent matrix represents the formation of 30 • and 80 • WEBs. The model successfully presents the formation of WEBs orientation; however, it relies on the assumption that material elastic-plastic properties remain homogenous, and that the diffusion coefficient is not affected by temperature or stress variations. however, it relies on the assumption that material elastic-plastic properties remain homogenous, and that the diffusion coefficient is not affected by temperature or stress variations. During RCF loading, the dislocation density can also play a role to drive carbon. As the ferrite phase arises due to the localised plastic deformation of martensite, it contains a significant dislocation density, as reported by Fu et al. [100], and is considered a Cottrell atmosphere [101]. The higher dislocation density can favourably drag carbon atoms and form a strong interaction between nearby dislocation and carbon atoms. During bearing operation, a significant amount of dislocation gliding occurs, which can act as a driving force for carbon migration within the solid solution. Kang et al. [47] presented an analytical model of martensite tempering assisted by dislocation glide during rolling contact fatigue. It was reported that during RCF, the carbon glide according to their atmospheres, enabling carbon to be transported from the matrix to nearby carbides, resulting in carbide thickening within the DER ferrite. The transport of carbon atoms is calculated by considering the interaction between the carbon and the dislocations, binding the energy of the carbon with the carbides.
Based on the dislocation gliding mechanism, Fu et al. presented strain-induced martensite decay [102] and stress-induced carbide precipitation [48] models to estimate the initiation and growth of DERs and WEBs, respectively. Later on, a unified model [103] for microstructural alterations was presented, suggesting the carbon migration distance as a major aspect of their formation. For dark etching regions, the carbon travels hundreds of nanometres of distance to precipitate with the pre-existing carbides. This thickening of the carbides within the DER zone was confirmed by atomic probe tomography (APT) investigations [102]. For white etching band formation, the carbon migrates from the ferrite matrix and travels several micrometres to develop lenticular carbides [100]. The overall schematic to model microstructural alterations (DER, WEB) is demonstrated in Figure 12. The generalised equations for carbon flux J equilibria in microstructural alterations are given by [103] d r , l During RCF loading, the dislocation density can also play a role to drive carbon. As the ferrite phase arises due to the localised plastic deformation of martensite, it contains a significant dislocation density, as reported by Fu et al. [100], and is considered a Cottrell atmosphere [101]. The higher dislocation density can favourably drag carbon atoms and form a strong interaction between nearby dislocation and carbon atoms. During bearing operation, a significant amount of dislocation gliding occurs, which can act as a driving force for carbon migration within the solid solution. Kang et al. [47] presented an analytical model of martensite tempering assisted by dislocation glide during rolling contact fatigue. It was reported that during RCF, the carbon glide according to their atmospheres, enabling carbon to be transported from the matrix to nearby carbides, resulting in carbide thickening within the DER ferrite. The transport of carbon atoms is calculated by considering the interaction between the carbon and the dislocations, binding the energy of the carbon with the carbides.
Based on the dislocation gliding mechanism, Fu et al. presented strain-induced martensite decay [102] and stress-induced carbide precipitation [48] models to estimate the initiation and growth of DERs and WEBs, respectively. Later on, a unified model [103] for microstructural alterations was presented, suggesting the carbon migration distance as a major aspect of their formation. For dark etching regions, the carbon travels hundreds of nanometres of distance to precipitate with the pre-existing carbides. This thickening of the carbides within the DER zone was confirmed by atomic probe tomography (APT) investigations [102]. For white etching band formation, the carbon migrates from the ferrite matrix and travels several micrometres to develop lenticular carbides [100]. The overall schematic to model microstructural alterations (DER, WEB) is demonstrated in Figure 12. The generalised equations for carbon flux J d equilibria in microstructural alterations are given by [103]  Equations (2) and (3) can be used to calculate carbon diffusion from zone I to II, where N represents the stress cycles, r represents the thickness of carbide, l represents the thickening of the lenticular carbides, V represents the volume of the entire system, and C denotes carbon concentration per unit volume. Abdullah et al. [94,104] presented extensive experimental data on the formation of DERs and WEBs, along with some suggestions to improve existing analytical models. Figure 13 shows the numerical prediction of the formation of DERs with progressive RCF cycles at (a) 4 GPa and (b) 5 GPa contact stress, along with 40 °C, 100 °C, and 160 °C operating temperatures. It was reported that the dislocation gliding model overestimates the formation of microstructural alterations. The main reason for deviation from the experimental data was suggested to be the fact that at extreme loading conditions (near the tempering conditions of bearing steel), the inclusion of carbon diffusivity becomes imperative. It is suggested that the dislocation gliding model is a simplification of carbon migration within the martensitic matrix. According to Fu's model [102], during DERs formation, carbon migrates from the martensite lattice to the pre-existing tempered state to form θ-Fe3C, whereas the initial carbides remain unchanged. APT results have revealed that three forms of carbides coexist in the DER zone, e.g., θ-Fe3C, η-Fe2C, and ε-Fe2.4C [102]. Moreover, the lenticular carbides formed during WEBs have relatively less carbon content, as compared to θ-Fe3C [105]. Following extensive experimental investigations, a semi-empirical model for predicting WEBs [106] was presented based on the saturation density of LABs and HABs under a range of contact pressures and stress cycles. Based on the growth pattern of WEBs, it was suggested that the WEBs formation is mainly driven by the Equations (2) and (3) can be used to calculate carbon diffusion from zone I to II, where N represents the stress cycles, r p represents the thickness of carbide, l LC represents the thickening of the lenticular carbides, V o represents the volume of the entire system, and C V denotes carbon concentration per unit volume.
Abdullah et al. [94,104] presented extensive experimental data on the formation of DERs and WEBs, along with some suggestions to improve existing analytical models. Figure 13 shows the numerical prediction of the formation of DERs with progressive RCF cycles at (a) 4 GPa and (b) 5 GPa contact stress, along with 40 • C, 100 • C, and 160 • C operating temperatures. It was reported that the dislocation gliding model overestimates the formation of microstructural alterations. The main reason for deviation from the experimental data was suggested to be the fact that at extreme loading conditions (near the tempering conditions of bearing steel), the inclusion of carbon diffusivity becomes imperative. It is suggested that the dislocation gliding model is a simplification of carbon migration within the martensitic matrix. According to Fu's model [102], during DERs formation, carbon migrates from the martensite lattice to the pre-existing tempered state to form θ-Fe 3 C, whereas the initial carbides remain unchanged. APT results have revealed that three forms of carbides coexist in the DER zone, e.g., θ-Fe 3 C, η-Fe 2 C, and ε-Fe 2.4 C [102]. Moreover, the lenticular carbides formed during WEBs have relatively less carbon content, as compared to θ-Fe 3 C [105]. Following extensive experimental investigations, a semiempirical model for predicting WEBs [106] was presented based on the saturation density of LABs and HABs under a range of contact pressures and stress cycles. Based on the growth pattern of WEBs, it was suggested that the WEBs formation is mainly driven by the diffusion process. In a recent mechanistic study [107], it was reported that both LABs and HABs arise due to recrystallisation from energy build-up in the initial microstructure, which later transforms to the elongated ferrite grains via the grain rotation/coalescence recovery mechanism owing to subsurface plasticity. Mustafa et al. [17] suggested that HABs are formed in the densely formed areas of LABs. However, at elevated temperatures and loads, HABs prior to LABs have also been reported [88,94]. The early formation of HABs suggests that the reversal of the sequence of WEBs (e.g., HABs before LABs) does not depend upon applied load, but rather the temperature-load combination. diffusion process. In a recent mechanistic study [107], it was reported that both LABs and HABs arise due to recrystallisation from energy build-up in the initial microstructure, which later transforms to the elongated ferrite grains via the grain rotation/coalescence recovery mechanism owing to subsurface plasticity. Mustafa et al. [17] suggested that HABs are formed in the densely formed areas of LABs. However, at elevated temperatures and loads, HABs prior to LABs have also been reported [88,94]. The early formation of HABs suggests that the reversal of the sequence of WEBs (e.g., HABs before LABs) does not depend upon applied load, but rather the temperature-load combination. It is proposed to couple the existing dislocation gliding models [47,48,102,103] with a more sophisticated elastoplastic model in order to mimic the cyclic hardening response of bearing steel. It has been debated for decades whether the microstructural alterations are governed by the diffusion of carbon or subsurface plastic deformations. It is more likely that the formation mechanism for such structural alterations is affected by carbon diffusion through plastic deformation; however, their individual contributions to the formation mechanism of microstructural alterations is still debatable.

Overview and Conclusions
The complex heterogenous nature of bearing steel can be realised by chemical composition, non-metallic inclusion, cyclic plasticity, strain hardening/softening, residual stresses evolution and altered microstructure. Modelling the cyclic plasticity and predicting the formation of microstructural alterations become imperative for researchers to understand micromechanical and metallurgical response of bearing steel subjected to RCF. Considerable work has been done so far on the damage mechanics of subsurface-initiated failure to estimate the fatigue life of bearing components; however, limited research has It is proposed to couple the existing dislocation gliding models [47,48,102,103] with a more sophisticated elastoplastic model in order to mimic the cyclic hardening response of bearing steel. It has been debated for decades whether the microstructural alterations are governed by the diffusion of carbon or subsurface plastic deformations. It is more likely that the formation mechanism for such structural alterations is affected by carbon diffusion through plastic deformation; however, their individual contributions to the formation mechanism of microstructural alterations is still debatable.

Overview and Conclusions
The complex heterogenous nature of bearing steel can be realised by chemical composition, non-metallic inclusion, cyclic plasticity, strain hardening/softening, residual stresses evolution and altered microstructure. Modelling the cyclic plasticity and predicting the formation of microstructural alterations become imperative for researchers to understand micromechanical and metallurgical response of bearing steel subjected to RCF. Considerable work has been done so far on the damage mechanics of subsurface-initiated failure to estimate the fatigue life of bearing components; however, limited research has been reported on modelling the evolved microstructure response of the RCF affected zone.
The continuous accumulation of cyclic plasticity is modelled with the help of elasticlinear kinematic plastic (ELKP), linear/nonlinear kinematic hardening, or combined nonlinear isotropic/kinematic hardening (NIKH) models. The cyclic hardening coefficients for these material models are acquired and calibrated against extensive experimental data. The constitutive response of bearing steel can be characterised by the subsurface hardness measurements of RCF tested samples, which represent the flow stresses of the RCF-affected region. These evolved flow stresses can be incorporated into material models to mimic the cyclic hardening/softening of the bearing steel microstructure.
A more realistic elastic-plastic finite element model considers the nonlinearity of the bearing steel microstructure, along with cyclic plasticity and degradation of the microstructure. It is recommended to integrate the decay of retained austenite in the NIKH material model to evaluate the progression of residual stresses in the later stages of RCF. This will enable the accurate modelling of the deep zone residual stresses in the evolved stressaffected regions with continuously accumulating plastic strains and microstructural phase transformations. Based on the constitutive parameters obtained from RCF testing, a reliable multifaceted fatigue life prediction model can also be developed, incorporating the effects of parent matrix, carbides, inclusions, residual stresses, material evolution/degradation, and the resulting evolved subsurface stress fields.
The analytical and numerical models of microstructural alterations, reported in the literature, are based on carbon diffusion and/or carbide thickening due to the dislocationassisted carbon migration theory. The carbide thickening models predict the formation of DERs/WEBs in good agreement with the reported experimental data under nominal conditions. However, under extreme conditions, their formation mechanism is greatly influenced by the diffusion process. These analytical models also lack the ability to describe the unique directionalities of the lower angle and higher angle bands. Consequently, a comprehensive understanding is still needed to predict such microstructural alterations. Coupling the three-dimensional elastoplastic models with the dislocation gliding mechanism and integrating the diffusion process is a long-term goal in developing a state-of-the-art prognostic model for microstructural alterations.