The relentless pursuit of higher power density and extended service life in modern mechanical transmission systems places immense demands on their core component: the gear. Among the various failure modes threatening gear reliability, contact fatigue, manifesting as micropitting, pitting, and spalling, stands as a critical bottleneck. The conventional design approach, codified in standards like ISO 6336 and GB/T 3480, relies on calculating nominal contact stress and comparing it to a homogeneous material strength. While effective for initial sizing, this methodology often falls short in accurately predicting the complex fatigue behavior of modern case-hardened gears, where a gradient of mechanical properties—notably hardness and residual stress—profoundly influences the failure initiation site and life.
This inherent limitation underscores the industrial and academic imperative to develop more sophisticated predictive models. The key lies in deeply characterizing the correlation mechanism between the gradient load-bearing capacity of the carburized and quenched modified layer and the multiaxial stress-strain response induced by meshing contact. Such a model must transcend homogeneous material assumptions and explicitly account for the spatially varying material strength resulting from thermochemical treatments. This article, therefore, presents a comprehensive study on a coupled predictive model for contact fatigue risk in carburized cylindrical gear, integrating gear geometry, multiaxial contact mechanics, and material gradient effects into a unified framework.
Fundamental Contact Mechanics for Spur Gear Meshing
The analysis begins by establishing the mechanical framework for contact between two gear teeth. For a pair of spur gears with a contact ratio $\epsilon$ between 1 and 2, the load alternates between single and double tooth pairs along the line of action. At any instantaneous contact point, the geometry can be equivalently represented by two contacting cylinders with radii equal to the local curvatures of the gear tooth profiles.
Equivalent Cylinder Model and Hertzian Parameters
Based on the involute geometry and Hertzian elastic contact theory, the principal parameters for a given meshing point, characterized by its pressure angle $\alpha_{pm}$ or $\alpha_{qm}$, are derived. The governing equations for the equivalent contact are:
$$
F = \frac{9549 \cdot P_0}{n \cdot R_{pm} \cos(\alpha_{pm})}
$$
$$
b = \sqrt{ \frac{4F}{\pi B} \cdot \frac{\rho_{pm} \rho_{qm}}{\rho_{pm} + \rho_{qm}} \cdot \left( \frac{1-\nu_p^2}{E_p} + \frac{1-\nu_q^2}{E_q} \right) }
$$
$$
\rho_{pm} = \frac{m z_p \cos(\alpha) \tan(\alpha_{pm})}{2}; \quad \rho_{qm} = \frac{m z_q \cos(\alpha) \tan(\alpha_{qm})}{2}
$$
Where $F$ is the normal load (N), $P_0$ is transmitted power (kW), $n$ is rotational speed (rpm), $R_{pm}$ is the radius at the contact point on the pinion (m), $b$ is the semi-contact width (m), $B$ is the face width (m), $\rho_{pm}, \rho_{qm}$ are the radii of curvature at the contact point for pinion and gear respectively (m), $\nu$ and $E$ are Poisson’s ratio and Young’s modulus, $m$ is the module (m), $\alpha$ is the standard pressure angle, and $z_p, z_q$ are the numbers of teeth.
Stress Field in an Elastic Half-Plane
The contact problem is treated as a semi-infinite elastic body subjected to combined normal and tangential tractions over a finite surface region $[-a, a]$. For an arbitrary distribution of normal pressure $p(t)$ and tangential traction $q(t)$, the stress components at any subsurface point $(x, z)$ are given by the following singular integrals:
$$
\begin{aligned}
\sigma_x(x,z) &= -\frac{2z}{\pi} \int_{x_1}^{x_2} p(t) \frac{(x-t)^2}{[(x-t)^2+z^2]^2} dt – \frac{2}{\pi} \int_{x_1}^{x_2} q(t) \frac{(x-t)^3}{[(x-t)^2+z^2]^2} dt \\
\sigma_z(x,z) &= -\frac{2z^3}{\pi} \int_{x_1}^{x_2} \frac{p(t)}{[(x-t)^2+z^2]^2} dt – \frac{2z^2}{\pi} \int_{x_1}^{x_2} q(t) \frac{(x-t)}{[(x-t)^2+z^2]^2} dt \\
\tau_{xz}(x,z) &= -\frac{2z^2}{\pi} \int_{x_1}^{x_2} p(t) \frac{(x-t)}{[(x-t)^2+z^2]^2} dt – \frac{2z}{\pi} \int_{x_1}^{x_2} q(t) \frac{(x-t)^2}{[(x-t)^2+z^2]^2} dt
\end{aligned}
$$
Discrete Numerical Solution via Superposition
Direct analytical integration of the above equations for complex traction distributions is challenging. A robust numerical approach is adopted, discretizing the contact region into $n$ small rectangular elements, each carrying a uniform normal load $p_j$ and tangential load $q_j$. The stress field due to a single uniform rectangular load is known in closed form. For a normal load $p_j$ acting on an element centered at the origin with half-width $s$, the stresses at point $(x_{ij}, z_j)$ in a local coordinate system are:
$$
\begin{aligned}
\sigma_x^{p} &= -\frac{p_j}{\pi} \left[ \arctan\left(\frac{s-x_{ij}}{z_j}\right) + \arctan\left(\frac{s+x_{ij}}{z_j}\right) – \frac{z_j(s+x_{ij})}{r_1^2} + \frac{z_j(x_{ij}-s)}{r_2^2} \right] \\
\sigma_z^{p} &= -\frac{p_j}{\pi} \left[ \arctan\left(\frac{s-x_{ij}}{z_j}\right) + \arctan\left(\frac{s+x_{ij}}{z_j}\right) + \frac{z_j(s+x_{ij})}{r_1^2} – \frac{z_j(x_{ij}-s)}{r_2^2} \right] \\
\tau_{xz}^{p} &= -\frac{z_j^2 p_j}{\pi} \left( \frac{1}{r_1^2} – \frac{1}{r_2^2} \right)
\end{aligned}
$$
with $r_1^2 = (s-x_{ij})^2+z_j^2$ and $r_2^2 = (s+x_{ij})^2+z_j^2$. Similar closed-form expressions exist for a uniform tangential load $q_j$. By defining influence coefficient functions $T_{p,j-i}$ and $T_{q,j-i}$ that are independent of the load magnitude, the total stress at any point $i$ due to all $n$ loaded elements is obtained through linear superposition:
$$
\sigma_{x}^i = \sum_{j=1}^{n} \left( T_{p,j-i} \cdot p_j + T_{q,j-i} \cdot q_j \right)
$$
The same principle applies to $\sigma_z$ and $\tau_{xz}$. This method provides an efficient and accurate way to compute the complete subsurface multiaxial stress field for any given surface traction distribution on the cylindrical gear contact.
Characterization of Performance Gradients in Carburized Gears
The defining feature of a carburized and quenched cylindrical gear is the gradient in mechanical properties from the case to the core. Accurate mathematical representation of these gradients is essential for the fatigue model.
Hardness Gradient Model
The carbon concentration gradient after carburizing results in a corresponding hardness gradient. The measured hardness profile of a case-hardened cylindrical gear can be effectively modeled using a piecewise quadratic function anchored to key parameters: surface hardness ($HV_{surface}$), core hardness ($HV_{core}$), and the effective case hardening depth ($D_{CHD}$, typically measured at 550 HV).
$$
HV(z) =
\begin{cases}
a_a z^2 + b_a z + c_a, & 0 \leq z < D_{CHD} \\
a_b z^2 + b_b z + c_b, & D_{CHD} \leq z < z_{core} \\
HV_{core}, & z_{core} \leq z
\end{cases}
$$
Where $z$ is the depth from the surface, and coefficients $a, b, c$ are determined by fitting the model to measured data or standard profile shapes. This representation captures the steep gradient in the case and the transition to the core. The tensile strength $\sigma_{UTS}(z)$ is directly correlated to hardness via an approximately linear relationship established from material standards, allowing the strength gradient to be derived from the hardness gradient.
Residual Stress Gradient Model
The quenching process introduces significant residual stresses, typically compressive at the surface and tensile in the core. A statistically-derived predictive model for the residual stress profile $\sigma_{RS}(z)$ is employed, taking the form of a logistic function:
$$
\sigma_{RS}(z) = \sigma_D + \frac{\sigma_z – \sigma_D}{1 + e^{-k(z + \delta)}}
$$
Here, $\sigma_D$ is the maximum compressive residual stress (a negative value), $\sigma_z$ is the maximum tensile residual stress in the core, and $k$ and $\delta$ are fitting parameters controlling the slope and inflection point of the gradient curve. This model reliably predicts the characteristic profile, which is crucial as residual stresses directly superimpose on the applied contact stresses, altering the mean stress state experienced by the material.
| Parameter | Symbol | Value |
|---|---|---|
| Module | $m$ | 6.5 mm |
| Number of Teeth (Pinion/Gear) | $z_p / z_q$ | 24 / 25 |
| Pressure Angle | $\alpha$ | 20° |
| Face Width | $B$ | 30 mm |
| Surface Hardness | $HV_{surface}$ | ~700 HV |
| Core Hardness | $HV_{core}$ | ~450 HV |
| Effective Case Depth | $D_{CHD}$ | Variable |
| Young’s Modulus | $E$ | 210 GPa |
| Poisson’s Ratio | $\nu$ | 0.3 |
Multiaxial Fatigue Risk Assessment Using the Dang Van Criterion
To predict fatigue crack initiation under the complex, non-proportional multiaxial stress state in rolling-sliding contact, a critical plane approach is necessary. The Dang Van criterion is particularly suited for high-cycle rolling contact fatigue as it accounts for the beneficial effect of compressive hydrostatic stress.
Criterion Formulation
The Dang Van criterion postulates that fatigue crack initiation occurs on microscopic slip systems when a linear combination of the alternating shear stress $\Delta \tau_{\max}(\alpha, t)$ on a material plane and the contemporaneous hydrostatic stress $\sigma_H(t)$ reaches a critical value. The fatigue risk parameter, or the material’s resistance utilization factor $R_{FP}$, is defined for a plane oriented at angle $\alpha$ at time $t$ as:
$$
R_{FP}(\alpha, t) = \frac{\Delta \tau_{\max}(\alpha, t) + \kappa_D \cdot \sigma_H(t)}{\lambda}
$$
Where $\lambda$ is the shear fatigue limit of the material (related to its tensile fatigue limit $\sigma_{-1}$), and $\kappa_D$ is a material coefficient linking the sensitivity to hydrostatic stress, determined from fully reversed bending and torsion limits ($\sigma_{-1}$ and $\tau_{-1}$). For a given stress history at a material point, all possible material planes are examined. The plane experiencing the maximum value of $R_{FP}$ over a load cycle is identified as the critical plane, and its maximum $R_{FP}^{\max}$ value indicates the fatigue risk. Failure is predicted when $R_{FP}^{\max} \geq 1$.
Integration of Material Gradients
The model’s power for a cylindrical gear analysis lies in incorporating the spatial gradients. Both the material strength $\lambda(z)$ and the residual stress $\sigma_{RS}(z)$ are functions of depth $z$. The fatigue risk parameter thus becomes:
$$
R’_{FP}(\alpha, z, t) = \frac{\Delta \tau_{\max}(\alpha, z, t) + \kappa_D(z) \cdot [\sigma_H(z, t) + \sigma_{RS}(z)]}{\lambda(z)}
$$
This equation represents the fully coupled predictive model. For each subsurface point $(x, z)$, the discrete numerical method calculates the time-varying stress tensor $\boldsymbol{\sigma}(x,z,t)$ throughout a meshing cycle. This stress history is then processed through the gradient-aware Dang Van criterion to compute the local $R’_{FP}^{\max}(x,z)$. Mapping this value across the $(x, z)$ domain reveals the spatial distribution of fatigue risk within the cylindrical gear tooth volume.
Results, Analysis, and Parametric Studies
Applying the coupled model to the cylindrical gear pair described in Table 1, with a nominal contact stress of 1500 MPa at the pitch point and a coefficient of friction $\mu=0.05$ to simulate lubricated rolling-sliding, yields fundamental insights.
Subsurface Stress State
The computed stress fields confirm classical contact mechanics. Under pure normal loading, the principal stresses $\sigma_1, \sigma_2, \sigma_3$ are compressive near the surface, decaying with depth. The maximum shear stress $\tau_{\max}$, however, peaks not at the surface but in the subsurface region (at approximately $z \approx 0.3b$). The addition of tangential traction from friction slightly alters the magnitude and location of these extremes. The hydrostatic pressure $\sigma_H$ is highly compressive at the surface, while the alternating shear stress amplitude $\Delta\tau_{\max}$ exhibits its own distinct gradient.
Fatigue Risk Mapping Without Residual Stress
Initial calculation of $R_{FP}$ using only the hardness gradient (no residual stress) shows that the zone of maximum risk does not coincide with the point of maximum shear stress. Due to the rolling-sliding kinematics and the decaying material strength $\lambda(z)$, the highest $R_{FP}^{\max}$ value is located at a point offset from the contact center and slightly deeper than the $\tau_{\max}$ peak (e.g., around $x = -0.1$ mm, $z = 0.4$ mm with a value of 0.83). This demonstrates that the risk is governed by the interplay between the applied multiaxial stress and the local material strength, not by a single stress component.
The Critical Role of Residual Stresses
Incorporating the characteristic residual compressive stress profile fundamentally changes the risk landscape. The compressive $\sigma_{RS}(z)$ adds to the applied compressive $\sigma_H$, significantly increasing the denominator in the Dang Van equation. Consequently, the maximum $R_{FP}^{\max}$ value at the previously identified critical location drops substantially, for instance, from 0.83 to 0.71. This quantifies the potent beneficial effect of compressive residual stresses in retarding the initiation of subsurface-origin fatigue cracks (pitting) in a cylindrical gear. The risk contour plot shows a general depression of $R_{FP}$ values across the near-surface region.
Influence of Key Characteristic Parameters
The model enables systematic study of how various factors influence the contact fatigue risk of a cylindrical gear.
Friction Coefficient (Lubrication Condition): Increasing $\mu$ from 0 (pure rolling) to 0.2 (high sliding friction) simulates deteriorating lubrication. The risk map shifts dramatically: the location of maximum $R_{FP}^{\max}$ migrates both towards the surface and towards the trailing edge of the contact. Furthermore, its magnitude increases significantly.
Surface Hardness (Case Strength): Increasing $HV_{surface}$ from 600 HV to 750 HV directly elevates the near-surface material strength $\lambda(z)$. This results in a uniform reduction of $R_{FP}^{\max}$ values across the case depth without altering the spatial pattern of the risk zone. It confirms that increasing surface hardness is an effective way to boost the contact fatigue resistance of a cylindrical gear.
Contact Load: Increasing the nominal contact pressure from 1500 MPa to 1900 MPa expands the contact semi-width and intensifies the stress field. The model predicts the emergence of a region where $R_{FP}^{\max} > 1$, indicating a high probability of fatigue crack initiation. This aligns perfectly with the well-known exponential relationship between applied load and fatigue life.
| Parameter Variation | Effect on Max $R_{FP}^{\max}$ Location | Effect on Max $R_{FP}^{\max}$ Magnitude | Implied Failure Mode Shift |
|---|---|---|---|
| Friction $\mu$ Increase | Migrates towards surface & trailing edge | Significant Increase | Subsurface pitting risk increases; surface-initiated micropitting becomes more likely. |
| Surface Hardness Increase | Minimal change | Decrease | Overall fatigue resistance improves. |
| Residual Compressive Stress Increase | Minimal change | Decrease | Subsurface crack initiation is suppressed, extending life. |
| Contact Load Increase | Zone expands and intensifies | Sharp Increase (can exceed 1.0) | Dramatic acceleration of fatigue failure. |
Conclusion and Future Perspectives
This study has developed and demonstrated a comprehensive, coupled predictive model for assessing contact fatigue risk in carburized and quenched cylindrical gears. By integrating a discrete numerical method for calculating the multiaxial contact stress field with mathematical representations of hardness and residual stress gradients, and applying the Dang Van multiaxial fatigue criterion in a spatially resolved manner, the model provides a powerful tool for moving beyond homogeneous-material assumptions.
The key findings are: 1) The fatigue risk hotspot is not defined by a single stress component but by the complex interaction between the time-varying multiaxial stress state and the local, gradient-defined material strength. 2) Compressive residual stresses induced by carburizing and quenching play a critical role in reducing fatigue risk, quantified by the model as a significant drop in the $R_{FP}$ index. 3) Parametric studies clearly elucidate failure mode shifts: good lubrication promotes subsurface-origin pitting, while high friction shifts risk to the surface, favoring micropitting; increased surface hardness uniformly improves resistance, and higher loads exponentially accelerate failure risk.
The predictions from this model for the cylindrical gear show excellent correlation with empirical observations from gear testing, such as the distinct crack initiation sites for pitting and micropitting and the load-life relationship. Future work will focus on validating the model against extensive experimental S-N data for different cylindrical gear geometries and material grades, incorporating the effects of non-metallic inclusions, and extending the framework to model the competition between fatigue and wear in real operating conditions. This model establishes a solid foundation for the virtual prototyping and optimized design of high-performance, durable cylindrical gear transmissions.