An Analytical Method for Time-Varying Mesh Stiffness Calculation of Spur Gears Considering Friction Effects

Dynamic performance prediction and monitoring through vibration analysis is an effective approach for optimizing the design and diagnosing faults in gear transmission systems. Among the various excitations in gear dynamics, the time-varying mesh stiffness (TVMS), caused by the alternating number of tooth pairs in contact, is a primary source of parametric excitation. The frequency and amplitude of TVMS variations significantly influence the dynamic response of the system. Therefore, accurately calculating the TVMS is fundamental to reliable gear dynamic analysis.

Current mainstream methods for calculating TVMS are the potential energy method and the finite element method (FEM). The potential energy method, pioneered by Yang and Lin, calculates stiffness by considering the elastic energy stored in the meshing teeth, including Hertzian contact energy, bending energy, and axial compressive energy. Subsequent research enhanced this model by incorporating shear energy and later, the gear body flexibility. While the potential energy method is now mature and efficient, and FEM can handle complex contact scenarios, both approaches typically treat the problem statically and neglect the influence of a crucial dynamic parameter: friction.

In the meshing process of spur gears, friction is a complex, time-varying parameter influenced by tooth geometry, surface roughness, lubrication regime, load, temperature, curvature radius, and sliding velocity. This friction induces additional elastic deformation not accounted for in traditional models. The friction coefficient, and consequently the friction force, changes continuously during the mesh cycle, making the actual dynamic deformation more complex than what is predicted by frictionless models. Ignoring this effect may lead to inaccuracies in the calculated TVMS and, subsequently, in dynamic simulations.

This paper addresses this gap. An improved analytical model for calculating the TVMS of spur gears is proposed, which incorporates the influence of time-varying tooth friction. The foundation is a friction coefficient model derived from extensive ball-on-disk experiments under elastohydrodynamic lubrication (EHL) conditions. To apply this model accurately to spur gears, key time-varying parameters such as surface velocity, slide-to-roll ratio, entraining velocity, and Hertzian pressure are derived based on the precise geometry of involute spur gears. Crucially, the model accounts for the load-sharing variations between single and double tooth contact zones. This enhanced friction model is then integrated into the classical potential energy method framework. The proposed method calculates the TVMS under different operating loads, and its results are validated against finite element analysis. This work provides a new, more comprehensive theoretical basis for TVMS, paving the way for more accurate nonlinear dynamic analysis of gear systems that include friction effects.

Precise Tooth Profile Formulation

An accurate mathematical representation of the tooth profile, including both the involute flank and the fillet (transition curve), is essential for precise calculation of contact geometry and subsequent TVMS. The profile is generated based on the rack-cutter method, as illustrated in the principle of generating method. The coordinates of any point on the gear tooth profile after generation can be expressed as:

$$
\begin{cases}
x_i = (x_{i0} – r_b \varphi_i) \cos \varphi_i + (y_{i0} + r_b) \sin \varphi_i \\
y_i = -(x_{i0} – r_b \varphi_i) \sin \varphi_i + (y_{i0} + r_b) \cos \varphi_i
\end{cases}
\quad (i = 1, 2)
$$

Where the subscript $i=1$ denotes points on the involute segment and $i=2$ denotes points on the fillet (transition curve) segment. $ (x_{i0}, y_{i0}) $ are the coordinates of a point on the rack-cutter, $ r_b $ is the base circle radius, and $ \varphi_i $ is the roll angle during the generation process.

For the Involute Segment (i=1): The rack-cutter coordinates are defined by the pressure angle $\alpha$: $x_{10} = r_b \tan \alpha$, $y_{10} = r_b / \tan \alpha$. The roll angle $\varphi_1$ is related to the pressure angle $\alpha_{\varphi}$ at the specific point on the involute by: $\varphi_1 = \tan \alpha_{\varphi} – \tan \alpha$.
For the Fillet Segment (i=2): The rack-cutter coordinates correspond to the circular tip of the tool: $ (x_{20}, y_{20}) = (x_k + r_0 \sin \beta, \; y_k – r_0 \cos \beta) $, where $ (x_k, y_k) $ is the intersection point between the straight line and the circular tip of the rack-cutter, $ r_0 $ is the tip radius, and the relationship between $\beta$ and $\varphi_2$ is given by: $\beta = \arctan\left(\frac{x_k – r_0 \sin \alpha}{r_b \varphi_2 + r_b \tan \alpha – y_k + r_0 \cos \alpha}\right)$.

The generated involute profile is then rotated to its correct angular position for meshing analysis. The resulting precise profile of an involute spur gear and its associated force analysis is crucial for the subsequent steps.

Potential Energy Method Incorporating Friction

The total elastic energy stored in a pair of meshing spur gears under a normal load $ F_n $ consists of several components: Hertzian contact energy $ U_h $, bending energy $ U_b $, shear energy $ U_s $, axial compressive energy $ U_a $, and gear body foundation deflection energy $ U_f $. According to the potential energy principle, the relationship between each energy component and its corresponding stiffness $ K $ is:

$$
U_h = \frac{F_n^2}{2K_h}, \quad U_b = \frac{F_n^2}{2K_b}, \quad U_s = \frac{F_n^2}{2K_s}, \quad U_a = \frac{F_n^2}{2K_a}, \quad U_f = \frac{F_n^2}{2K_f}
$$

The total energy $ U $ for a gear pair is the sum of energies from both the driving and driven gears:

$$
U = \sum_{i=1}^{2} (U_{hi} + U_{bi} + U_{si} + U_{ai} + U_{fi}) = \frac{F_n^2}{2K_{total}}
$$

Therefore, the total mesh stiffness $ K_{total} $ for one meshing tooth pair is:

$$
\frac{1}{K_{total}} = \sum_{i=1}^{2} \left( \frac{1}{K_{bi}} + \frac{1}{K_{si}} + \frac{1}{K_{ai}} + \frac{1}{K_{fi}} \right) + \frac{1}{K_h}
$$

The Hertzian contact stiffness $ K_h $ and the foundation stiffness $ K_f $ can be obtained from established formulas. The Hertzian stiffness for two cylinders in contact (approximating the teeth) is:

$$
K_h = \frac{\pi E L}{4(1-\nu^2)}
$$

where $ E $ is Young’s modulus, $ \nu $ is Poisson’s ratio, and $ L $ is the face width.

The foundation stiffness accounting for the deflection of the gear body is given by a more complex formula involving the geometry near the tooth root. Based on Sainsot’s work, it can be expressed as:

$$
\frac{1}{K_f} = \frac{\cos^2 \alpha_{\varphi}}{EL} \left[ L^* \left( \frac{u_f}{S_f} \right)^2 + M^* \left( \frac{u_f}{S_f} \right) + P^* (1+Q^* \tan^2 \alpha_{\varphi}) \right]
$$

where $ \alpha_{\varphi} $ is the pressure angle at the load application point, $ u_f $ and $ S_f $ are geometric parameters defined from the tooth root, and $ L^*, M^*, P^*, Q^* $ are dimensionless constants dependent on the gear geometry.

The bending, shear, and axial compression energies are calculated by modeling the tooth as a non-uniform cantilever beam. When friction is considered, the normal force $ F_n $ resolves into components that include the friction force $ \mu F_n $. The direction of the friction force reverses at the pitch point. Referring to the gear meshing process and changing friction force diagram, the force components acting on the tooth at a distance $ x $ from the root are:

$$
\begin{aligned}
F_y &= \begin{cases}
F_n \cos \alpha_{\varphi} – \mu F_n \sin \alpha_{\varphi}, & \text{for } \alpha_{\varphi} < \alpha_{\text{pitch}} \\
F_n \cos \alpha_{\varphi} + \mu F_n \sin \alpha_{\varphi}, & \text{for } \alpha_{\varphi} \ge \alpha_{\text{pitch}}
\end{cases} \\
F_x &= \begin{cases}
F_n \sin \alpha_{\varphi} + \mu F_n \cos \alpha_{\varphi}, & \text{for } \alpha_{\varphi} < \alpha_{\text{pitch}} \\
F_n \sin \alpha_{\varphi} – \mu F_n \cos \alpha_{\varphi}, & \text{for } \alpha_{\varphi} \ge \alpha_{\text{pitch}}
\end{cases}
\end{aligned}
$$

The bending moment $ M $ at section $ x $ is: $ M = F_y \cdot (x_{\varphi} – x) – F_x \cdot y_{\varphi} $, where $ (x_{\varphi}, y_{\varphi}) $ are the coordinates of the load application point.

Using beam theory, the stored energies are:

$$
U_b = \int_{0}^{x_{\varphi}} \frac{M^2}{2EI_x} dx, \quad U_s = \int_{0}^{x_{\varphi}} \frac{1.2 F_y^2}{2GA_x} dx, \quad U_a = \int_{0}^{x_{\varphi}} \frac{F_x^2}{2EA_x} dx
$$

where $ I_x $ is the area moment of inertia, $ A_x $ is the cross-sectional area (both varying with $ x $), and $ G = E/[2(1+\nu)] $ is the shear modulus. The factor 1.2 accounts for the non-uniform shear stress distribution. The integrations are performed from the tooth root ($x=0$) to the load point ($x=x_{\varphi}$). The individual tooth stiffness components ($K_b, K_s, K_a$) are then derived by substituting these energy expressions back into $ K = F_n^2/(2U) $.

Time-Varying Friction Coefficient under EHL Conditions

Friction in lubricated spur gears is typically governed by elastohydrodynamic lubrication (EHL) regimes. An empirical formula for the EHL friction coefficient $ \mu $, derived from regression analysis of ball-on-disk test data, is used as the basis:

$$
\mu(SR, P, v_0, S, u_e, R) = e^{f(SR, P, v_0, S)} \cdot (u_e)^{b_8} \cdot (R)^{b_9}
$$

where the function $ f $ is given by:

$$
f(SR, P, v_0, S) = b_1 + b_2 \cdot (\log_{10}(v_0))^2 + b_3 \cdot S^2 + b_4 \cdot e^{SR} + b_5 \cdot e^{P} + (b_6 \cdot (\log_{10}(v_0))^2 + b_7 \cdot S^2) \cdot P
$$

The parameters are:

$ SR $: Slide-to-roll ratio.
$ P $: Maximum Hertzian contact pressure.
$ v_0 $: Kinematic viscosity of the lubricant at atmospheric pressure.
$ S $: Root mean square (RMS) surface roughness.
$ u_e $: Entraining (rolling) velocity.
$ R $: Equivalent radius of curvature at the contact point.

The regression coefficients are: $b_1 = -8.916465$, $b_2 = 1.03303$, $b_3 = 1.036077$, $b_4 = -0.354068$, $b_5 = 2.812084$, $b_6 = -0.100601$, $b_7 = 0.752755$, $b_8 = -0.390958$, $b_9 = 0.620305$.

To apply this formula to meshing spur gears, the time-varying parameters $ SR, P, u_e, R $ must be expressed as functions of the gear geometry and operating conditions. The calculations are performed along the line of action (LOA).

Velocities: The tangential velocities of the pinion ($p$) and gear ($g$) at the pitch point are $v_p = r_{bp} \omega_p$ and $v_g = r_{bg} \omega_g$, where $r_b$ is the base radius and $\omega$ is the angular velocity. Along the LOA, the velocities in the direction normal to the teeth are equal: $v_{LOA} = v_p = v_g$. The tangential surface velocities at a generic contact point with pressure angle $\alpha_{\varphi}$ are:
$$
v_{t,p} = \frac{v_{LOA}}{\cos \alpha_{\varphi,p}}, \quad v_{t,g} = \frac{v_{LOA}}{\cos \alpha_{\varphi,g}}
$$
These can be related to the pitch point velocities: $v_{t,p} = v_p \cdot \tan \alpha_{\varphi,p}$, $v_{t,g} = v_g \cdot \tan \alpha_{\varphi,g}$.

Sliding and Entraining Velocities:
$$
v_s = |v_{t,p} – v_{t,g}| = v_{LOA} \cdot |\tan \alpha_{\varphi,p} – \tan \alpha_{\varphi,g}|
$$
$$
u_e = \frac{v_{t,p} + v_{t,g}}{2} = \frac{v_{LOA}}{2} \cdot (\tan \alpha_{\varphi,p} + \tan \alpha_{\varphi,g})
$$

Slide-to-Roll Ratio:
$$
SR = \frac{v_s}{u_e} = 2 \cdot \frac{ |\tan \alpha_{\varphi,p} – \tan \alpha_{\varphi,g}| }{ \tan \alpha_{\varphi,p} + \tan \alpha_{\varphi,g} }
$$

Equivalent Radius of Curvature: The radii of curvature for pinion and gear are $ \rho_p = r_{bp} \tan \alpha_{\varphi,p} $ and $ \rho_g = r_{bg} \tan \alpha_{\varphi,g} $. The equivalent radius is:
$$
R = \frac{\rho_p \cdot \rho_g}{\rho_p + \rho_g}
$$

Hertzian Contact Pressure: The maximum pressure for a line contact is:
$$
P = \sqrt{ \frac{ F_n / L }{ \pi R } \cdot \frac{1}{ \frac{1-\nu_p^2}{E_p} + \frac{1-\nu_g^2}{E_g} } }
$$
Crucially, the normal load per unit length $ F_n/L $ is not constant. In the double tooth contact zone, the total torque is shared between two pairs, while in the single tooth contact zone, one pair carries the full load. This load variation must be factored into the calculation of $ P $ at each meshing position. A comparison of friction coefficients under different loads shows a significant increase with higher torque/pressure.

Table 1: Geometric Parameters and Material Properties of the Example Spur Gears
Parameter	Value
Young’s Modulus, $E$ (GPa)	206
Poisson’s Ratio, $\nu$	0.3
Pressure Angle, $\alpha$ (deg)	20
Face Width, $L$ (mm)	20
Number of Teeth, $z_p / z_g$	45 / 45
Module, $m$ (mm)	3
Surface Roughness, $S$ ($\mu m$, RMS)	1.6
Addendum Coefficient	1
Dedendum/Clearance Coefficient	0.25
Contact Ratio, $\epsilon$	1.7358

Results and Analysis

Using the described methodology, the TVMS for the spur gear pair defined in Table 1 was calculated under different pinion input torques $T_p$ (100 N·m, 500 N·m, 1000 N·m) and a constant speed. The meshing position is normalized from $ \epsilon – 1 $ (the previous tooth pair’s exit) to 1 (the next tooth pair’s entry), where $ \epsilon $ is the contact ratio.

Single Tooth Pair Stiffness

The time-varying stiffness for a single meshing tooth pair is shown in the corresponding figure. Key observations are:

Load Dependence: As the transmitted torque increases, the friction coefficient rises, leading to a reduction in the calculated mesh stiffness. Compared to the frictionless case, the root mean square (RMS) value of the single-pair stiffness decreased by 1.86%, 2.24%, and 2.44% for torques of 100 N·m, 500 N·m, and 1000 N·m, respectively.
Pitch Point Behavior: At the pitch point ($\alpha_{\varphi} = \alpha$), the theoretical sliding velocity is zero ($SR=0$), resulting in a friction coefficient of zero. Therefore, the stiffness value at this point is identical to the frictionless calculation.
Transition Points: At the boundaries between single and double tooth contact zones (critical positions), the stiffness exhibits a sudden jump. This is due to the abrupt change in the load carried by the individual tooth pair as another pair enters or leaves contact, which affects both the Hertzian pressure and the associated friction force.

Gear Pair Mesh Stiffness

The total TVMS for the gear pair is obtained by summing the stiffness of all tooth pairs in contact simultaneously. The results for the total TVMS under different torques are presented. The influence of friction is amplified in this combined stiffness. The RMS value of the total mesh stiffness decreased by 2.22%, 2.65%, and 2.86% for the respective torque levels. This confirms that friction has a substantial impact on the primary excitation parameter in gear dynamics.

Finite Element Validation

To verify the proposed analytical method, a finite element model was constructed for the gear pair under $T_p = 1000$ N·m. The model included multiple teeth to simulate the contact transition. Contact regions were finely meshed, and friction was applied according to the time-varying model. The mesh stiffness was determined by relating the applied torque to the relative angular displacement of the gears. A comparison of stiffness values at several discrete mesh positions shows excellent agreement, as summarized in Table 2.

Table 2: Stiffness Comparison: Potential Energy Method vs. Finite Element Method (Tp=1000 N·m)
Meshing Position (Normalized)	Potential Energy Method (N/m)	Finite Element Method (N/m)	Relative Error (%)
0	6.164 × 10⁸	6.10 × 10⁸	0.97
0.2	6.144 × 10⁸	6.07 × 10⁸	1.2
0.4	3.695 × 10⁸	3.67 × 10⁸	0.75
0.6	3.678 × 10⁸	3.66 × 10⁸	0.48

The minor relative errors validate the reliability and accuracy of the proposed analytical method for calculating the TVMS of spur gears with friction.

Conclusions

This paper addresses the limitation of conventional time-varying mesh stiffness models that neglect tooth friction. An improved analytical method based on the potential energy principle is developed, incorporating a time-varying friction coefficient derived for EHL conditions and adapted for spur gears. The key contributions and findings are:

Integrated Friction Model: The proposed method successfully integrates a realistic, load- and kinematics-dependent friction model into the TVMS calculation framework for spur gears. It accounts for the reversal of friction force direction at the pitch point and the load-sharing variation in single/double contact zones.
Significant Friction Influence: Friction has a measurable and non-negligible effect on TVMS. As the transmitted load increases, the friction coefficient rises, leading to a decrease in the overall mesh stiffness. The reduction in the RMS value of the total TVMS can approach 3% under high torque for the studied case.
Characteristic Stiffness Behavior: The calculated TVMS exhibits the expected characteristics: it matches the frictionless value at the pitch point (zero sliding), shows discontinuities at the transition points between contact zones, and varies smoothly elsewhere.
Method Validity and Generality: The results are in excellent agreement with finite element analysis, validating the method. While the specific EHL friction formula was used here, the proposed analytical framework is general and can accommodate other friction models representing different lubrication states or surface conditions.

By providing a more physically accurate model for the TVMS of spur gears, this work establishes a stronger foundation for subsequent nonlinear dynamic analysis of gear systems. Considering friction in the stiffness excitation can lead to more precise predictions of vibration response, which is crucial for advanced gear design, performance optimization, and condition monitoring. The method is readily applicable to the dynamic analysis of other gear types, such as helical or bevel gears, with appropriate modifications to the geometric and kinematic relationships.