In our study, we propose an improved method based on the potential energy approach to calculate the time-varying meshing stiffness of helical gears, aiming to enhance both the computational efficiency and accuracy in dynamic analysis of gear transmission systems. Helical gears are widely used in high-precision machinery such as printing presses, where their dynamic behavior directly influences print quality and operational stability. By accurately modeling the time-varying meshing stiffness, we gain deeper insights into the system’s dynamic responses, which facilitates optimal gear design and performance improvement.
1. Improved Potential Energy Method for Helical Gears
Traditional potential energy methods treat gear teeth as cantilever beams on the base circle, but this simplification leads to errors when the root circle differs from the base circle. To overcome this, we modify the calculation by considering the actual root geometry and incorporate the Smith slicing theory for helical gears. The helical gear is divided into \(N\) thin slices along the axial direction, each approximated as a spur gear with a small helix angle. The meshing force angle \(\alpha_1(y)\) for each slice is given by:
$$
\alpha_1(y) = \alpha_1 + \frac{y}{l} (\alpha – \alpha_1)
$$
where \(\alpha\) is the pressure angle, \(\alpha_1\) is the angle between the meshing force and the perpendicular to the tooth centerline, and \(l\) is the contact line length. The tooth profile modification is introduced via the corrected distance \(h’_x\) from the tooth centerline at a distance \(x\) from the root:
$$
h’_x =
\begin{cases}
h_x, & r_x \leq r_b \\
h_x + C_a \left[ \frac{r_x – r_b}{L_a} – \frac{1}{2\pi} \sin\left(2\pi\frac{r_x – r_b}{L_a}\right) \right], & r_b < r_x \leq r_a
\end{cases}
$$
where \(C_a\) is the maximum tip relief, \(L_a\) is the modification length, \(r_x\) is the radius at distance \(x\), \(r_a\) is the addendum circle radius, and \(r_b\) is the base circle radius. The bending stiffness of a single tooth slice, \(k’_b\), is integrated as:
$$
k’_b =
\begin{cases}
\displaystyle \sum_{i=1}^{N} \frac{1}{\frac{3}{EB h’_x} \int_{\alpha_1}^{\alpha_2} \frac{\left[1 + \frac{h’_x}{d}\right]^3 \cos\alpha – \cos\alpha_1\left(1 + \frac{h’_x}{d}\right)^2 \cos\alpha}{\cos\alpha}\, d\alpha}, & r_b \geq r_f \\[10pt]
\displaystyle \sum_{i=1}^{N} \frac{1}{\frac{3}{EB h’_x} \int_{\alpha_1}^{\alpha_2} \frac{\left[1 + \frac{h’_x}{d}\right]^3 \cos\alpha – \cos\alpha_1\left(1 + \frac{h’_x}{d}\right)^2 \cos\alpha}{\cos\alpha}\, d\alpha}, & r_b < r_f
\end{cases}
$$
Similarly, the shear stiffness \(k’_s\) and axial compressive stiffness \(k’_a\) are expressed as:
$$
k’_s =
\begin{cases}
\displaystyle \sum_{i=1}^{N} \frac{1}{\frac{6(1+\nu)}{5EBh’_x} \int_{\alpha_1}^{\alpha_2} \frac{\cos\alpha_1 \cos^2\alpha}{\cos\alpha}\, d\alpha}, & r_b \geq r_f \\[10pt]
\displaystyle \sum_{i=1}^{N} \frac{1}{\frac{6(1+\nu)}{5EBh’_x} \int_{\alpha_1}^{\alpha_2} \frac{\cos\alpha_1 \cos^2\alpha}{\cos\alpha}\, d\alpha}, & r_b < r_f
\end{cases}
$$
$$
k’_a =
\begin{cases}
\displaystyle \sum_{i=1}^{N} \frac{1}{\frac{1}{EB} \int_{\alpha_1}^{\alpha_2} \frac{\sin\alpha_1 \cos\alpha}{2h’_x}\, d\alpha}, & r_b \geq r_f \\[10pt]
\displaystyle \sum_{i=1}^{N} \frac{1}{\frac{1}{EB} \int_{\alpha_1}^{\alpha_2} \frac{\sin\alpha_1 \cos\alpha}{2h’_x}\, d\alpha}, & r_b < r_f
\end{cases}
$$
The Hertzian contact stiffness \(k’_h\) and gear body stiffness \(k’_f\) are given by standard formulas. The total single tooth pair stiffness \(k\) is then obtained by summing the contributions of both gears and the contact. The overall time-varying meshing stiffness of the helical gear pair is obtained by superimposing the contributions of all simultaneously engaged tooth pairs, as illustrated in the composite stiffness curve.
The number of slices \(N\) significantly affects accuracy and computational cost. We compared results for \(N = 100, 200, 500, 800, 1000, 2000\). For \(N \geq 1000\), the stiffness curve becomes smooth and stable, so we chose \(N = 1000\) for all subsequent calculations.
2. Validation of the Improved Method
To validate our improved potential energy method, we compare its results with those obtained from the finite element method (FEM), the ISO 6336-1 standard, and the traditional potential energy method. The FEM model was built in SolidWorks and meshed in HyperMesh with 6,111,016 nodes and 553,960 elements, then solved in ABAQUS. The contact forces and elastic deformations were extracted to compute meshing stiffness.
The comparison results are summarized in the table below:
| Method | Max (MN·m⁻¹) | Min (MN·m⁻¹) | Mean (MN·m⁻¹) | Error vs ISO (%) |
|---|---|---|---|---|
| ISO 6336-1 | 1036 | 1036 | 1036 | 0 |
| Improved Potential Energy | 1011 | 988.3 | 998.4 | –3.67 |
| Finite Element Method | 1086 | 1034 | 1054 | +1.73 |
| Traditional Potential Energy | 1148 | 1127 | 1138 | +9.84 |
The improved method shows excellent agreement with the ISO standard (within 5% error) and closely follows the FEM results. Our method significantly reduces computational time compared to FEM while correcting the overestimation of traditional potential energy. The improved potential energy method thus provides a fast, accurate, and practical tool for calculating the time-varying meshing stiffness of helical gears.
3. Dynamic Modeling of Helical Gear Transmission System
We developed a nine-degree-of-freedom bending-torsion-axial coupled dynamic model for the helical gear transmission system of a printing press, comprising the plate cylinder gear (P), blanket cylinder gear (B), and impression cylinder gear (I). The generalized coordinate vector is:
$$
\{q\}^T = \{y_P, z_P, \theta_P, y_B, z_B, \theta_B, y_I, z_I, \theta_I\}
$$
Considering time-varying meshing stiffness, backlash, transmission error, relative displacement, and meshing damping, the differential equations of motion are:
$$
\begin{cases}
m_P \ddot{y}_P + c_{Py} \dot{y}_P + k_{Py} y_P = F_{y1} \\
m_P \ddot{z}_P + c_{Pz} \dot{z}_P + k_{Pz} z_P = -F_{z1} \\
I_P \ddot{\theta}_P = T_P – r_{bP} F_{n1} \\
m_B \ddot{y}_B + c_{By} \dot{y}_B + k_{By} y_B = -F_{y1} – F_{y2} \\
m_B \ddot{z}_B + c_{Bz} \dot{z}_B + k_{Bz} z_B = F_{z1} + F_{z2} \\
I_B \ddot{\theta}_B = -T_B + r_{bB} F_{n1} – r_{bB} F_{n2} \\
m_I \ddot{y}_I + c_{Iy} \dot{y}_I + k_{Iy} y_I = F_{y2} \\
m_I \ddot{z}_I + c_{Iz} \dot{z}_I + k_{Iz} z_I = -F_{z2} \\
I_I \ddot{\theta}_I = -T_I + r_{bI} F_{n2}
\end{cases}
$$
where \(F_{yi} = F_{ni} \cos\beta_b\), \(F_{zi} = F_{ni} \sin\beta_b\), and the normal meshing force for each gear pair is:
$$
F_{ni} = c_i \dot{\delta}_i + k_i(t) f(\delta_i), \quad i = 1,2
$$
The relative displacements \(\delta_1\) and \(\delta_2\) along the line of action for the P-B and B-I pairs are:
$$
\delta_1 = (y_P – y_B)\cos\beta_b + (z_P – z_B)\sin\beta_b + r_{bP}\theta_P – r_{bB}\theta_B – e_1(t)
$$
$$
\delta_2 = (y_B – y_I)\cos\beta_b + (z_B – z_I)\sin\beta_b + r_{bB}\theta_B – r_{bI}\theta_I – e_2(t)
$$
The time-varying mesh stiffness functions \(k_1(t)\) and \(k_2(t)\) obtained from the improved potential energy method are fitted with a third-order Fourier series:
$$
k_1(t) = k_{m1} + \sum_{j=1}^{3} \left[ a_{1j} \cos(j\omega_1 t) + b_{1j} \sin(j\omega_1 t) \right]
$$
$$
k_2(t) = k_{m2} + \sum_{j=1}^{3} \left[ a_{2j} \cos(j\omega_2 t) + b_{2j} \sin(j\omega_2 t) \right]
$$
4. Dynamic Analysis and Results
We solved the dynamic model in MATLAB using the parameters listed in the table below (identical for all three gears).
| Parameter | Plate Cylinder Gear P | Blanket Cylinder Gear B |
|---|---|---|
| Number of teeth | 88 | 88 |
| Module (mm) | 3.5 | 3.5 |
| Pressure angle (°) | 15 | 15 |
| Face width (mm) | 50 | 50 |
| Helix angle (°) | 18 | 18 |
| Speed (rpm) | 170 | 170 |
| Material | 40Cr | 40Cr |
| Elastic modulus (GPa) | 211 | 211 |
| Poisson’s ratio | 0.277 | 0.277 |
| Mass (kg) | 25 | 25 |
| Moment of inertia (kg·m²) | 0.40 | 0.40 |
| Backlash (μm) | 0.1 | 0.1 |
4.1 Dynamic Meshing Force
The dynamic meshing forces of the P-B and B-I gear pairs are shown in time and frequency domains. The P-B pair fluctuates between 1971.2 N and 3445.4 N with a peak-to-peak value of 1474.2 N and a mean of 2700.9 N. The B-I pair fluctuates between 1972.8 N and 3443.7 N with a peak-to-peak of 1470.9 N and a mean of 2701.1 N. Both pairs exhibit nearly identical trends due to similar gear parameters, with slight differences arising from phase shift. The theoretical mean force is 2832.2 N, giving a relative error of about 4.63%, validating our dynamic model. Frequency analysis shows dominant peaks at 249.4 Hz and 22.6 Hz, with the B-I pair having slightly lower amplitudes.
4.2 Radial Vibration Analysis
Radial vibrations (y-direction) of the three gears are summarized below:
| Parameter | Gear P | Gear B | Gear I |
|---|---|---|---|
| Displacement range (mm) | 1.13×10⁻⁴ ~ 2.30×10⁻⁴ | –5.43×10⁻⁶ ~ –3.32×10⁻⁶ | –2.49×10⁻⁴ ~ –1.11×10⁻⁴ |
| Velocity range (mm/s) | –9.09×10⁻² ~ 9.09×10⁻² | –1.64×10⁻³ ~ 1.64×10⁻³ | –1.09×10⁻¹ ~ 1.09×10⁻¹ |
| Acceleration range (mm/s²) | –134 ~ 134 | –2.26 ~ 2.26 | –89.4 ~ 89.4 |
Gear B (blanket cylinder) exhibits the smallest radial vibration, while gears P and I show similar amplitudes but opposite directions. This indicates that the blanket cylinder gear experiences balanced forces from both sides, leading to reduced radial motion.
4.3 Axial Vibration Analysis
Axial vibrations (z-direction) are presented in the table below:
| Parameter | Gear P | Gear B | Gear I |
|---|---|---|---|
| Displacement range (mm) | –1.96×10⁻⁶ ~ –1.04×10⁻⁶ | 3.95×10⁻⁵ ~ 7.59×10⁻⁵ | –1.79×10⁻⁴ ~ –1.12×10⁻⁴ |
| Velocity range (mm/s) | –5.03×10⁻⁴ ~ 5.03×10⁻⁴ | –2.93×10⁻² ~ 2.93×10⁻² | –5.08×10⁻⁴ ~ 5.08×10⁻⁴ |
| Acceleration range (mm/s²) | –0.93 ~ 0.93 | –48.6 ~ 48.6 | –0.81 ~ 0.81 |
Gear B shows the largest axial vibrations, while gears P and I have similar and much smaller axial motions. This is because gear B is subjected to axial forces from both adjacent gears in opposite directions, increasing its axial excitation. The axial displacement of gear P is minimal, indicating stable performance.
5. Conclusion
We proposed an improved potential energy method for computing the time-varying meshing stiffness of helical gears, which accounts for the root circle effect and tooth profile modifications. The method yields results within 5% of the ISO standard and close to FEM, while significantly reducing computational burden. By integrating this stiffness into a nine-DOF dynamic model of a printing press helical gear transmission, we analyzed dynamic meshing forces and vibrations. The blanket cylinder gear shows smaller radial vibration but larger axial vibration compared to the plate and impression cylinder gears. These findings provide valuable insights for optimizing the design and performance of helical gear systems in high-precision applications.