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

Traditional Hertzian contact solution for determining subsurface elastic stress fields for two bodies in contact is applicable only for homogeneous materials [7]. Bearing L10 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

Stresses developed at the contact of two elastic solids were first analyzed by Hertz in early 1880s. He proposed that contact area is, in general, elliptical and for the purpose of calculating local deformations, both bodies can be approximated as elastic half spaces loaded over a small elliptical region of their surfaces [12]. Based on Hertz analysis, ellipsoidal compressive stress distribution in the contact area can be written as
$σ=3Q2πab[1−(xa)2−(yb)2]12$
(1)
where 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]
$a=a*[3Q2∑ρ((1−ξI2)EI+(1−ξII2)EII)]13$
(2)

$b=b*[3Q2∑ρ((1−ξI2)EI+(1−ξII2)EII)]13$
(3)

where $ξ$ and E are Poisson's ratio and elastic modulus of bodies I and II, respectively; $∑ρ$ is the summation of principal curvatures of two bodies in contact; $a* and b*$ are nondimensional parameters defined by curvature difference of two bodies [7].

It is well known that bearings generally fail due to rolling contact fatigue arising from continuously varying subsurface shear stresses [7]. Three principal shearing stresses are generally used for bearing life analysis: the orthogonal shearing stress, $τ0$ (based on Lundberg–Palmgren model [13]); the von Mises equivalent octahedral shearing stress, $τoct$ (based on Ioannides and Harris model [14]); and the maximum shearing stress, $τmax$ (based on Zaretsky's model [15]). All these shearing stresses vary in proportion with the maximum Hertz pressure $S$ experienced by two bodies as
$τ=k×S$
(4)

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)≈$ 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.5b, 0.72b, 0.76b, respectively. Corresponding depths for roller bearings are 0.5b, 0.79b, 0.79b, 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

It is well known that the bearing fatigue life, L, is inversely proportional to equivalent dynamic load $Fe$ and also maximum Hertzian contact pressure S as [11]
$L∝1[Fe]p∝1Sn$
(5)

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 = 3p) [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 = 3p = 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 μm and 1.5 mm$. 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.

As discussed earlier, to study the influence of these graded properties of raceway material on the peak contact pressure under elliptical contact loading conditions, a wide range of operating conditions, as detailed in Table 3, were considered in the finite element analysis (FEA). The 3D ball–raceway and ball–channel models and axisymmetric ball-plate model (as shown in Fig. 1) were simulated for nominal loads ranging from 100 N to 13,500 N. At each load, simulations were performed for silicon nitride and steel balls in contact with M50-NiL and P-675 steel raceway/plate material. The resulting peak contact pressures were in the range of 1.5–3.1 GPa, which is the typical operating range for majority of industrial ball-bearing applications. Additionally, case depths of 0.5 mm and 1.5 mm were used in some of the simulations (model nos. 29–33 in Table 3) with different possible gradations in elastic modulus from surface to core region. Let $Esurface$ and $Ecore$ denote surface and core elastic moduli, respectively, of the raceway/plate/channel materials. Let δ be the case depth and d be the fraction of the elastic modulus variation from $Esurface$ to $Ecore$ over this case depth defined as
$d=(Esurface−EcoreEsurface)$
(6)

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

The composite elastic modulus of case layer decreases linearly from surface to core material. In order to use Hertz equations for such graded material, effective elastic modulus of the case layer must be determined. This effective elastic modulus can also be used in bearing fatigue life prediction equations, which are mainly applicable for through-hardened steels/homogeneous materials. To determine this modulus, an equivalent contact pressure approach will be used with traditional Hertz equations. Using Eq. (1), peak contact pressure can be expressed as
$S=3Q2πab$
(7)
Substituting expressions for a and b (i.e., Eqs. (2) and (3)) in Eq. (7), we get
$S=(3Q)132πa*b*[12∑ρ((1−ξI2)EI+(1−ξII2)EII)]−23$
(8)

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

Let $Eball=EI$ be the elastic modulus of ball and $ETH=EII$ be the elastic modulus of through-hardened steel raceway. $STH$ represents the peak contact pressure for the raceway with uniform elastic modulus of $ETH$. Let $SCH$ represent peak contact pressure for graded case-hardened raceway material. In each of the 33 test cases, contact pressures for case-hardened and through-hardened raceway under similar operating conditions were determined and are included in Table 3. Let $Eeffective$ represent the effective elastic modulus of the case layer, which predicts peak contact pressure of $SCH$ from traditional Hertz Eq. (8). Therefore, using Eq. (8), for identical operating conditions, these terms can be related as
$STHSCH=[(1Eball+1ETH)(1Eball+1Eeffective)]−23$
(9)

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

In this section, the aim is to determine mathematical relationship between effective elastic modulus of case layer, and geometry, material, elastic modulus gradation and load parameters of ball–raceway contact. Let f be the function, which represents relationship between these parameters, which can be represented as
$Eeffective=f(Eball,Esurface,d,δ,Q,ξI,ξII,∑ρ,F(ρ))$
(10)
Surrogate/meta models were constructed to determine an approximate functional form for 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
$1E*=(1−ξI2)Eball+(1−ξII2)Esurface$
(11)
The nondimensional model design parameter K is defined as
$K=1δ(2QE*∑ρ)13$
(12)

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.

matlab software program was used for the regression analysis to determine curve fit equation based on the 33 training data points presented in Tables 3 and 4. Values of the nondimensional model design parameter 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
$Eeffective=c1Esurface+c2Ecore1+c3dK$
(13)

where $c1, c2, and c3$ 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−d)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.

## 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

•
• Ecore =

elastic modulus of raceway/plate/channel material at core

•
• Eeffective =

effective elastic modulus of case hardened bearing steel

•
• Esurface =

surface elastic modulus of raceway/plate/channel material

•
• E0 =

coefficient term in simple power law and exponential variation of elastic modulus

•
• $E*$ =

effective elastic modulus of two bodies in contact

•
• $EI and EII$ =

elastic modulus of body I and II

•
• fi =

inner raceway groove radius to ball diameter ratio

•
• fo =

outer raceway groove radius to ball diameter ratio

•
• $Fe$ =

•
• $Fr$ =

•
• $F(ρ)$ =

curvature difference of two bodies in contact

•
• $F(ρ)i$ =

curvature difference at inner raceway contact

•
• $F(ρ)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 =

•
• n =

stress-life exponent

•
• N =

rpm

•
• p =

•
• 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$ =

•
• $ro$ =

•
• $R2$ =

coefficient of determination

•
• $Radj2$ =

•
• 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

•
• $ξI and ξII$ =

Poisson's ratio of body I and II

•
• $ρ1I and ρ2I$ =

principal curvatures of body I

•
• $ρ1II and ρ2II$ =

principal curvatures of body II

•
• σ =

Hertz pressure in the contact region

•
• $∑ρ$ =

curvature sum of two bodies in contact

•
• $∑ρi$ =

curvature sum at inner raceway contact

•
• $∑ρo$ =

curvature sum at outer raceway contact

•
• $τmax$ =

maximum shear/Tresca stress

•
• $τoct$ =

maximum octahedral shear stress

•
• $τ0$ =

orthogonal shear stress

•
• $ϕi$ =

inner contact osculation

•
• $ϕ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.

Using standard bearing macrogeometry relations, we can determine the osculation values for the inner and outer raceways as
$fi=riD=0.52; φi=12fi=0.962$

$fo=roD=0.52; φo=12fo=0.962$
For inner raceway-ball contact, curvature sum, and curvature difference were calculated using principal radii of curvatures for two bodies $ρ1I=ρ2I$  =  $(2/D)$, $ρ1II$  =  $(2/di)$ and $ρ2II$  =  $(1/fiD)$, respectively,
$∑ρi=ρ1I+ρ2I+ρ1II+ρ2II=0.2018 mm−1F(ρ)i=(ρ1I−ρ2I)+(ρ1II−ρ2II)∑ρ=0.93997$
(A1)
Similarly for outer raceway-ball contact, curvature sum and curvature difference can be calculated as
$∑ρo=0.1378 mm−1$

$F(ρ)o=0.91209$

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

##### Case A: Elastic Modulus of Raceway Material = 200 GPa.
To determine ball-raceway contact dimensions, matlab subroutines were developed to solve for complete elliptical integrals of first and second kind, and load distribution integrals [7]. Maximum load experienced by ball-raceway contact was found to be $Qmax$  = 4527.88 N. At this load, inner and outer raceway contact dimensions were determined to be
$ai=2.591 mm, bi=0.2772 mm, ao=2.5054 mm, bo=0.3424 mm$
Using traditional Hertz solutions, peak compressive stress experienced by inner and outer raceways were determined as $Smaxi$ = 3011.59 MPa and $Smaxo$ = 2521.42 MPa. Basic dynamic capacities for inner and outer raceways were determined using following equations [7]:
$C=93.2[2f2f−1]0.41(1∓γ)1.39(1±γ)13[γ cos α]0.3D1.8Z−13$
(A2)
where $γ=(D cos α)/dm$. Using Eq. (A2), basic dynamic capacity for inner and outer raceway contact was found to be $Qci$ = 7054.95 N and $Qco$ = 13,958.14 N, respectively. Cubic mean equivalent radial load [7] experienced by rotating inner and outer raceways is $Qei$ = 2489.76 N and $Qeo$ = 2604.82 N, respectively. In this analysis, Lundberg–Palmgren model is used to compare fatigue lives obtained from elastic modulus variations
$L=[QcQe]3$
(A3)
Therefore, based on dynamic capacity and equivalent radial load, we can calculate fatigue lives for both inner and outer raceways using Eq. (A3) as
$Li=[QciQei]3=22.75 million cycles$

$Lo=[QcoQeo]3=153.87 million cycles$
Combined life of bearing can be obtained from following equation [7]:
$L10=(Li−1.11+Lo−1.11)−0.9=20.485 million cycles=189.672 h$
(A4)
##### Case B: Elastic Modulus of Raceway Material = 180 GPa.
Now, let us assume that elastic modulus of the bearing raceway material is decreased by 10% from 200 GPa to 180 GPa. Also, it will be assumed that Poisson's ratio of the material and elastic modulus of the ball remains constant. Using matlab subroutine, maximum load ($Qmax$) experienced by ball-raceway contact was determined to be 4521.44 N. At this load, inner and outer contact dimensions were found to be
$ai=2.6364 mm, bi=0.2821 mm, ao=2.5497 mm, bo=0.3485 mm$
Peak compressive Hertz stress experienced by inner and outer raceways were $Smaxi$ = 2904.18 MPa and $Smaxo$ = 2430.79 MPa, respectively. Using Eq. (5) and stress life exponent of 9, we can relate fatigue lives with peak hertz stresses experienced by the raceways in cases A and B as
$LBLA=[SmaxASmaxB]9$
(A5)

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.

## References

References
1.
Giannakopoulos
,
A. E.
, and
Suresh
,
S.
,
1997
, “
Indentation of Solids With Gradients in Elastic Properties—Part I: Point Force
,”
Int. J. Solids Struct.
,
34
(
19
), pp.
2357
2392
.
2.
Klecka
,
M. A.
,
Subhash
,
G.
, and
Arakere
,
N. K.
,
2013
, “
Microstructure-Property Relationships in M50-NiL and P675 Case-Hardened Bearing Steels
,”
Tribol. Trans.
,
56
(
6
), pp.
1046
1059
.
3.
Suresh
,
S.
,
2001
, “
Graded Materials for Resistance to Contact Deformation and Damage
,”
Science
,
292
(
5526
), pp.
2447
2451
.
4.
Giannakopoulos
,
A. E.
, and
Suresh
,
S.
,
1997
, “
Indentation of Solids With Gradients in Elastic Properties—Part II: Axisymmetric Indentors
,”
Int. J. Solids Struct.
,
34
(
19
), pp.
2393
2428
.
5.
Liu
,
T.-J.
,
Wang
,
Y.-S.
, and
Zhang
,
C.
,
2008
, “
Axisymmetric Frictionless Contact of Functionally Graded Materials
,”
Arch. Appl. Mech.
,
78
(
4
), pp.
267
282
.
6.
Guler
,
M. A.
, and
Erdogan
,
F.
,
2004
, “
,”
Int. J. Solids Struct.
,
41
(
14
), pp.
3865
3889
.
7.
Harris
,
T. A.
,
1991
,
Rolling Bearing Analysis
, 3rd ed.,
Wiley
,
New York
, Chaps. 2, 5, and 6.
8.
Zaretsky
,
E. V.
,
1992
,
STLE Life Factors for Rolling Bearings
,
Society of Tribologists and Lubrication Engineers
,
Park Ridge, IL
.
9.
Bhattacharyya
,
A.
,
Londhe
,
N. D.
,
Arakere
,
N. K.
, and
Subhash
,
G.
,
2017
, “
A New Approach Towards Life Prediction of Case Hardened Bearing Steels Subjected to Rolling Contact Fatigue
,”
ASTM Int. J. Mater. Perform. Charact.
,
6
(
4
), pp.
1
22
.
10.
Weinzapfel
,
N.
,
,
F.
, and
Bakolas
,
V.
,
2011
, “
A 3D Finite Element Model for Investigating Effects of Material Microstructure on Rolling Contact Fatigue
,”
Tribol. Lubr. Technol.
,
67
(
1
), pp.
17
19
.
11.
Zaretsky
,
E. V.
,
Vlcek
,
B. L.
, and
Hendricks
,
R. C.
,
2005
, “
Effect of Silicon Nitride Balls and Rollers on Rolling Bearing Life
,” NASA Glenn Research Center, Cleveland, OH, Report No.
NASA/TM-2005-213061
.
12.
Johnson
,
K. L.
,
2001
,
Contact Mechanics
,
Cambridge University Press
,
Cambridge, UK
, Chap. 4.
13.
Lundberg
,
G.
, and
Palmgren
,
A.
,
1947
, Dynamic Capacity of Rolling Bearings (Acta Polytechnica. Mechanical Engineering Series), Vol. 1, Generalstabens Litografiska Anstalts Förl., Stockholm, Sweden.
14.
Ioannides
,
E.
, and
Harris
,
T. A.
,
1985
, “
A New Fatigue Life Model for Rolling Bearings
,”
ASME J. Tribol.
,
107
(
3
), pp.
367
378
.
15.
Zaretsky
,
E. V.
,
Poplawski
,
J. V.
, and
Peters
,
S. M.
,
1996
, “
Comparison of Life Theories for Rolling Element Bearings
,”
Tribol. Trans
.,
39
(
2
), pp.
237
248
.
16.
Zaretsky
,
E. V.
,
2012
, “
Rolling Bearing Steels—A Technical and Historical Perspective
,”
J. Mater. Sci. Technol.
,
28
(
1
), pp.
58
69
.
17.
Thomas
,
H. R.
, and
Hoersch
,
V. A.
,
1930
, “
Stresses Due to the Pressure of One Elastic Solid Upon Another
,” Engineering Experiment Station, University of Illinois at Urbana Champaign, Urbana, IL, Bulletin No.
212
.
18.
Parker
,
R. J.
, and
Zaretsky
,
E. V.
,
1972
, “
Reevaluation of the Stress-Life Relation in Rolling-Element Bearings
,” NASA Lewis Research Center, Cleveland, OH, Report No.
NASA-TN-D-6745
.
19.
Londhe
,
N. D.
,
Arakere
,
N. K.
, and
Haftka
,
R. T.
,
2015
, “
Reevaluation of Rolling Element Bearing Load-Life Equation Based on Fatigue Endurance Data
,”
Tribol. Trans.
,
58
(
5
), pp.
815
828
.
20.
Harris
,
T. A.
, and
McCool
,
J. I.
,
1996
, “
On the Accuracy of Rolling Bearing Fatigue Life Prediction
,”
ASME J. Tribol.
,
118
(
2
), pp.
297
310
.
21.
Londhe
,
N. D.
,
2014
, “
Modification of Standard Load Life Equation Used for Designing Rolling Element Bearings
,”
Master's thesis
, University of Florida, Gainesville, FL.
22.
Harris
,
T. A.
,
Skiller
,
J.
, and
Spitzer
,
R. F.
,
1992
, “
On the Fatigue Life of M50 NiL Rolling Bearings
,”
Tribol. Trans.
,
35
(
4
), pp.
731
737
.
23.
Schneider
,
M. J.
, and
Chatterjee
,
M. S.
,
2013
, “
Introduction to Surface Hardening of Steels
,” Steel Heat Treating Fundamentals and Processes, ASM Handbook, Vol. 4A, ASM International, Geauga County, OH, pp.
259
267
.
24.
Kim
,
N. H.
,
2015
,
Introduction to Nonlinear Finite Element Analysis
,
Springer
, New York, Chap. 5.
25.
Errichello
,
R.
,
Bundy
,
R.
, and
Eckert
,
R.
,
2013
, “
Investigations of Bearing Failures Associated With White Etching Areas (WEAs) in Wind Turbine Gearboxes
,”
Tribol. Trans.
,
56
(
6
), pp.
1069
1076
.
26.
Philip
,
T. V.
,
1986
, “
New Bearing Steels Beats Speed and Heat
,”
Power Transm. Des.
,
28
(
6
), pp.
43
46
.