In modern mechanical power transmission systems, the demand for higher power density and operational reliability places immense stress on critical components like gears. Among the prevalent failure modes, contact fatigue—manifesting as pitting, micro-pitting, and spalling—remains a primary bottleneck limiting the lifespan and performance of gear drives. For cylindrical gears subjected to surface-hardening treatments like carburizing and quenching, the failure mechanism is intricately linked to the complex interplay between the applied multiaxial contact stresses and the material’s inherent graded properties. The conventional design approach, as codified in standards like ISO 6336 and GB/T 3480, primarily relies on Hertzian contact theory to calculate nominal surface stress against a uniform material strength. However, this method falls short in capturing the significant influence of the gradient in hardness and residual stress introduced by the carburizing process. These gradients fundamentally alter the material’s load-bearing capacity through the depth, making subsurface-initiated fatigue a dominant risk. Therefore, developing a predictive model that deeply couples the mechanical stress response from gear meshing with the physically accurate representation of the material’s gradient strength is paramount for accurate fatigue life assessment and design optimization of high-performance cylindrical gears.
This study aims to establish a coupled mathematical model for predicting the contact fatigue risk of carburized and quenched spur (cylindrical) gears. The model integrates the explicit analytical solution for contact stresses with a discrete numerical framework to efficiently solve the singular integral problems inherent in elastic half-plane theory. This computational core is then coupled with a physically-based representation of the hardness and residual stress gradients. Fatigue risk is evaluated using the Dang Van multiaxial fatigue criterion, which is particularly suited for high-cycle rolling contact fatigue where both shear stress amplitude and hydrostatic pressure are critical. The final model provides a spatial map of the material’s load capacity utilization ratio, identifying critical risk zones for crack initiation. The influence of key operational and material parameters—such as friction coefficient, surface hardness, contact load, and the gradient characteristics themselves—is systematically investigated to elucidate their contribution to the fatigue risk of the treated cylindrical gear.
Theoretical Foundation: Gear Meshing and Contact Stress Analysis
Equivalent Cylindrical Gear Meshing Model
The contact between two involute teeth of a spur cylindrical gear pair can be accurately modeled locally as the contact between two equivalent cylinders. Their radii are equal to the radii of curvature of the gear tooth profiles at the specific point of contact along the line of action. For a gear pair with a contact ratio between 1 and 2, the contact point moves from the root to the tip of the driving gear, passing through the pitch point. The equivalent radius of curvature ($\rho$) and the subsequent Hertzian contact parameters vary continuously during this motion. At any meshing position defined by the operating pressure angle $\alpha_{op}$, the geometrical and load parameters can be derived.
The curvature radii for the pinion (p) and gear (q) at the contact point are:
$$
\rho_p = \frac{m z_p \cos \alpha}{2} \tan \alpha_{op}, \quad \rho_g = \frac{m z_g \cos \alpha}{2} \tan (\alpha_t – \alpha_{op})
$$
where $m$ is the module, $z_p$ and $z_g$ are the number of teeth, $\alpha$ is the standard pressure angle, and $\alpha_t$ is the transverse pressure angle at the operating pitch circle. The normal load per unit face width $F_N/B$ is calculated from the transmitted torque. According to Hertzian theory, the half-width of the contact ellipse (for spur gears, a rectangular line contact) $b$ is:
$$
b = \sqrt{ \frac{4 F_N}{\pi B} \cdot \frac{\rho_p \rho_g}{\rho_p + \rho_g} \cdot \left( \frac{1-\nu_p^2}{E_p} + \frac{1-\nu_g^2}{E_g} \right) }
$$
where $E$ and $\nu$ are the Young’s modulus and Poisson’s ratio, respectively. The maximum Hertzian contact pressure $p_0$ is:
$$
p_0 = \sqrt{ \frac{F_N / B}{\pi \cdot \left( \frac{1-\nu_p^2}{E_p} + \frac{1-\nu_g^2}{E_g} \right)} \cdot \frac{\rho_p + \rho_g}{\rho_p \rho_g} }
$$
For the purpose of analyzing a critical location, such as near the pitch point where pure rolling is assumed (though sliding exists in practice due to deformations), the parameters can be calculated. The following table summarizes the key parameters for a sample case used in subsequent analysis.
| Parameter | Symbol | Value |
|---|---|---|
| Module | $m$ | 6.5 mm |
| Pinion Teeth | $z_p$ | 24 |
| Gear Teeth | $z_g$ | 25 |
| Pressure Angle | $\alpha$ | 20° |
| Face Width | $B$ | 30 mm |
| Young’s Modulus | $E$ | 210 GPa |
| Poisson’s Ratio | $\nu$ | 0.3 |
| Applied Contact Pressure | $p_0$ | 1.5 GPa |
| Coefficient of Friction | $f$ | 0.05 |
Elastic Half-Plane Contact Model with Tangential Traction
The stress field within the body of a cylindrical gear tooth subjected to combined normal and tangential surface loads is derived from the theory of elasticity for a half-plane. Consider a normal pressure distribution $p(t)$ and a tangential traction distribution $q(t)$ acting over a surface interval from $x_1$ to $x_2$. 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} \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 \\
\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} \frac{q(t)(x-t)}{[(x-t)^2 + z^2]^2} dt \\
\tau_{xz}(x,z) &= -\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}
$$
For a simple Hertzian pressure distribution $p(x) = p_0 \sqrt{1 – (x/b)^2}$ and a related tangential traction $q(x) = f \cdot p(x)$, these integrals are complex to solve analytically in closed form for the entire field. A discrete numerical approach is therefore employed for efficient and flexible computation.
Discrete Numerical Method via Superposition of Rectangular Load Elements
The contact pressure and traction distributions are discretized into a series of $n$ rectangular load elements of constant width $\Delta x$. Each element $j$ carries a constant normal load $p_j$ and tangential load $q_j$. The stress field induced by a single rectangular normal load element of half-width $s$ and magnitude $p_j$, centered at the origin of a local coordinate system, has an explicit analytical solution:
$$
\begin{aligned}
\sigma_x^{p}(x,z) &= -\frac{p_j}{\pi} \left[ \arctan\left(\frac{s-x}{z}\right) + \arctan\left(\frac{s+x}{z}\right) – \frac{z(s+x)}{(s-x)^2+z^2} + \frac{z(x-s)}{(s+x)^2+z^2} \right] \\
\sigma_z^{p}(x,z) &= -\frac{p_j}{\pi} \left[ \arctan\left(\frac{s-x}{z}\right) + \arctan\left(\frac{s+x}{z}\right) + \frac{z(s+x)}{(s-x)^2+z^2} – \frac{z(x-s)}{(s+x)^2+z^2} \right] \\
\tau_{xz}^{p}(x,z) &= -\frac{p_j z^2}{\pi} \left( \frac{1}{(s-x)^2+z^2} – \frac{1}{(s+x)^2+z^2} \right)
\end{aligned}
$$
Similarly, the stresses due to a constant tangential load element $q_j$ are:
$$
\begin{aligned}
\sigma_x^{q}(x,z) &= -\frac{q_j}{\pi} \left( \frac{z^2}{(s+x)^2+z^2} – \frac{z^2}{(s-x)^2+z^2} + 2 \ln \frac{r_2}{r_1} \right) \\
\sigma_z^{q}(x,z) &= -\frac{q_j z^2}{\pi} \left( \frac{1}{(s-x)^2+z^2} – \frac{1}{(s+x)^2+z^2} \right) \\
\tau_{xz}^{q}(x,z) &= -\frac{q_j}{\pi} \left[ \arctan\left(\frac{s-x}{z}\right) + \arctan\left(\frac{s+x}{z}\right) – \frac{z(s+x)}{(s-x)^2+z^2} + \frac{z(x-s)}{(s+x)^2+z^2} \right]
\end{aligned}
$$
where $r_1^2 = (s-x)^2 + z^2$ and $r_2^2 = (s+x)^2 + z^2$. The total stress state at any subsurface point $i$ is obtained by superimposing the contributions from all load elements $j$, after transforming the local coordinates of each element to the global coordinate system. This superposition can be expressed efficiently using influence coefficient matrices. For the $x$-direction stress:
$$
\sigma_x^i = \sum_{j=1}^{n} \left[ T_{p, j \to i} \cdot p_j + T_{q, j \to i} \cdot q_j \right]
$$
Here, $T_{p, j \to i}$ and $T_{q, j \to i}$ are the pre-computed influence coefficients linking the load on element $j$ to the stress at point $i$, derived from the explicit formulas above. This method allows for rapid calculation of the complete multiaxial stress tensor $\boldsymbol{\sigma}(x,z)$ for any arbitrary discretized surface loading, which is essential for the subsequent fatigue analysis.
Characterization of Material Property Gradients
Hardness Gradient Model
The carburizing and quenching process enriches the surface of the cylindrical gear with carbon, leading to a hardened case with a distinct gradient in hardness from the surface to the core. This gradient is a direct consequence of the carbon concentration profile. The hardness $HV(z)$ as a function of depth $z$ can be effectively modeled. A piecewise quadratic function provides a good fit to empirical data and common prediction formulas (e.g., Thomas’s formula):
$$
HV(z) =
\begin{cases}
a_a z^2 + b_a z + c_a, & 0 \leq z < \text{D}_{CHD} \\
a_b z^2 + b_b z + c_b, & \text{D}_{CHD} \leq z < z_{core} \\
HV_{core}, & z_{core} \leq z
\end{cases}
$$
where $\text{D}_{CHD}$ is the effective case hardening depth (e.g., depth to 550 HV), $z_{core}$ is the transition depth to the homogeneous core, $HV_{core}$ is the core hardness, and the coefficients $a, b, c$ are determined to ensure continuity and match the specified surface hardness $HV_{surface}$ at $z=0$. The hardness gradient is a primary indicator of the material’s local yield and fatigue strength.
Residual Stress Gradient Model
The heat treatment also induces significant residual stresses, predominantly compressive at the surface, which are highly beneficial for retarding fatigue crack initiation. A phenomenological model based on a logistic function effectively captures the typical residual stress profile $\sigma_{RS}(z)$:
$$
\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 asymptotic residual stress in the core (often tensile or near zero), $k$ controls the gradient steepness, and $\delta$ defines the depth of the maximum compression. This model fits well with experimental data from cylindrical gear specimens and allows for parametric studies.
Conversion to Fatigue Strength Parameters
To employ a fatigue criterion, the local hardness must be converted into relevant material strength parameters. For high-cycle fatigue under multiaxial stresses, the fully reversed bending fatigue limit $\sigma_{-1}$ and the fully reversed torsional fatigue limit $\tau_{-1}$ are required. Empirical relationships established in standards like ISO 6336-5 provide a basis for conversion. A linear approximation between Vickers hardness (HV) and ultimate tensile strength ($R_m$) is commonly used:
$$
R_m(z) \approx \alpha \cdot HV(z) + \beta \quad \text{(in MPa)}
$$
Subsequently, the fatigue limits can be estimated as fractions of the tensile strength:
$$
\sigma_{-1}(z) \approx 0.5 \cdot R_m(z), \quad \tau_{-1}(z) \approx 0.577 \cdot \sigma_{-1}(z) \approx 0.2885 \cdot R_m(z)
$$
These relationships establish the depth-dependent material strength field $\lambda(z)$ used in the Dang Van criterion, thereby creating a direct link between the metallurgical gradient and the mechanical fatigue resistance of the cylindrical gear.
Multiaxial Fatigue Risk Assessment via Dang Van Criterion
The Dang Van criterion is a stress-invariant based, multiaxial high-cycle fatigue criterion well-suited for rolling/sliding contact problems. It postulates that fatigue damage occurs in critical shear planes within grains when a linear combination of the local shear stress amplitude and the hydrostatic stress exceeds a material-dependent threshold. For a given material point undergoing a time-varying stress tensor $\boldsymbol{\sigma}(t)$, the criterion is evaluated as follows.
First, the time-varying hydrostatic stress $\sigma_H(t)$ is calculated:
$$
\sigma_H(t) = \frac{1}{3} \text{tr}(\boldsymbol{\sigma}(t)) = \frac{\sigma_{xx}(t) + \sigma_{yy}(t) + \sigma_{zz}(t)}{3}
$$
Second, for each potential material plane defined by its normal vector, the shear stress vector $\vec{\tau}(t)$ is determined. The shear stress amplitude $\Delta \tau /2$ for that plane is half the range of the magnitude of $\vec{\tau}(t)$ over a load cycle. The Dang Van fatigue parameter $\text{DV}(t)$ for that plane at time $t$ is:
$$
\text{DV}(t) = \frac{\Delta \tau}{2} + \alpha_{DV} \cdot \sigma_H(t)
$$
The material constants $\alpha_{DV}$ and the shear fatigue limit $\tau_{DV}$ (often related to $\tau_{-1}$) are determined from the two basic fatigue limits:
$$
\alpha_{DV} = \frac{3 \cdot (\tau_{-1} – \sigma_{-1}/2)}{\sigma_{-1}}, \quad \tau_{DV} = \tau_{-1}
$$
The critical plane is the one where the maximum value of $\text{DV}(t)$ over time, denoted $\text{DV}_{max}$, is found. Fatigue risk is then quantified by a Risk Factor ($RF$), also called the Load Capacity Utilization Factor:
$$
RF = \frac{\text{DV}_{max}}{\tau_{DV}}
$$
A value of $RF \geq 1$ indicates a high probability of fatigue crack initiation. To incorporate the residual stress field $\sigma_{RS}(z)$, it is treated as a static stress component added to the stress tensor from the external load. Therefore, the hydrostatic stress in the criterion becomes $\sigma_H^{total}(t) = \sigma_H^{external}(t) + \sigma_{RS}(z)/3$. The final, gradient-aware risk factor for the cylindrical gear material point is:
$$
RF'(x,z) = \frac{ \max_t \left[ \frac{\Delta \tau(x,z,t)}{2} + \alpha_{DV}(z) \cdot \left( \sigma_H^{ext}(x,z,t) + \frac{\sigma_{RS}(z)}{3} \right) \right] }{\tau_{DV}(z)}
$$
This formulation allows for a spatially resolved (2D map in the x-z plane) prediction of the most vulnerable subsurface regions within the cylindrical gear tooth flank.
Results and Discussion: Model Predictions for Cylindrical Gears
Subsurface Stress Fields
Applying the discrete numerical model to the sample cylindrical gear parameters with a Hertzian pressure of 1.5 GPa and a friction coefficient of 0.05 generates the complete subsurface stress field. The principal stresses $\sigma_1, \sigma_2, \sigma_3$ (where $\sigma_1 \geq \sigma_2 \geq \sigma_3$) are derived from the calculated stress tensor. All three principal stresses are compressive near the surface, with $\sigma_3$ (the most negative) having the greatest magnitude directly under the load center. The evolution of the maximum shear stress $\tau_{max} = (\sigma_1 – \sigma_3)/2$ is particularly instructive. Its contour plot reveals that the maximum value does not occur at the surface but at a subsurface depth approximately 0.3-0.4 mm below the surface, aligning with the classical theory for rolling contact. The presence of tangential traction (friction) slightly shifts the location and magnitude of this maximum. Plotting these stresses along the centerline ($x=0$) versus depth clearly shows the peak in shear stress, identifying the classical orthogonal shear stress reversal zone as a potential site for subsurface-initiated pitting in the cylindrical gear.
Effect of Material Gradients on Fatigue Thresholds
The incorporation of the hardness and residual stress gradients profoundly alters the fatigue risk map compared to a homogeneous material assumption. First, the material shear strength parameter $\tau_{DV}(z)$ decreases with depth following the hardness gradient. The hydrostatic stress $\sigma_H(t)$, which is compressive at the surface, combines with the compressive residual stress $\sigma_{RS}(z)$ to create a more beneficial (more compressive) total hydrostatic state in the case region. For the baseline case, the calculation of $RF'(x,z)$ without residual stress shows a maximum risk factor of approximately 0.83 located at a subsurface point offset from the centerline ($x \approx -0.1$ mm, $z \approx 0.4$ mm). This offset is induced by the sliding friction component. When the residual stress profile is added, the maximum $RF’$ value drops significantly to about 0.71 at a similar location. This demonstrates the potent beneficial effect of compressive residual stresses in raising the effective fatigue limit of the carburized cylindrical gear layer, effectively “shielding” the material from the applied shear stresses. The gradient in strength ($\tau_{DV}(z)$) ensures that even as stresses decay with depth, the material’s resistance decays as well, often keeping the critical risk zone in the subsurface region rather than at the very surface or much deeper core.
Parametric Sensitivity Analysis
The coupled model is an excellent tool for investigating the influence of key parameters on the contact fatigue risk of cylindrical gears.
1. Friction Coefficient: Increasing the coefficient of friction $f$ from 0 (pure rolling) to 0.2 (high sliding) dramatically changes the risk profile. The table below summarizes the findings:
| Friction Coeff. ($f$) | Max RF’ Location ($x$, $z$) | Max RF’ Value | Primary Implication |
|---|---|---|---|
| 0.00 | (0.0, ~0.45 mm) | 0.65 | Risk is symmetric, deepest. |
| 0.05 (Baseline) | (-0.1, ~0.40 mm) | 0.71 | Risk shifts slightly upstream of rolling. |
| 0.10 | (-0.2, ~0.25 mm) | 0.82 | Risk moves towards surface and further upstream. |
| 0.20 | (-0.3, ~0.15 mm) | 1.05 | Risk zone migrates near-surface, RF>1 indicates high failure probability. |
This migration is critical: under good lubrication (low $f$), the critical shear stress is subsurface, favoring classic pitting. Under poor lubrication (high $f$), the surface and near-surface shear stresses increase markedly, making surface-originated micro-pitting and scuffing more likely failure modes for the cylindrical gear.
2. Surface Hardness: Increasing the surface hardness $HV_{surface}$ from 600 to 750 HV, while keeping the case depth constant, directly increases the material strength parameter $\tau_{DV}(z)$ in the near-surface region. The maximum $RF’$ value decreases proportionally (e.g., from 0.78 to 0.65 for a given load), but the spatial location of the risk zone remains largely unchanged. This confirms that surface hardening is an effective means to increase the safety margin against contact fatigue, but it does not alter the fundamental mechanism of where failure is most likely to initiate in the cylindrical gear.
3. Applied Contact Pressure (Load): Increasing the maximum Hertzian pressure $p_0$ from 1.5 GPa to 1.9 GPa has a nonlinear and severe impact. Not only does the contact half-width $b$ increase, but the subsurface stress magnitudes scale proportionally. For the higher load, the computed $RF’$ maximum exceeds 1.0, clearly predicting a high probability of fatigue failure within a finite number of cycles. This aligns perfectly with the standard $S-N$ curve behavior and gear endurance testing, where increased load drastically reduces the expected life of the cylindrical gear.
4. Case Depth and Gradient Steepness: Varying the effective case depth $\text{D}_{CHD}$ changes the gradient profile. A shallower case (e.g., 1.0 mm) results in a steeper hardness drop-off. This can cause the location of maximum $RF’$ to coincide with a region where the strength $\tau_{DV}(z)$ is decreasing rapidly, potentially creating a sharper risk peak. A deeper, more gradual case (e.g., 2.5 mm) provides a thicker “armor” of high-strength material, typically lowering and broadening the risk distribution, offering more damage tolerance to the cylindrical gear.
Conclusion
This study has developed and demonstrated a comprehensive, physics-based predictive model for assessing the contact fatigue risk in carburized and quenched spur cylindrical gears. The model successfully couples a high-fidelity elastic contact stress analysis, solved via an efficient discrete numerical method, with a realistic characterization of the material’s graded properties (hardness and residual stress). The application of the Dang Van multiaxial fatigue criterion translates this coupled stress-strength state into a spatially resolved risk factor ($RF’$), providing a powerful tool for identifying potential failure initiation sites—be they subsurface or near-surface.
The key findings are: First, the model accurately captures the transition in failure mode with lubrication condition; low friction keeps the critical risk zone subsurface (classic pitting), while high friction migrates it towards the surface (micro-pitting). Second, increasing surface hardness directly enhances the safety margin by elevating the material’s fatigue strength gradient. Third, compressive residual stresses induced by carburizing and quenching provide a significant and quantifiable benefit by superimposing a favorable hydrostatic stress, effectively lowering the calculated risk factor. Fourth, increased applied load predictably escalates the risk factor, potentially driving it beyond the failure threshold. The predictions from this model show strong qualitative and quantitative agreement with established empirical knowledge and experimental observations regarding the fatigue behavior of cylindrical gears.
This approach moves beyond the limitations of standardized nominal stress methods by explicitly accounting for the gradient nature of high-performance gear materials. It provides a valuable framework for the virtual prototyping and optimization of cylindrical gears, enabling engineers to tailor heat treatment parameters (case depth, surface hardness) and operational conditions (load, lubrication) to achieve target reliability and life goals, thereby addressing a critical bottleneck in the design of advanced power transmission systems.