The analytical expressions currently available for Hertzian contact stresses are applicable only for homogeneous materials and not for case-hardened bearing steels, which have inhomogeneous microstructure and graded elastic properties in the subsurface region. Therefore, this article attempts to determine subsurface stress fields in ball bearings for graded materials with different ball and raceway geometries in contact. Finite element models were developed to simulate ball-on-raceway elliptical contact and ball-on-plate axisymmetric contact, to study the effects of elastic modulus variation with depth due to case hardening. Ball bearings with low, moderate, and heavy load conditions are considered. The peak contact pressure for case-hardened steel is always more than that of through-hardened steel under identical geometry and loading conditions. Using equivalent contact pressure approach, effective elastic modulus is determined for case-carburized steels, which will enable the use of Hertz equations for different gradations in elastic modulus of raceway material. Nonlinear regression tools are used to predict effective elastic modulus as a weighted sum of surface and core elastic moduli of raceway material and design parameters of ball–raceway contact area. Mesh convergence study and validation of equivalent contact pressure approach are also provided. Implications of subsurface stress variation due to case hardening on bearing fatigue life are discussed.

## Introduction

A number of applications in tribology, geology, optoelectronics, biomechanics, fracture mechanics, and nanotechnology involve components with compositionally graded materials [1]. Transmission components such as rolling element bearings and gears are heat treated using carburization and nitriding processes to enhance resistance to damage and deformation under repeated contact. These surface treatment processes result in gradient in carbide and nitride concentration in the subsurface region, leading to spatial variation of mechanical properties such as elastic modulus, yield strength, and hardness [2]. The mechanical response of materials with spatial gradients in composition and microstructure has been a subject of many experimental and analytical studies by number of researchers [1–6]. Giannakopoulos and Suresh [1,4] analyzed axisymmetric graded half space problem for a point load, and for conical, spherical, and flat indenters. In their study, frictionless contact, constant Poisson's ratio, and Young's modulus variation with depth either as simple power law (*E*(*z*) = *E*_{0}*z ^{k}*

^{′}, 0 ≤

*k*′ < 1) or as exponential variation (

*E*(

*z*) =

*E*

_{0}

*e*) were assumed. Liu et al. [5] solved axisymmetric frictionless contact problem of functionally graded materials using transfer matrix method and Hankel integral transform technique. With the help of numerical techniques, they obtained solutions for contact pressure and contact region for indenters of various geometries. Guler and Erdogan [6] analyzed contact problem of graded metal/ceramic coatings for planar elasticity (plane stress and plane strain) conditions. They considered continuously varying thermo-mechanical properties to analyze influence of material inhomogeneity constant, the coefficient of friction, and various length-scale parameters on the critical stresses that determine fatigue and fracture of coating. All the solutions from previous studies [1,3–6] are applicable only to limited cases such as axisymmetric and planar contact problems. However, in most of the industrial applications, especially for ball bearings, elliptical contact exists between ball and raceway. These conditions require contact analysis of elastically graded materials under three-dimensional (3D) loading. Moreover, previous studies primarily rely on extension of Boussinesq's method, whereas practices in bearing industry have evolved to use of Hertz equations for stress, deformation, and fatigue life calculations [7,8]. Bhattacharyya et al. [9] analyzed three-dimensional point contact using finite element method, but main focus of their study was contact stresses under elastic–plastic loading conditions. Using 3D finite element models and Voronoi tessellations, Weinzapfel et al. [10] proposed framework to account for influence of microstructure topologies on rolling contact fatigue of bearing steels. In their work, they used predefined Hertz stress field; therefore, it is mainly applicable for through-hardened steels. Due to inhomogeneous composition, exact nature of Hertz stress field is unknown for case-hardened steels. Therefore, our goal was to study contact problem of case-hardened ball bearing raceways using Hertz equations under three-dimensional as well as axisymmetric loading conditions.

^{αz}Traditional Hertzian contact solution for determining subsurface elastic stress fields for two bodies in contact is applicable only for homogeneous materials [7]. Bearing *L*_{10} life, defined as number of revolutions for which 90% of bearings survive at a given load, is a strong function of peak Hertzian contact pressure and the peak orthogonal subsurface shear stress [11]. Variations in contact pressure and subsurface stress field due to elastic modulus variation in the case layer of case-hardened steel are expected to introduce a significant correction to bearing life. Therefore, the primary objective of this work is to analyze the influence of elastic modulus gradation on the peak Hertzian contact pressure and subsurface stress-field experienced by the ball–raceway contact. For this, detailed mesh density/convergence study of the 3D finite element models of the ball–raceway contact with elastic modulus gradient is performed. Discussions regarding procedures adopted for the development of these models along with the challenges encountered in finite element method simulations of ball-bearing contacts are also presented.

## Hertz Theory of Contact Mechanics

*Q*is the normal load experienced by two bodies in contact in

*X–Y*plane;

*a*and

*b*represent semimajor and semiminor axes of the elliptical contact area.

*S*will be used to represent maximum compressive stresses experienced by two bodies in contact. Hertz also provided expressions for semimajor axis

*a*and semiminor axis

*b*as [7]

where $\xi $ and *E* are Poisson's ratio and elastic modulus of bodies *I* and *II*, respectively; $\u2211\rho $ is the summation of principal curvatures of two bodies in contact; $a*\u2009and\u2009b*$ are nondimensional parameters defined by curvature difference of two bodies [7].

where the proportionality constant *k* is equal to 0.25, 0.28, and 0.32 for the orthogonal, octahedral, and the maximum shearing stress, respectively, for commonly observed contact conditions (i.e., $(b/a)\u2248$ 0.1064). For roller bearings, these values are 0.25, 0.29, and 0.30, respectively [16]. For ball bearings (point contact), the subsurface depths where maximum orthogonal, octahedral, and subsurface shear stress observed are 0.5*b*, 0.72*b*, 0.76*b*, respectively. Corresponding depths for roller bearings are 0.5*b*, 0.79*b*, 0.79*b*, respectively [16]. Analytical expressions Eqs. (1)–(4) were derived based on the assumption that subsurface properties such as elastic modulus and hardness remain constant [17]. However, for case-hardened bearing steels such as M-50 NiL and P-675, a linear gradation of carbide volume fraction with depth, and consequently a linear variation in elastic modulus with depth [2] have been observed. Due to this gradient, traditional Hertz solution is not applicable; therefore, we resort to a 3D finite element model of the ball–raceway contact with modulus gradient built into the raceway.

## Sensitivity of Bearing Fatigue Life to Elastic Modulus Variations

*L*, is inversely proportional to equivalent dynamic load $Fe$ and also maximum Hertzian contact pressure

*S*as [11]

where *p* and *n* represent load-life and stress-life exponents, respectively. For point contact, stress-life exponent is three times the load-life exponent (*n* = 3*p*) [11]. Based on bearing endurance data available in 1950s, Lundberg and Palmgren [13] evaluated load-life exponent *p* = 3 and therefore *n* = 9. The study by Parker and Zaretsky [18] indicates that for vacuum-processed steels, *p* = 4 and *n* = 3*p* = 12. These results are based on fatigue data collected from five-ball fatigue tester with maximum Hertzian stress in the range of 4.5–6 GPa. However, reevaluation studies by Londhe et al. [19], based on actual bearing fatigue lives reported by bearing and aircraft engine manufacturers [20], indicate that load-life exponent *p* for ball bearings is 4.1 and the corresponding stress-life exponent is then *n* = 3**p* = 12.3. These results are based on Bayesian statistics approach rather than traditional minimization of root-mean-square approach [18,20] used for determining unknown parameters. As seen in Eq. (5), because fatigue life is proportional to the contact pressure raised to large exponent (9, 12, or 12.3), life estimates are sensitive to even small changes in contact pressure. Even subtle reductions in peak contact pressure due to elastic modulus gradient can significantly improve fatigue life prediction of ball bearings.

Analytical study presented in the Appendix shows that peak Hertz contact pressure experienced by ball bearings is fairly sensitive to variations in elastic modulus of raceway material [21]. For example, a 10% reduction in elastic modulus for raceway material results in a 3.6% drop in the peak contact pressure due to reduced contact stiffness. This results in a 38.7% correction in bearing fatigue life prediction with stress-life exponent (*n*) of 9. These numerical results are based on analysis of 209 single row radial deep groove ball bearings (DGBB) with inner and outer raceway diameters of 52.291 mm and 77.706 mm, respectively. Groove radii for both raceways and ball diameter were 6.6 mm and 12.7 mm, respectively. Bearing was loaded radially with 8900 N force on nine rolling elements. Similarly, with stress-life exponents of 12 and 12.3, the expected correction in fatigue life is 54.7% and 56.4%, respectively. Therefore, for accurate prediction of fatigue life of case-carburized bearing steels, Hertz equations must be corrected to account for varying elastic modulus from the surface to core region of the carburized raceway material.

At the microstructural level, case-hardened steels contain carbide precipitates surrounded by steel matrix [2]. These carbides are usually stiffer than steel matrix. Therefore, composite elastic modulus is used to represent equivalent elastic modulus of carbide inclusions and steel matrix at each case depth. For M-50 NiL steel, composite elastic modulus decreases linearly from about 228 GPa to 202 GPa, 11.4% variation over 2 mm case depth [2]. Similar data for P-675 indicates that the composite elastic modulus varies linearly from 224 GPa to 202 GPa, about 9.82% variation, over 1 mm case depth. As discussed earlier, finite element models were constructed to study the dependency of peak Hertzian stress on linearly varying elastic modulus of a case-hardened raceway material. A wide range of variations in elastic moduli and case depths were considered. Finite element models also included different geometries, which resulted in circular and elliptical contact loads. Simulations include steel-on-steel and ceramic-on-steel contact of varying ball sizes. Results from these simulations are presented in following sections.

## Finite Element Analysis

In the present study, all the finite element models were developed using abaqus/standard software with nonlinear contact analysis. These models are shown in Fig. 1. Three-dimensional ball-on-raceway contact, as shown in Fig. 1(a), was used to simulate elliptical contact area observed inside 6309 DGBB. Geometrical dimensions for the model (as per Harris et al. [22]) are specified in Table 1. Simulations included contact between a silicon nitride ball on steel raceway and steel ball on steel raceway using one-quarter symmetry. Symmetric boundary conditions were used along planes perpendicular to elliptical contact area and passing through semimajor and semiminor axis. Similarly, a ball on plate model with a circular contact was also developed using axisymmetric boundary conditions as shown in Fig. 1(b). This model was simulated for silicon nitride and steel balls of 12.7 mm and 25.4 mm diameter on a steel substrate. Details about the geometries used in each simulation are summarized in Table 1. A sample three-dimensional ball-inside-channel model (as per Zaretsky et al. [15]) with infinite bore diameter was also considered in this analysis. As shown in Fig. 1(c), in this model, symmetric boundary conditions were used as well to save computational time. In all of the models, bottom edge/surface of the plate, raceway, and channel was fixed.

It should be noted that both through-hardened and case-hardened variants of steel were used in both 3D and axisymmetric models. For through-hardened steels, uniform elastic modulus of 200 GPa was used for ball, raceway, plate, and channel material. Since, in open literature, elastic modulus variation data are only available for M50-NiL and P-675 case-hardened steels [2], majority of simulations were performed for these materials. However, hypothetical samples of case-carburized steels were also considered in some of the simulations to generalize the results. These steels can be considered to be processed by pack carburizing, gas carburizing, liquid carburizing, and vacuum/low pressure carburizing methods. These carburization techniques are commonly used by bearing manufacturers [23]. Typical case depths attainable from these manufacturing methods can vary between $50\u2009\mu m\u2009\u2009and\u2009\u20091.5\u2009mm$. The case hardness from these processes is reported in the range of 50–63 Rockwell C hardness (about 20% variation). Corresponding concentrations of carbide volume fraction and elastic moduli variations of the case layer are approximated to follow a similar trend [2]. Generally, for aerospace main shaft bearings, a case depth of 2 mm is common with 10–20% linear decrease in case layer elastic modulus.

*δ*be the case depth and

*d*be the fraction of the elastic modulus variation from $Esurface$ to $Ecore$ over this case depth defined as

Figure 2 depicts elastic modulus variation for M50-NiL material in comparison to through-hardened steel (denoted as TH) as a function of subsurface depth.

While the above gradations were considered only for raceway/plate/channel material, the elastic moduli of the balls were considered to have no gradation, i.e., 320 GPa for silicon nitride balls and 200 GPa for steel balls in all of these models. The spatial variation of elastic modulus in the desired domain of raceway/plate/channel was implemented in ABAQUS via a user-defined material subroutine “user-defined material subroutine,” by assigning numerical values to each material point of the discretization element. In some of the simulations, gradient in elastic modulus was set up using temperature-dependent material properties and gradation in temperature over the case layer of raceway/plate/channel material. Details about the discretization element and mesh techniques used are discussed below.

Computational modeling of nonconformal Hertzian contact between the ball and raceway (or two elastic bodies) is challenging despite significant advances in computational resources. Hence, to obtain a balance between computational cost and accuracy of the FE simulations, mesh discretization scheme, as shown in Fig. 1, was developed for all of the models. For the ball–raceway model, fine mesh was used near the center of elliptical contact area. In this region, smallest dimensions of 3D solid element for raceway geometry were 145.4 *μ*m × 7.2 *μ*m × 3.07 *μ*m (for semimajor axis, semiminor axis, and subsurface depth, respectively). For quarter section of ball geometry, near the center of contact area smallest dimension of 3D solid element were 61.79 *μ*m × 53.73 *μ*m × 24.35 *μ*m. For the raceway, because of steeper gradient, very fine mesh elements were used along semiminor axis *b*. This level of mesh refinement was necessary to determine accurate contact pressure, contact dimensions, and subsurface shear stresses for geometries with uniform elastic modulus of raceway material. Since contact pressure is highly sensitive to the curvature of geometries, quadratic elements were used instead of linear elements [24]. Structured meshing technique along with hexahedral elements (C3D20) was used for the entire FEA of ball–raceway and ball-channel models. With this mesh refinement, approximate time for one 3D ball–raceway contact simulation was 48 h on 32 GB RAM fast processor computer.

The ball-plate model shown in Fig. 1(b) used very refined axisymmetric mesh. At the center of the contact area, smallest dimensions of the two-dimensional element for raceway geometry were 11.11 *μ*m × 6.86 *μ*m and for ball geometry they were 39.90 *μ*m × 34.14 *μ*m. Similar to ball–raceway model, structured meshing technique with quadratic quadrilateral elements (CAX8) were used.

The FEA results were validated against analytical solutions available in Harris [7] and Thomas and Hoersch [17], for the case of through hardened or uniform elastic modulus materials, for test configurations shown in Table 2. The goal of this analysis was to determine the mesh refinement level that would result in acceptable convergence of contact pressure and subsurface stresses for cases where analytical solutions were known. Finite element analysis results shown in Table 2 for maximum contact pressure (*S*) are within 1% of the analytical solution. Table 2 also indicates maximum values of critical subsurface shear stresses and the corresponding depths at which they are observed. For axisymmetric ball-plate model, the errors between analytical and FEA solutions for critical subsurface stresses are less than 2%. The errors in the predicted contact pressure from analytical and computational results were under 3% for all the simulations.

After determining reasonable mesh refinement level, next step is to identify design parameters of interest. As per Hertz theory, these parameters are principal curvatures of two bodies, their material properties, and the normal load acting on the contact. For linearly graded materials such as case-hardened bearing steels, three more design parameters, surface elastic modulus $Esurface$, fractional variation in elastic modulus from surface to core, i.e., *d* (as per Eq. (6)), and case depth *δ*, were defined. Comprehensive FE study was performed by varying each of these design parameters. Table 3 shows details of the 33 different simulated test cases used in this study. Out of these, ten correspond to ball–raceway contact inside 6309 DGBB, 21 correspond to ball-plate model, and two correspond to ball-channel model. Poisson's ratio of 0.3 was kept constant for all the materials. As discussed earlier, most of the test cases were run for a case layer of up to 2 mm depth.

## Results and Discussion

The maximum contact stress (*S*) experienced between the ball and raceway/plate/channel in each of the 33 test configurations are shown in Table 3. For the case-hardened raceway, it is denoted as $SCH$ and that experienced by through-hardened raceway under identical ball material, geometry, and loading conditions is denoted as $STH$. Test cases 1–28 mainly include M50-NiL and P-675 case-carburized bearing steels and test cases 29–33 include different possible variants of steels from carburization process, to generalize the analysis results. From the results presented in Table 3, it can be concluded that peak contact pressure experienced by case-hardened raceway will always be different than that experienced by through-hardened raceway even if other design parameters are kept constant. For example, for test case 2, steel ball rolling over steel raceway, under 1279.1 N load, the peak contact pressure varies by up to 3.61% from STH = 1.663 GPa to SCH = 1.723 GPa. Using Eq. (5), this 3.61% variation in peak pressure will correspond to 37.57% correction in life with *n* = 9. With *n* = 12 or 12.3, the corrections in the life will be even higher, i.e., 53% and 54.64%, respectively. Compared to through-hardened raceway, the peak pressure for carburized raceway is higher in this case due to higher elastic modulus of the material in the near surface region of the raceway. Higher elastic modulus results in smaller deformed area in the contact and hence increased peak contact pressure. Figure 3 shows contact pressure variation along semiminor axis for ball-channel models, i.e., test cases 32 and 33. These curves are obtained by fitting Eq. (1) through contact pressure values at each node along semiminor axis. The coefficients of determination for all the curve fits are 0.99. This indicates that even for case-carburized steel, contact pressure variation along semiminor axis can be approximated using equation similar to traditional Hertz Eq. (1). From Fig. 3, it can also be seen that reduced contact width results in increased peak contact pressure (by about 3%) for case-hardened steel compared to through-hardened steel. With stress-life exponents of 9, 12, and 12.3, this will result in 31.56%, 44.16%, and 45.48% correction in predicted fatigue life, respectively.

Therefore, this underlines the necessity of the fact that gradations in the elastic modulus of case-hardened steel must be considered for fatigue life predictions of bearing raceways. Comparing the contact pressures experienced by case-hardened steel and through-hardened steels, it may appear that through-hardened bearing steels will have higher life than case-hardened bearing steels (as per Eq. (5)). However, in reality, case-carburized bearing steels are known to outperform their through-hardened counterparts under identical operating conditions [22,25]. This improved rolling contact fatigue performance of case-carburized bearing steels is mainly attributed to the presence of residual compressive stresses in circumferential and axial directions of bearing inner rings and fine carbide microstructure in the subsurface region [26]. Analysis on the influence of residual compressive stresses on Hertzian contact stresses in the presence of graded material properties will be reported in a separate article. The main goal of the present analysis is to analyze peak contact pressure for elastically graded bearing raceways and provide simple correction to Hertz equations. For this, nonlinear regression analysis was performed as explained in Statistical Analysis section. For all the test cases presented in Table 3, maximum orthogonal, Tresca, and von Mises stresses were analyzed in the subsurface region of raceway material. Regression estimates for proportionality constant *k* (in Eq. (4)) were 0.247, 0.316, and 0.59, respectively, for ball–raceway contact inside 6309 DGBB. These estimates of *k* for graded material are almost identical to that of uniform elastic modulus/through-hardened material. Therefore, no significant difference was observed in proportionality constant *k* for case-hardened and through-hardened bearing steels. Next mesh convergence study was undertaken to study the influence of mesh discretization schemes on the peak contact pressure values observed in finite element models.

## Mesh Convergence Study

As discussed in the Finite Element Analysis section, reasonable mesh size was determined for three-dimensional and axisymmetric models to minimize the computational cost of finite element simulations. A mesh convergence study was performed to check the influence of mesh size on contact pressure variations for ball–raceway model presented in Table 2. A total of 180,980 second-order 3D solid elements were used, out of which 16,432 were used to discretize the ball geometry and 164,548 were used to discretize the raceway geometry. For convergence study, total mesh size was reduced by an approximate factor of two in each step. Therefore, total numbers of second-order 3D solid elements in four mesh discretization schemes were 180,980, 112,849, 65,144 and 16,365, respectively. Models with through-hardened and M50-NiL steel properties for raceway were simulated using this four mesh discretization schemes. Results obtained for mesh convergence study are shown in Fig. 4, which shows that the variations in peak contact pressure from this four mesh discretization schemes are insignificant. The variation in peak contact pressure for model with 65,144 no. of elements and model with 180,980 no. of elements is less than 0.8% for both the materials. For all the 12 test cases corresponding to the 3D analysis, errors in maximum contact pressures from FE simulations and analytical solutions were less than 3%. Therefore, these confirm that contact pressures have converged satisfactorily. In the present analysis, total of 180,980 no. of elements were used for the models shown in Table 3. Similar analysis was performed for axisymmetric ball-plate model with ball diameter of 12.7 mm. Total number of CAX8 elements selected were 9231, 20,334, 45,400, and 77,475. Each mesh discretization scheme was used to simulate contact of steel ball on steel plate and steel ball on M50-NiL plate. Peak contact pressures observed are plotted in Fig. 5 as function of total number of elements. This plot confirms that peak contact pressures have converged for axisymmetric models.

## Effective Elastic Modulus of Case-Hardened Bearing Steels

For case-carburized raceway, all the terms, except elastic modulus, are constant in Eq. (8).

Equation (9) can be used to determine the effective elastic modulus such that peak Hertz contact pressure for through-hardened raceway material with elastic modulus of $Eeffective$ is $SCH$. It should be noted that Eq. (8) was determined after substituting expressions for semimajor axis and semiminor axis in terms of material properties. Therefore, derived Eq. (9) implicitly accounts for influence of contact dimensions in determining effective elastic modulus of case-hardened steel. Thus, with $Eeffective$, the analytical Hertzian solutions can be used for determining peak contact pressures, subsurface stress fields, and contact dimensions for case-hardened steels such as M50-NiL and P-675. Using Eq. (9), effective elastic modulus of case-carburized steels was determined for each of the 33 test cases. These values are presented in Table 4, which reveal that effective elastic modulus of case-hardened steel is a function of the design parameters of ball–raceway contact area. Especially, geometries with the same material exhibit different effective elastic modulus of case layer under different loads, as can be seen from test cases 2 and 8. For larger dimensions of the ball–raceway contact area, the effective elastic modulus of case-layer is significantly lower than the surface elastic modulus. This is because Hertzian stresses are distributed over larger subsurface volume; hence, effect of gradations in elastic modulus is felt more at the surface. From Table 4, it can also be seen that sharper gradients in elastic modulus, i.e., smaller case depths *δ* or higher gradation *d*, results in significantly lower effective elastic modulus for case-hardened steel than its value at the surface (test cases 29–33). Therefore, these results indicate that effective elastic modulus of case-hardened bearing steel is not only dependent on gradation parameters but is also a function of operating conditions of the bearing raceways. Next step is to determine reasonable surrogate to predict effective elastic modulus of case-hardened steel depending on design parameters of interest, to enable the use of Hertz equations for stress-fatigue life calculations.

## Statistical Analysis

*f*be the function, which represents relationship between these parameters, which can be represented as

*f*in order to establish relationship between design parameters and effective elastic modulus of carburized steel. FE simulation results are numerical approximations with some inherent errors commonly termed as noise. Regression and kriging are two popular choices for surrogate construction, but for noisy data points, regression surrogate is better choice because it can predict approximate global trend in the region that is far away from the current design space. Kriging surrogates are generally preferred for noise free data without any numerical errors. Therefore, regression techniques were used to determine surrogate for

*f*. Certain nondimensional parameters common to Hertzian contacts are used in this analysis [12]. $E*$ is the effective elastic modulus of two bodies in contact defined as

*K*is defined as

Curvature difference of two bodies in contact is nondimensional term; therefore, it is not considered in the definition of *K*. Moreover, its affect was found to be captured by regression coefficients.

*K*were determined for each test case. It was found that linear regression techniques do not provide satisfactory surrogate based on these 33 training data points. However, it was observed that nonlinear regression techniques provide a reasonable surrogate to predict $Eeffective$ as a function of $Esurface$, $Ecore$,

*d,*and

*K*. Using the simplest nonlinear form, this surrogate can be represented as

where $c1,\u2009\u2009c2,\u2009and\u2009c3$ are the constants determined using least-squares fit methods. It should be noted that in Eq. (13), core elastic modulus can be represented as $Ecore$ = $(1\u2212d)Esurface$. The nonlinear regression predicts coefficient of determination, i.e., $R2$ = 0.992 and $Radj2$ = 0.992 for this model. The estimates of coefficients in Eq. (13) are $c1$ = 0.95023, $c2$ = 0.050402, and $c3$ = 0.49188. The P-values for the corresponding estimates are close to 0, indicating that they are statistically significant. The estimated root-mean-square error based on Eq. (13) and the 33 data points presented in Tables 3 and 4 is 1.01 GPa. Due to very high coefficient of determination of close to 1, this seminumerical solution can be considered as correction to Hertz analytical equations for linear elastically graded materials such as case-hardened bearing steels. Next section represents validation of this method to predict peak contact pressure for case-carburized bearing steels.

## Validation Study

To validate this approach of equivalent peak contact pressures, FEA simulations were performed again with effective elastic modulus for raceway/plate/channel material for all the 33 sample test cases presented in Table 3. The results obtained for this cross-validation of peak contact pressures are also presented in Table 4, which shows that this approach works very well for up to 10% and 20% gradations in elastic modulus over 2 mm and 1 mm case depth. For these cases, the errors in actual peak contact pressures experienced by case-hardened steels and the ones obtained using effective elastic modulus of raceway materials are less than 1%. For axisymmetric ball-plate models, all the errors are less than 0.2%, even for up to 9.1% gradation in elastic modulus over 0.5 mm case depth and 16.67% gradation over 1.5 mm case depth (as per test cases # 29–31). Comparing the results for models 2 and 8, we can see that for identical geometry and material properties (i.e., ball–raceway contact with M50-NiL raceway), the effective elastic modulus of case-hardened steel is lower at higher loads. For 1279.1 N load, it is 223.1 GPa, whereas for 2800 N load, it is 221.74 GPa. Therefore, account of gradations in elastic modulus of carburized steels is particularly useful for ball–raceway contact area with larger dimensions and at higher loads. Results for test case #33 indicate that effective elastic modulus is significantly lower than surface elastic modulus for extreme gradations in elastic modulus, i.e., up to 13.04% gradation over 1.5 mm case depth. For ceramic ball-steel raceway contact, the effective elastic modulus is just 215.66 GPa compared to 230–200 GPa variation in elastic modulus over 1.5 mm case layer used in this model. Even for such an extreme design case, validation error is still less than 0.3%. For this test case, contact pressure variation along semiminor axis is also plotted in Fig. 3. From this figure, it can be inferred that contact pressure profile for case-carburized steel with graded material properties is nearly same as that for through-hardened steel with elastic modulus of $Eeffective$. Figure 6 shows plot of 33 data points where peak contact pressures obtained using effective elastic modulus are plotted against peak contact pressures obtained using elastically graded raceway/plate/channel material. As expected, the trend line has an approximate slope of 1 with coefficient of determination: $R2$ = 1, which confirms validity of this method.

## Summary

Many modern bearing steels such as M50-NiL and P-675 are case-hardened, to develop a high fracture toughness core with lower carbon content and hardened surface region with higher carbon content. Due to this gradation in carbide microstructure from surface to core, there exists gradient in material properties such as elastic modulus, hardness, and yield strength over the case layer. Due to high stress-life exponents of 9, 12, or 12.3, bearing fatigue life is very sensitive to variations of peak Hertz pressure experienced by ball–raceway contact. Therefore, determination of accurate Hertz contact pressure experienced by ball–raceway contact is necessary for accurate fatigue life prediction and reliable bearings design. Hertz theory is mainly applicable for homogeneous materials such as through-hardened steels. It cannot be directly used for inhomogeneous materials such as case-carburized bearing steels. In this study, comprehensive finite element analyses were performed to study peak contact pressure for elastically graded materials for different variations of ball–raceway contact design parameters. Because of the inverse elastic gradient, the peak contact pressure experienced by case-hardened steel is different than that experienced by through-hardened steel, under identical geometry and loading. The concept of effective elastic modulus is introduced, which may enable the use of Hertz equations for case-carburized bearing steels. Results from a wide range of simulations are used to determine an accurate regression equation for an effective elastic modulus of case-hardened steel. Using this surrogate, effective elastic modulus can be determined as weighted sum of surface and core elastic moduli of case-hardened bearing steel and geometrical and loading parameters. Even though in present work, gradations in elastic modulus were assumed to be from 230 to 200 GPa or 240 to 200 GPa, this analysis can be adopted for any possible gradations in elastic modulus over different possible case depths. This approach can significantly simplify stress-life analysis for case-carburized bearing steels, as complex, time-consuming, 3D finite element simulations can be avoided.

## Acknowledgment

Nikhil Londhe would like to thank Dr. Anup Pandkar for insightful discussions on this topic.

## Funding Data

The National Science Foundation (NSF) GOALI (Grant No. 1434708).

## Nomenclature

- $a$ =
semimajor axis of elliptical contact area;

*i*and*o*denotes inner and outer raceway contact, respectively - $a*$ =
dimensionless semimajor axis of the contact ellipse

- $b$ =
semiminor axis of elliptical contact area;

*i*and*o*denotes inner and outer raceway contact, respectively - $b*$ =
dimensionless semiminor axis of the contact ellipse

- $cj$ =
nonlinear regression equation coefficients for

*j*ϵ [1, 3] *d*=percent (%) drop in elastic modulus from $Esurface$ to $Ecore$

*D*=ball diameter

- $di$ =
inner raceway diameter

- $dm$ =
bearing pitch diameter

- $do$ =
outer raceway diameter

- $Eball$ =
elastic modulus of ball material

*E*_{core}=elastic modulus of raceway/plate/channel material at core

*E*_{effective}=effective elastic modulus of case hardened bearing steel

*E*_{surface}=surface elastic modulus of raceway/plate/channel material

*E*_{0}=coefficient term in simple power law and exponential variation of elastic modulus

- $E*$ =
effective elastic modulus of two bodies in contact

- $EI\u2009and\u2009EII$ =
elastic modulus of body I and II

*f*=_{i}inner raceway groove radius to ball diameter ratio

*f*=_{o}outer raceway groove radius to ball diameter ratio

- $Fe$ =
equivalent dynamic load

- $Fr$ =
applied radial load

- $F(\rho )$ =
curvature difference of two bodies in contact

- $F(\rho )i$ =
curvature difference at inner raceway contact

- $F(\rho )o$ =
curvature difference at outer raceway contact

*k*=proportionality constant for maximum shear stress and peak contact pressure relation

*K*=nondimensional model design parameter

*k*′ =depth exponent in power-law variation in elastic modulus

*L*=bearing fatigue life, millions of revolutions or millions of stress cycles

- $Li$ =
fatigue life of inner raceway, millions of revolutions or millions of stress cycles

- $Lo$ =
fatigue life of outer raceway, millions of revolutions or millions of stress cycles

- $L10$ =
bearing fatigue life in millions of revolutions corresponding to 90% survival probability

- LF =
elastic modulus gradation life factor

*n*=stress-life exponent

*N*=rpm

*p*=load-life exponent

*Q*=normal load experienced by the contact

- $Qc$ =
basic dynamic capacity of raceways;

*i*and*o*denote inner and outer raceway contacts, respectively - $Qe$ =
cubic mean equivalent radial load experienced by raceway;

*i*and*o*denote inner and outer raceway contacts, respectively - $Qmax$ =
maximum load experienced by ball-raceway contact

- $ri$ =
inner raceway groove radius

- $ro$ =
outer raceway groove radius

- $R2$ =
coefficient of determination

- $Radj2$ =
adjusted coefficient of determination

*s*=probability of survival

*S*=maximum compressive stress/peak contact stress

- $SCH$ =
maximum contact pressure for graded/ case-hardened raceway

- $Smaxi$ =
maximum compressive hertz stress experienced by inner raceway

- $Smaxo$ =
maximum compressive hertz stress experienced by outer raceway

- $STH$ =
maximum contact pressure with uniform elastic modulus/through-hardened raceway material

*x*=*x*-coordinate*y*=*y*-coordinate*z*=subsurface depth below the center of contact area

*Z*=number of rolling elements/balls in ball bearings

*α*=bearing contact angle

*δ*=case depth

- $\xi I\u2002and\u2002\xi II$ =
Poisson's ratio of body I and II

- $\rho 1I\u2002and\u2002\rho 2I$ =
principal curvatures of body I

- $\rho 1II\u2002and\u2002\rho 2II$ =
principal curvatures of body II

*σ*=Hertz pressure in the contact region

- $\u2211\rho $ =
curvature sum of two bodies in contact

- $\u2211\rho i$ =
curvature sum at inner raceway contact

- $\u2211\rho o$ =
curvature sum at outer raceway contact

- $\tau max$ =
maximum shear/Tresca stress

- $\tau oct$ =
maximum octahedral shear stress

- $\tau 0$ =
orthogonal shear stress

- $\varphi i$ =
inner contact osculation

- $\varphi o$ =
outer contact osculation

### Appendix

To demonstrate the influence of elastic modulus variations on the bearing fatigue lives, a 209 single row deep groove ball bearing example from Londhe [21] will be used. Detailed information about its component geometries and material properties is provided in Table 5.

After determining the geometrical properties of the bearing, contact stresses and deformations can be determined.

##### Case A: Elastic Modulus of Raceway Material = 200 GPa.

##### Case B: Elastic Modulus of Raceway Material = 180 GPa.

Solving Eq. (A5), we get fatigue life of the inner and outer raceways of the bearing in case B as, $Li$ = 31.55 million cycles and $Lo$ = 213.92 million cycles, respectively. From Eq. (A4), combined fatigue life of the bearing in case B is 28.404 million cycles. Therefore, if we decrease elastic modulus of the raceway material by 10% then, the combined bearing fatigue life is corrected by 38.66% from case A to case B. Hence, this analysis confirms that bearing fatigue lives are highly sensitive to the variations in elastic modulus of raceways and there is need to account for gradient in elastic modulus for accurate prediction of fatigue lives for case-hardened bearing steels. It should be noted that in this analysis, elastic modulus of 180 GPa is used just as a representative case. The actual variations in elastic modulus for case-carburized steels are from 230 to 200 GPa over different possible case depths as considered in FEA.