The dynamic performance, vibration, and acoustic emission of gear transmission systems are fundamentally governed by their time-varying mesh stiffness (TVMS). This parameter serves as a primary internal excitation source. Accurate modeling of TVMS is therefore critical for reliable dynamic response prediction, health monitoring, and remaining useful life estimation of gear systems. Among various factors affecting mesh stiffness, tooth surface wear, a gradual material loss due to contact fatigue and sliding friction, progressively alters the gear tooth profile and consequently the load distribution and deflection characteristics. While simplified uniform or linear wear models have been used, they often deviate from the actual non-uniform wear patterns observed in spur gears. This article presents a comprehensive numerical methodology for calculating the mesh stiffness of spur gears that explicitly accounts for realistic, non-uniform tooth surface wear. The approach synergistically integrates the generating method for precise tooth profile definition, the potential energy method for stiffness calculation, and the Archard wear model for wear depth prediction.
The foundation of an accurate stiffness model lies in a precise geometric representation of the gear tooth. The profile of a spur gear manufactured by a rack cutter consists of two distinct segments: the involute profile and the root fillet (or transition curve). The generating method principle is employed to derive parametric equations for these segments using the rolling angle of the imaginary rack cutter as the common variable.
Figure 1 illustrates the coordinate systems. The cutter coordinate system $Px_0y_0$ is fixed to the rack cutter, with its origin at the pitch point $P$. The gear coordinate system $Oxy$ has its origin at the gear center $O$, located at a distance $r = mz/2$ from $P$, where $m$ is the module and $z$ is the number of teeth.
The transition curve, generated by the tip rounding of the rack cutter, is described by:
$$ \begin{cases} x_1(\varphi) = (r – x_c – \rho_0 \cos \gamma)\cos \varphi + (x_c \tan \gamma + \rho_0 \sin \gamma)\sin \varphi \\ y_1(\varphi) = (r – x_c – \rho_0 \cos \gamma)\sin \varphi – (x_c \tan \gamma + \rho_0 \sin \gamma)\cos \varphi \end{cases} $$
where $\rho_0$ is the tip radius of the cutter, $(x_c, y_c)$ are the coordinates of the tip circle center in the cutter system, and $\gamma = \arctan[(r\varphi – y_c)/x_c]$. The valid range for the rolling angle $\varphi$ on the transition curve is $\varphi_E \le \varphi \le \varphi_D$, where $\varphi_E$ is found from the intersection with the root circle of radius $r_f$, and $\varphi_D = (y_c + x_c \cot \alpha)/r$.
The involute profile, generated by the straight-sided flank of the rack cutter, is given by:
$$ \begin{cases} x_2(\varphi) = \left[ r – \frac{1}{2}(r\varphi – y_w)\sin 2\alpha \right] \cos \varphi + (r\varphi – y_w)\cos^2\alpha \sin \varphi \\ y_2(\varphi) = \left[ r – \frac{1}{2}(r\varphi – y_w)\sin 2\alpha \right] \sin \varphi – (r\varphi – y_w)\cos^2\alpha \cos \varphi \end{cases} $$
Here, $y_w = \frac{m}{2}(\pi/2 + 2x \tan \alpha)$ represents half of the tooth space width on the pitch line, $x$ is the profile shift coefficient, and $\alpha$ is the pressure angle (typically 20°). The valid range is $\varphi_B \le \varphi \le \varphi_D$, where $\varphi_B$ is determined from the intersection with the addendum circle of radius $r_a$.
Mesh Stiffness Calculation via the Potential Energy Method
The total mesh stiffness of a pair of spur gears is computed by considering the elastic energy stored in the gear teeth during loading. According to the potential energy method, the total energy $U_{total}$ comprises Hertzian contact energy $U_h$, bending energy $U_b$, shear energy $U_s$, axial compressive energy $U_a$, and the fillet foundation energy $U_f$. For a meshing force $F$, the energy components and their corresponding stiffnesses $k_i$ are related by $U_i = F^2/(2k_i)$.
The Hertzian contact stiffness for two parallel cylinders is constant and given by:
$$ k_h = \frac{\pi E L}{4(1-\nu^2)} $$
where $E$ is Young’s modulus, $L$ is the face width, and $\nu$ is Poisson’s ratio.
The fillet foundation stiffness accounts for the deflection of the gear body below the tooth. A widely accepted formula is:
$$ \frac{1}{k_f} = \frac{\cos^2 \beta_c}{EL} \left\{ L^*\left(\frac{u_f}{S_f}\right)^2 + M^*\left(\frac{u_f}{S_f}\right) + P^*(1 + Q^*\tan^2 \beta_c) \right\} $$
Here, $\beta_c$ is the load angle relative to the tooth centerline, $u_f$ is the distance from the load application point to the root circle along the line of force, and $S_f$ is half the thickness of the root chord. The coefficients $L^*, M^*, P^*, Q^*$ are polynomial functions of the gear’s geometrical parameters.
The bending, shear, and axial compressive stiffnesses are derived by modeling the tooth as a non-uniform cantilever beam. Their inverses are calculated by integrating the strain energy along the tooth profile from the root to the contact point $C$ with coordinates $(x_C, y_C)$:
$$ \frac{1}{k_b} = \int_{\varphi_E}^{\varphi_D} \frac{3[\cos\beta_C (x_C – x_1) – y_C \sin\beta_C]^2}{2EL y_1^3} \frac{dx_1}{d\varphi} d\varphi + \int_{\varphi_C}^{\varphi_D} \frac{-3[\cos\beta_C (x_C – x_2) – y_C \sin\beta_C]^2}{2EL y_2^3} \frac{dx_2}{d\varphi} d\varphi $$
$$ \frac{1}{k_s} = \int_{\varphi_E}^{\varphi_D} \frac{1.2 \cos^2 \beta_C}{2GL y_1} \frac{dx_1}{d\varphi} d\varphi + \int_{\varphi_C}^{\varphi_D} \frac{-1.2 \cos^2 \beta_C}{2GL y_2} \frac{dx_2}{d\varphi} d\varphi $$
$$ \frac{1}{k_a} = \int_{\varphi_E}^{\varphi_D} \frac{\sin^2 \beta_C}{2EL y_1} \frac{dx_1}{d\varphi} d\varphi + \int_{\varphi_C}^{\varphi_D} \frac{-\sin^2 \beta_C}{2EL y_2} \frac{dx_2}{d\varphi} d\varphi $$
In these equations, $G=E/[2(1+\nu)]$ is the shear modulus, and $\beta_C = \varphi_C – \alpha_0$, where $\varphi_C$ is the rolling angle at the contact point and $\alpha_0$ is the standard pressure angle. The integrals are split between the transition curve (using $x_1, y_1$) and the involute curve (using $x_2, y_2$).
The single-tooth-pair mesh stiffness $k_{pair}$ is then obtained by considering all energy components for both the driving and driven spur gears:
$$ k_{pair} = \frac{1}{ \frac{1}{k_h} + \sum_{i=1}^{2}\left( \frac{1}{k_{f,i}} + \frac{1}{k_{b,i}} + \frac{1}{k_{s,i}} + \frac{1}{k_{a,i}} \right) } $$
where $i=1,2$ denotes the driving and driven gear, respectively.
For a gear pair with a contact ratio $\varepsilon$ between 1 and 2, the total mesh stiffness $K(t)$ varies periodically between single and double tooth engagement. It can be represented as a piecewise function over one mesh period $T$:
$$ K(t) = \begin{cases}
k(m(t)) + k(m(t)+T), & 0 \le m(t) < (\varepsilon-1)T \\
k(m(t)), & (\varepsilon-1)T \le m(t) < T
\end{cases} $$
where $m(t) = \mod(t, T)$ is the elapsed time within the current mesh cycle, and $k(\cdot)$ is the single-tooth-pair stiffness function. This function is often approximated by a first-order Fourier series expansion for analytical convenience:
$$ k(t) = k_0 \left( 1 + \epsilon_{at} \cos(\omega_t t) + \epsilon_{bt} \sin(\omega_t t) \right) $$
where $k_0$ is the average stiffness, and $\epsilon_{at}, \epsilon_{bt}$ are harmonic coefficients.
Modeling Non-Uniform Tooth Surface Wear
Tooth wear in spur gears is inherently non-uniform because the contact pressure, sliding velocity, and curvature vary along the path of contact. The Archard wear model provides a fundamental framework for predicting wear volume. The model states that the wear volume $V$ is proportional to the normal load $F_n$, the sliding distance $S$, and a dimensionless wear coefficient $k$, and inversely proportional to the material hardness $H$: $V = k F_n S / H$.
For a specific point on the tooth profile of a driving spur gear, the normal wear depth $h_I$ accumulated over a running time $t$ can be expressed as:
$$ h_I = S \, n_1 t \, I_h \, W_M W_L W_P $$
where $S=2a\lambda_1$ is the total sliding distance per engagement cycle, $a$ is the Hertzian contact half-width, $\lambda_1$ is the specific sliding coefficient for the driving gear, $n_1$ is the rotational speed (rpm), and $I_h = k/H$ is the wear rate. The coefficients $W_M$, $W_L$, and $W_P$ account for surface modification, lubrication condition, and load fluctuation effects, respectively. For standard conditions (no special treatment, basic lubricant, steady load), these factors are set to unity: $W_M=W_L=W_P=1$.
Assuming identical material properties for both mating spur gears, the wear depth at the same pressure angle location is equal. The wear depth is a function of the pressure angle $\alpha_k$ at the meshing point on the involute profile, which can be calculated from the profile coordinates:
$$ \alpha_k = \arccos\left( \frac{r_b}{\sqrt{x_2^2(\varphi_k) + y_2^2(\varphi_k)}} \right) $$
where $r_b = r \cos \alpha_0$ is the base circle radius.
Numerical simulation using the Archard model reveals a characteristic wear pattern: maximum wear occurs near the root and tip of the tooth where sliding is predominant, while theoretically zero wear occurs at the pitch point (pure rolling). The relationship between wear depth $\Delta h$ and pressure angle $\alpha_k$ for a given number of meshing cycles $N$ can be accurately fitted by a cubic polynomial. For a specific case with parameters listed in Table 1 and $N=10^6$ cycles, the fitted relationship is:
$$ \Delta h(\alpha_k) = | A_N \alpha_k^3 + B_N \alpha_k^2 + C_N \alpha_k + D_N | $$
The coefficients $A_N, B_N, C_N, D_N$ are linear functions of the meshing cycles $N$, reflecting the progressive nature of wear.
| Parameter | Driving/Driven Gear |
|---|---|
| Number of Teeth, $z$ | 40 |
| Module, $m$ (mm) | 3 |
| Pressure Angle, $\alpha_0$ (°) | 20 |
| Face Width, $L$ (mm) | 20 |
| Young’s Modulus, $E$ (GPa) | 206 |
| Poisson’s Ratio, $\nu$ | 0.3 |
| Addendum Coefficient, $h_a^*$ | 1.0 |
| Dedendum Coefficient, $c^*$ | 0.25 |
| Contact Ratio, $\varepsilon$ | 1.7135 |
| Coefficient | Expression as a Function of Cycles $N$ |
|---|---|
| $A_N$ | $1.3452 \times 10^{-6} \, N$ |
| $B_N$ | $-1.8301 \times 10^{-6} \, N$ |
| $C_N$ | $9.2300 \times 10^{-7} \, N$ |
| $D_N$ | $-1.5629 \times 10^{-7} \, N$ |
The worn tooth profile is obtained by offsetting the original involute profile inward, normal to the surface at each point, by the calculated wear depth $\Delta h$. The coordinates of the worn profile $(x’_2, y’_2)$ are:
$$ \begin{cases} x’_2 = x_2 – \Delta h \sin \beta_k \\ y’_2 = y_2 – \Delta h \cos \beta_k \end{cases} $$
where $\beta_k$ is the angle of the surface normal at the pre-wear contact point. This modified profile is then used in the potential energy integrals to compute the mesh stiffness of the worn spur gears.
Results, Analysis, and Validation
First, the validity of the proposed modeling approach for healthy (unworn) spur gears is established. The mesh stiffness calculated using the presented energy method, which employs the exact generating profile equations, is compared with results from other established methods. Using the parameters from Table 1, Figure 2 shows the comparison over one complete mesh cycle. The results from the presented method show excellent agreement with those from refined finite element analyses and an improved energy method that models the tooth as a cantilever from the root circle. In contrast, a classical energy method that approximates the entire tooth profile as an involute from the base circle overestimates the stiffness. This discrepancy occurs because when the root circle is larger than the base circle (common for gears with $z > 41$ for $\alpha=20^\circ$), the classical method incorrectly uses a thinner, more flexible tooth section, leading to a lower calculated stiffness—our method correctly uses the more rigid geometry defined by the actual transition curve, yielding higher and more accurate stiffness values. This validates the geometric precision of our model.
The core analysis investigates the impact of non-uniform surface wear on the mesh stiffness of spur gears. The wear depth is calculated for different cumulative meshing cycles $N$, representing progressive stages of gear life. The corresponding time-varying mesh stiffness (TVMS) is then computed. Figure 3 illustrates the TVMS for an unworn gear pair and for pairs with increasing levels of wear ($N = 3\times10^6, 5\times10^6, 7\times10^6$ cycles). A clear trend is observed: as wear progresses, the overall mesh stiffness decreases. The reduction is not uniform across the mesh cycle. The characteristic “double-peak” pattern within the double-tooth-contact region becomes less pronounced with wear.
The reduction in stiffness is quantified by calculating the relative reduction $\Delta K/K_0$ at different angular positions of the driving gear, where $K_0$ is the stiffness of the unworn gear. Figure 4 plots this relative reduction for different wear levels. Key observations include:
- The stiffness reduction is not constant; it varies along the path of contact.
- The reduction is generally more significant in the double-tooth-contact regions compared to the single-tooth-contact region for the same wear level.
- Even for substantial wear ($N=7\times10^6$, with maximum wear depths of ~0.18 mm at the root and ~0.08 mm at the tip), the maximum stiffness reduction is relatively modest, around 1.6%. This indicates that early and moderate wear primarily alters the load-sharing relationship between tooth pairs rather than drastically weakening an individual tooth’s stiffness.
This nuanced understanding is crucial for condition monitoring, as small but systematic changes in mesh stiffness can be an early indicator of surface degradation in spur gears.
| Wear State (Cycles, N) | Max. Wear Depth (mm) | Peak Mesh Stiffness, $K_{max}$ (N/m) | Reduction vs. Unworn |
|---|---|---|---|
| Unworn (0) | 0.000 | 2.018 × 109 | 0.0% |
| Moderate (3×106) | 0.078 | 2.001 × 109 | 0.84% |
| Heavy (5×106) | 0.130 | 1.990 × 109 | 1.39% |
| Severe (7×106) | 0.182 | 1.986 × 109 | 1.59% |
The relationship between wear depth and operating conditions for spur gears can be generalized. The wear rate $I_h$ is central. The normal load $F_n$ at a meshing point is related to the transmitted torque $T$ and the base radius $r_b$:
$$ F_n = \frac{T}{r_b \cos\alpha_k} $$
The contact half-width $a$ for two cylindrical surfaces in line contact is:
$$ a = \sqrt{ \frac{4 F_n}{\pi L} \frac{(1-\nu_1^2)/E_1 + (1-\nu_2^2)/E_2}{(1/R_1) + (1/R_2)} } $$
where $R_1$ and $R_2$ are the radii of curvature of the mating teeth at the contact point. For spur gears, these radii change continuously along the path of contact. The specific sliding coefficient $\lambda_1$ for the driving gear is defined as $(v_{s1} – v_{s2}) / v_{s1}$, where $v_{s}$ are the sliding velocities. Substituting these into the wear depth equation provides a comprehensive framework for predicting wear under various loads, speeds, and material properties for spur gears.
Conclusion and Future Perspectives
This article has presented an integrated numerical methodology for accurately calculating the time-varying mesh stiffness of spur gears while accounting for realistic, non-uniform tooth surface wear. The key elements of the methodology are:
- Precise Geometric Modeling: Using the generating method principle to derive exact parametric equations for both the involute and root fillet segments of the spur gear tooth.
- Comprehensive Stiffness Calculation: Applying the potential energy method to sum contributions from Hertzian contact, bending, shear, axial compression, and fillet foundation deflections for both mating gears.
- Physics-Based Wear Prediction: Employing the Archard wear model to calculate location-dependent wear depth as a function of the meshing point’s pressure angle, leading to a realistic non-uniform wear profile.
The analysis of results leads to two primary conclusions concerning the behavior of spur gears. First, tooth surface wear follows a predictable non-uniform pattern, with maxima at the root and tip regions, a minimum near the pitch point, and a cubic relationship with the pressure angle. Second, this non-uniform wear induces a progressive reduction in the overall mesh stiffness of the gear pair. The reduction is non-uniform across the mesh cycle and is more pronounced in double-tooth-contact regions. However, for the early to moderate wear stages analyzed, the stiffness reduction is relatively small (under 2%), suggesting that initial wear effects are more significant for altering dynamic load distribution and vibration signatures than for drastically compromising static load capacity.
This modeling framework provides a crucial theoretical foundation for high-fidelity dynamic modeling of gear systems where surface degradation is a concern. It enables the simulation of how wear evolves over time and how that evolution feedbacks into changing system dynamics, which is vital for predictive maintenance algorithms. Future work can extend this model in several meaningful directions for spur gears:
- Incorporating the effect of wear debris and changing lubrication regimes on the wear coefficient $k$ and friction.
- Coupling the wear-stiffness model with a dynamic model to simulate the two-way interaction between wear progression and system vibrational response.
- Extending the methodology to include profile modifications (tip and root relief) and their interaction with wear.
- Validating the model against extensive experimental data from run-to-failure tests on spur gearboxes.
By bridging precise geometry, contact mechanics, and tribological wear, this approach offers a powerful tool for advancing the design, analysis, and health management of robust and reliable spur gear transmissions.