The operational reliability and longevity of power transmission systems are fundamentally dependent on the performance of their gear components. Among the various failure modes, contact fatigue, manifesting as micropitting, pitting, and spalling, remains a critical bottleneck limiting the power density and service life of modern machinery. Standard calculation methods, such as those outlined in ISO 6336 and GB/T 3480, provide a foundational framework based on nominal contact stress and material strength. However, these approaches often treat material properties as homogeneous, overlooking the profound influence of gradient characteristics introduced by surface-hardening processes like carburizing and quenching. For high-performance cylindrical gears, the subsurface region possesses a complex gradient in hardness and residual stress, which fundamentally alters the stress response and fatigue crack initiation mechanisms under cyclic contact loading. Therefore, a predictive model that deeply integrates the gradient mechanical properties of the case-hardened layer with the multiaxial stress state induced by gear meshing is essential for accurate fatigue life assessment and robust design.
This study aims to establish a coupled mathematical model for predicting the contact fatigue risk of case-hardened spur cylindrical gears. The model explicitly incorporates gradient representations of hardness and residual stress, interfaces them with a detailed elastic contact stress analysis, and evaluates fatigue risk using a multiaxial fatigue criterion. We begin by establishing the gear meshing contact model, followed by the characterization of property gradients. Subsequently, the Dang Van multiaxial fatigue criterion is integrated to calculate the material load utilization factor. Finally, a comprehensive parametric analysis is conducted to elucidate the influence of key operational and material parameters on the predicted fatigue risk.
1. Gear Meshing and Elastic Contact Stress Analysis
1.1 Equivalent Contact Model for Spur Cylindrical Gears
For a pair of spur cylindrical gears with a contact ratio \(\epsilon\) between 1 and 2, the meshing process alternates between single and double tooth pairs along the line of action. At any contact point during meshing, the local geometry can be accurately represented by the contact of two equivalent cylinders. The principal radii of curvature at the contact point are derived from the basic parameters of the involute profile. According to the Hertzian theory for elastic contact of cylinders, the semi-contact width \(b\) and the nominal contact parameters can be expressed as a function of the operating conditions and gear geometry.
The load per unit face width \(F\) is given by:
$$ F = \frac{9,549 P_0}{n R_{pm} \cos \alpha_{pm}} $$
where \(P_0\) is the transmitted power, \(n\) is the rotational speed (rpm), \(R_{pm}\) is the radius at the contact point on the pinion, and \(\alpha_{pm}\) is the operating pressure angle at that point.
The contact semi-width \(b\) is calculated as:
$$ 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) } $$
where \(B\) is the face width, \(\rho_{pm}\) and \(\rho_{qm}\) are the radii of curvature for the pinion and gear, \(\nu\) is Poisson’s ratio, and \(E\) is the Young’s modulus. The subscripts \(p\) and \(q\) denote the pinion and gear, respectively. The radii of curvature are:
$$ \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 \(m\) is the module, \(z\) is the number of teeth, and \(\alpha\) is the standard pressure angle.
1.2 Elastic Half-Plane Contact Formulation
The subsurface stress field resulting from surface tractions on an elastic half-plane forms the basis for contact mechanics analysis. For a general distribution of normal pressure \(p(t)\) and tangential traction \(q(t)\) over a surface region from \(x_1\) to \(x_2\), the stress components at any subsurface point \((x, z)\) are given by the following integral equations, where the \(x\)-axis lies on the surface and the \(z\)-axis points into the material:
$$
\begin{aligned}
\sigma_x &= -\frac{2z}{\pi} \int_{x_1}^{x_2} \frac{p(t) (x-t)^2}{[(x-t)^2 + z^2]^2} dt – \frac{2}{\pi} \int_{x_1}^{x_2} \frac{q(t) (x-t)^3}{[(x-t)^2 + z^2]^2} dt \\[6pt]
\sigma_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} \frac{q(t) (x-t)}{[(x-t)^2 + z^2]^2} dt \\[6pt]
\tau_{xz} &= -\frac{2z^2}{\pi} \int_{x_1}^{x_2} \frac{p(t) (x-t)}{[(x-t)^2 + z^2]^2} dt – \frac{2z}{\pi} \int_{x_1}^{x_2} \frac{q(t) (x-t)^2}{[(x-t)^2 + z^2]^2} dt
\end{aligned}
$$
These equations represent singular integrals, and analytical solutions are only available for simple load distributions. For complex load patterns, such as those from real gear contact with friction, numerical methods are indispensable.
1.3 Discrete Numerical Solution Method
A highly efficient numerical approach involves discretizing the contact interface 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 rectangular uniform load can be solved analytically. By superimposing the contributions from all discrete elements, the total stress at any evaluation point \(i\) can be computed rapidly.
For instance, the stress in the \(x\)-direction at point \(i\) due to all elements is:
$$ \sigma_x^i = \sum_{j=1}^{n} \left( T_{p,j-i} \cdot p_j + T_{q,j-i} \cdot q_j \right) $$
where \(T_{p,j-i}\) and \(T_{q,j-i}\) are influence functions derived from the analytical solution for a rectangular load patch. These functions depend only on the relative geometry between the load element \(j\) and the evaluation point \(i\), making the computation extremely efficient. Similar expressions are derived for \(\sigma_z^i\) and \(\tau_{xz}^i\).
This method allows for the precise calculation of the complete multiaxial subsurface stress tensor for any given surface load distribution, including the combined effect of normal Hertzian pressure and tangential friction traction, which is crucial for analyzing cylindrical gears under rolling-sliding contact.
2. Characterization of Gradient Properties in Case-Hardened Cylindrical Gears
2.1 Hardness Gradient Model
The carburizing and quenching process generates a hardened case on the surface of cylindrical gears. The carbon concentration gradient leads to a corresponding gradient in hardness, which is the primary contributor to the enhanced load-carrying capacity. The hardness profile, typically measured as Vickers hardness (HV) versus depth from the surface, can be characterized by key parameters: surface hardness \(HV_{surf}\), core hardness \(HV_{core}\), and the effective case hardening depth \(D_{CHD}\) (e.g., depth to 550 HV).
An empirical model, such as the Thomas formula, effectively captures this gradient. The relationship can be expressed as a piecewise function, often approximated by quadratic polynomials in different zones. For computational implementation within our predictive model, we utilize such a formula to generate the hardness gradient \(HV(z)\):
$$ 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}
$$
The coefficients \(a, b, c\) are determined by fitting the boundary conditions (surface hardness, core hardness, and continuity at \(D_{CHD}\)).
The following table summarizes the effect of different case depths on the gradient steepness for a typical gear steel (e.g., 18CrNiMo7-6).
| Effective Case Depth, \(D_{CHD}\) (mm) | Surface Hardness, \(HV_{surf}\) | Core Hardness, \(HV_{core}\) | Gradient Steepness Description |
|---|---|---|---|
| 1.2 | 700 | 450 | Very steep transition |
| 1.8 | 700 | 450 | Moderately steep |
| 2.4 | 700 | 450 | Gradual transition |
2.2 Residual Stress Gradient Model
Carburizing and quenching also induces a characteristic residual stress profile, typically compressive at the surface, transitioning to tensile in the core, and then back to near-zero. This residual stress field significantly affects fatigue performance by superimposing on the applied contact stresses. A statistical prediction model, such as the one developed by FZG, provides a reliable sigmoidal function to describe the gradient \(\sigma_{RS}(z)\):
$$ \sigma_{RS}(z) = \sigma_D + \frac{\sigma_z – \sigma_D}{1 + e^{-k(z + \delta)}} $$
where \(\sigma_D\) is the maximum compressive residual stress (negative value), \(\sigma_z\) is the maximum tensile residual stress (positive value), and \(k\) and \(\delta\) are fitting parameters controlling the shape and inflection point of the curve. This model is integrated into our fatigue risk calculation.
2.3 Strength Parameter Conversion
To evaluate material strength under multiaxial fatigue criteria, parameters like the fully reversed bending fatigue limit \(\sigma_{-1}\) and the torsional fatigue limit \(\tau_{-1}\) are required. As direct measurement of these gradients is impractical, they are derived from the hardness gradient. Empirical relationships established in standards like ISO 6336-5 provide conversion. A linear fit is often sufficiently accurate for high-strength steels:
$$ \sigma_{UTS}(z) \approx \alpha \cdot HV(z) + \beta $$
where \(\sigma_{UTS}(z)\) is the ultimate tensile strength gradient, and \(\alpha\) and \(\beta\) are material-specific constants. The fatigue limits can then be estimated as a fraction of \(\sigma_{UTS}(z)\) (e.g., \(\sigma_{-1} \approx 0.45 \sigma_{UTS}\) for bending, \(\tau_{-1} \approx 0.577 \sigma_{-1}\) for shear based on the von Mises relation). This establishes the crucial link between the measurable hardness gradient and the material strength parameters needed for fatigue assessment in cylindrical gears.
3. Dang Van Multiaxial Fatigue Criterion for Risk Prediction
The Dang Van criterion is a high-cycle fatigue (HIF) criterion based on the concept of a critical plane where crack initiation is likely to occur. It postulates that failure happens when a linear combination of the shear stress amplitude and the hydrostatic stress on a material plane exceeds a critical value related to the material’s fatigue limits. This criterion is particularly well-suited for rolling-sliding contact in cylindrical gears because it accounts for the non-proportional, multiaxial stress history and the significant role of hydrostatic pressure.
For a given stress history at a material point, the criterion is evaluated on numerous material planes. For each plane with normal vector defined by angle \(\alpha\) (in 2D analysis), the shear stress amplitude \(\Delta \tau_{max}(\alpha, t)\) over a load cycle and the instantaneous hydrostatic stress \(\sigma_H(t)\) are calculated. The Dang Van damage parameter is:
$$ \text{RFP}(\alpha, t) = \frac{\Delta \tau_{max}(\alpha, t) + \kappa_D \cdot \sigma_H(t)}{\lambda} $$
where:
- \(\text{RFP}\) is the Risk Factor or material load utilization factor. An RFP ≥ 1 indicates a high risk of fatigue crack initiation.
- \(\kappa_D\) and \(\lambda\) are material parameters derived from the fully reversed bending and torsion fatigue limits: \(\kappa_D = \frac{3\tau_{-1} – \sigma_{-1}}{\sigma_{-1}}\) and \(\lambda = \tau_{-1}\).
The critical plane is identified as the one where the maximum value of \(\text{RFP}(\alpha, t)\) over time \(t\) and angle \(\alpha\) is found. This maximum value, \(\text{RFP}_{max}\), is the fatigue risk indicator for that material point.
Integration with Gradient Properties: To account for the case-hardening of cylindrical gears, the criterion is modified. The material parameters \(\kappa_D\) and \(\lambda\) become functions of depth \(z\) through their dependence on \(\sigma_{-1}(z)\) and \(\tau_{-1}(z)\). Furthermore, the residual stress \(\sigma_{RS}(z)\) is introduced as a mean stress component, effectively modifying the hydrostatic stress term. The enhanced criterion becomes:
$$ \text{RFP}'(\alpha, z, t) = \frac{\Delta \tau_{max}(\alpha, z, t) + \kappa_D(z) \cdot [\sigma_H(z, t) + \sigma_{RS}(z)]}{\lambda(z)} $$
This formulation allows the fatigue risk prediction model to reflect the spatially varying material strength and the beneficial effect of compressive residual stresses throughout the case-hardened layer of the cylindrical gears.
4. Model Implementation and Parametric Influence Analysis
We implement the coupled model for a standard spur gear pair. The key parameters are listed below.
| 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 |
| Young’s Modulus | \(E\) | 210 GPa |
| Poisson’s Ratio | \(\nu\) | 0.3 |
| Core Hardness | \(HV_{core}\) | 450 HV |
| Effective Case Depth | \(D_{CHD}\) | 2.1 mm |
The analysis focuses on the pitch point contact, with a nominal contact pressure of 1500 MPa and a friction coefficient of \(\mu = 0.05\) to simulate lubricated rolling-sliding conditions.
4.1 Subsurface Stress State and Fatigue Risk Maps
The discrete numerical method calculates the complete stress tensor field. The principal stresses (\(\sigma_1, \sigma_2, \sigma_3\)) and the maximum shear stress \(\tau_{max}\) are derived. As expected, high compressive stresses dominate near the surface, while the maximum shear stress peaks in the subsurface region (around 0.3 mm depth for this load case).
Applying the gradient-enhanced Dang Van criterion yields a spatial map of the maximum material load utilization factor \(\text{RFP}’_{max}(x,z)\). The analysis reveals a critical insight: the location of maximum risk (\(\text{RFP}’_{max} \approx 0.71\)) is not at the surface nor directly under the contact center. It is located in the subsurface region, offset laterally from the center due to the rolling-sliding action, and at a depth influenced by the competing effects of decreasing applied shear stress and decreasing material strength (\(\lambda(z)\)). This predicted risk zone aligns with the typical subsurface origin of classical pitting fatigue in cylindrical gears.
4.2 Influence of Friction Coefficient
Friction, representing the tangential traction, significantly alters the risk profile. The following table and analysis summarize its effect.
| Friction Coeff. (\(\mu\)) | Max RFP’ Value | Risk Zone Location (Depth) | Primary Implication |
|---|---|---|---|
| 0.00 (Pure Roll) | 0.68 | Deep Subsurface (~0.4 mm) | Classical subsurface pitting risk. |
| 0.05 | 0.71 | Mid Subsurface (~0.4 mm) | Increased risk, slight upward shift. |
| 0.10 | 0.78 | Shallow Subsurface (~0.25 mm) | Higher risk, zone migrates towards surface. |
| 0.20 | > 1.0 | Near-Surface / Surface | Very high risk for surface-initiated micropitting. |
The analysis clearly demonstrates that increasing friction (simulating lubricant breakdown or poor lubrication) not only increases the magnitude of the fatigue risk factor but also causes the critical risk zone to migrate from the deep subsurface towards the surface. This provides a mechanistic explanation for the transition from subsurface-origin pitting (under good lubrication) to surface-origin micropitting (under poor lubrication) commonly observed in cylindrical gears.
4.3 Influence of Surface Hardness and Load
Surface Hardness: Increasing the surface hardness, representing a more effective case hardening process, directly increases the material strength parameter \(\lambda(z)\) in the near-surface region. This results in a uniform reduction of the RFP’ value across the case depth without changing the location of the risk zone. For example, increasing surface hardness from 600 HV to 750 HV can reduce the max RFP’ by approximately 15-20%, significantly enhancing the contact fatigue resistance of the cylindrical gears.
Applied Load (Contact Pressure): Increasing the transmitted load, and consequently the Hertzian contact pressure, has a dramatic effect. It increases the magnitude of all stress components. Since the material strength gradient remains unchanged, the RFP’ value increases proportionally more. At a critical load level (e.g., ~1900 MPa in this model), the max RFP’ surpasses 1.0, indicating a high probability of fatigue failure. This aligns perfectly with the concept of a bending strength limit and explains the load-dependent finite life behavior of cylindrical gears in contact fatigue testing.
4.4 Influence of Residual Stress
The incorporation of the compressive residual stress gradient \(\sigma_{RS}(z)\) is a key advantage of this model. The residual stress contributes to the hydrostatic stress term \(\sigma_H\) in the Dang Van equation. Since \(\kappa_D\) is typically positive for metals, the compressive (negative) residual stress reduces the value of the numerator \(\Delta \tau_{max} + \kappa_D (\sigma_H + \sigma_{RS})\), thereby reducing the RFP’. In the presented case, the inclusion of residual stress reduced the max RFP’ from 0.83 (without residual stress) to 0.71 (with residual stress). This quantifies the significant beneficial effect of compressive residual stresses in delaying crack initiation in case-hardened cylindrical gears.
5. Conclusion
This study has developed and demonstrated a comprehensive, physics-based predictive model for assessing the contact fatigue risk of case-hardened spur cylindrical gears. The model successfully integrates several advanced aspects: a numerically efficient elastic contact stress solver, empirical models for hardness and residual stress gradients, and the Dang Van multiaxial fatigue criterion modified for gradient properties. The primary conclusions are as follows:
- The model provides a fast, deterministic method for predicting the location and severity of contact fatigue risk in spur cylindrical gears, moving beyond homogeneous material assumptions.
- The gradient of material properties is crucial. The hardness gradient defines the spatial variation of fatigue strength, while the compressive residual stress gradient provides a significant beneficial effect by reducing the effective driving force for fatigue crack initiation throughout the case depth.
- Parametric studies offer deep mechanistic insights:
- Increasing friction coefficient causes the fatigue risk zone to migrate from the subsurface towards the surface, theoretically explaining the transition from pitting to micropitting failure modes.
- Increasing surface hardness uniformly lowers the fatigue risk factor, enhancing overall durability.
- Increasing applied load can drive the risk factor above the critical threshold, leading to predicted finite-life failure, consistent with experimental S-N curve behavior.
- The predictions regarding failure location (subsurface vs. surface) and the influence of load and lubrication align with well-established empirical observations from gear testing, validating the model’s relevance for practical engineering analysis of cylindrical gears.
This coupled gradient-fatigue risk model serves as a powerful tool for the virtual prototyping and design optimization of high-performance cylindrical gears. It enables engineers to evaluate the impact of heat treatment parameters (case depth, hardness), manufacturing-induced residual stresses, and operating conditions (load, lubrication) on contact fatigue life, thereby facilitating the development of more reliable and power-dense gear transmissions. Future work may extend the model to include microstructural effects, explicit crack propagation modeling, and application to helical and other types of cylindrical gears through appropriate geometric equivalency transformations.