In the field of power transmission systems, the dynamic performance of gear pairs is a critical factor influencing system stability, noise generation, and operational lifespan. Among various gear types, the helical gear is widely favored in precision machinery, such as printing presses, due to its smoother engagement, higher load capacity, and reduced acoustic emissions compared to spur gears. The primary internal excitation within a gear transmission system originates from the time-varying meshing stiffness (TVMS), which is inherently linked to the periodic change in the number of contact teeth and the elastic deformation of the tooth structure during rotation. Accurate modeling of TVMS is therefore fundamental for any meaningful dynamic analysis, fault diagnosis, or design optimization of gear systems.
Traditional methods for calculating meshing stiffness can be broadly categorized into analytical methods, such as the Potential Energy Method (PEM), and numerical methods, primarily the Finite Element Method (FEM). While FEM offers high accuracy by simulating complex stress and deformation fields, it is computationally intensive, time-consuming, and less convenient for integration into system-level dynamic models or parametric studies. The standard PEM, which models the tooth as a non-uniform cantilever beam rooted at the base circle, provides a more efficient analytical solution. However, its classic formulation has notable limitations: it does not accurately account for the tooth fillet foundation when the base circle and root circle radii differ, and its application to helical gears requires careful consideration of the varying line of contact along the face width.
To address these shortcomings and enhance the efficiency and fidelity of dynamic analysis for gear systems, we propose an Improved Potential Energy Method (IPEM) specifically tailored for helical gears. This method incorporates a refined cantilever beam model that correctly considers the tooth root geometry relative to the root circle. It then employs Smith’s slice theory to decompose the three-dimensional helical gear contact problem into a series of two-dimensional spur gear slices, enabling the calculation of stiffness contributions along the instantaneous contact line. Furthermore, our model integrates the effects of profile modification. We validate the accuracy of the IPEM by comparing its results with those from FEM, the international standard ISO 6336-1, and the traditional PEM. Subsequently, a coupled bending-torsion-axial 9-degree-of-freedom (9-DOF) dynamic model of a triple-helical gear system, representative of a printing unit’s滚筒 drive, is established. Using the TVMS obtained from the IPEM as a key input, we perform dynamic simulations to investigate the meshing forces and vibration characteristics of the gear pairs.
The core of our analytical approach is the derivation of the meshing stiffness for a single pair of helical gear teeth. The total effective mesh stiffness \( k \) for a tooth pair is derived from the series combination of the Hertzian contact stiffness and the individual tooth stiffness components (bending, shear, axial compression, and fillet foundation stiffness) from both the driving and driven gears. The classical formula is:
$$ \frac{1}{k} = \frac{1}{k_h} + \sum_{i=1}^{2} \left( \frac{1}{k_{b,i}} + \frac{1}{k_{a,i}} + \frac{1}{k_{s,i}} + \frac{1}{k_{f,i}} \right) $$
where \( k_h \) is the Hertzian contact stiffness, and \( k_b, k_a, k_s, k_f \) represent the bending, axial compressive, shear, and fillet foundation stiffnesses, respectively, for gear \( i \) (1 for driver, 2 for driven). The traditional PEM calculates these stiffness components based on a cantilever beam model starting from the base circle. Our improvement lies in correctly defining the integration limits and tooth profile geometry relative to the actual root circle (\( r_f \)). We consider two cases:
Case 1: Base circle radius \( r_b \) \geq Root circle radius \( r_f \): The integration for energy calculation starts from the root circle.
Case 2: Base circle radius \( r_b \) < Root circle radius \( r_f \): The integration starts from the base circle, as the segment between the base and root circles is not part of the active tooth profile acting as a cantilever.
For a helical gear, the line of contact is inclined. To handle this, the gear is virtually sliced into \( N \) thin disks along its face width using Smith’s theory. Each slice is treated as a spur gear with a slightly modified pressure angle due to the helix angle \( \beta \). The load is assumed to be distributed along the instantaneous line of contact. The modified pressure angle \( \alpha'(y) \) at a distance \( y \) along the face width is:
$$ \alpha'(y) = \alpha + (\alpha_1 – \alpha) \frac{y}{l} $$
where \( \alpha \) is the standard pressure angle, \( \alpha_1 \) is the angle at the end of the contact line, and \( l \) is the total length of the contact line on the tooth surface.
Furthermore, we incorporate profile modification, typically a linear tip relief, into the calculation. The distance from the tooth centerline to a point on the profile at a radius \( r_x \) becomes \( h’_x \), which is a function of the modification amount \( C_a \) and length \( L_a \). The stiffness components are then calculated by integrating the strain energy along the tooth height and summing the contributions from all slices. The bending stiffness \( k’_b \), for instance, is calculated as follows for the two cases (illustrative formula for Case 1, \( r_b \geq r_f \)):
$$ \frac{1}{k’_b} = \sum_{j=1}^{N} \int_{0}^{d} \frac{\left[ \cos \alpha’_1 (r_x \sin \alpha’_1 + h’_x \cos \alpha’_1) – r_x \sin \alpha’_1 \right]^2}{2EI_x} \, dx $$
where \( E \) is Young’s modulus, \( I_x \) is the area moment of inertia of the tooth cross-section at distance \( x \) from the root, \( d \) is the effective cantilever length, and \( \alpha’_1 \) is evaluated for the \( j\)-th slice. Similar integral expressions are derived for the shear stiffness \( k’_s \) and axial compression stiffness \( k’_a \).
The Hertzian contact stiffness \( k’_h \) and the fillet foundation stiffness \( k’_f \) are less sensitive to the slicing and are calculated using established formulas:
$$ k’_h = \frac{\pi E B}{4(1-\nu^2)} $$
$$ \frac{1}{k’_f} = \frac{\cos^2 \alpha}{EB} \left[ L^* \left( \frac{u_f}{S_f} \right)^2 + M^* \left( \frac{u_f}{S_f} \right) + P^* (1 + Q^* \tan^2 \alpha) \right] $$
where \( B \) is face width, \( \nu \) is Poisson’s ratio, and \( L^*, M^*, P^*, Q^*, u_f, S_f \) are geometric parameters defined by the tooth and rim geometry.
The single-tooth-pair stiffness \( k_{pair}(t) \) is thus obtained by substituting \( k’_b, k’_s, k’_a, k’_f \) into the series formula. The total TVMS \( k_m(t) \) for the helical gear pair is the sum of the stiffnesses of all tooth pairs in simultaneous contact, which varies periodically with the gear rotation and the contact ratio \( \epsilon \).
The selection of the number of slices \( N \) is crucial for balancing accuracy and computational cost. Our analysis shows that for the helical gear parameters in this study, results converge stably with \( N = 1000 \), eliminating spurious fluctuations seen with lower slice counts (e.g., N=100, 200).
| Parameter | Symbol | Value (Example Gear) |
|---|---|---|
| Number of Teeth | \( z \) | 88 |
| Module | \( m_n \) | 3.5 mm |
| Pressure Angle | \( \alpha_n \) | 15° |
| Helix Angle | \( \beta \) | 18° |
| Face Width | \( B \) | 50 mm |
| Young’s Modulus | \( E \) | 211 GPa |
| Poisson’s Ratio | \( \nu \) | 0.3 |
To validate the IPEM, we compared its calculated TVMS for a helical gear pair with results from 3D Finite Element Analysis (FEA) in ABAQUS, the ISO 6336-1 standard formula (which provides a constant average value), and the traditional PEM. The FEA model consisted of over 550,000 elements to ensure accuracy. The comparison is summarized in the table below:
| Calculation Method | Max Stiffness (MN/m) | Min Stiffness (MN/m) | Mean Stiffness (MN/m) | Error vs. ISO Mean |
|---|---|---|---|---|
| ISO 6336-1 | 1036 | 1036 | 1036 | 0% |
| Improved PEM (IPEM) | 1011 | 988.3 | 998.4 | -3.67% |
| Finite Element Analysis | 1086 | 1034 | 1054 | +1.73% |
| Traditional PEM | 1148 | 1127 | 1138 | +9.84% |
The results demonstrate that the IPEM yields a mean stiffness value very close to both the FEA result (within 5.3%) and the ISO standard value. More importantly, the time-varying waveform from IPEM closely matches the trend of the FEA result, accurately capturing the periodic fluctuations due to changing contact conditions. In contrast, the traditional PEM overestimates the stiffness significantly because it incorrectly models the cantilever root. The IPEM thus offers a highly efficient and sufficiently accurate alternative to FEA for obtaining TVMS, making it ideal for dynamic simulations.
With the validated TVMS, we proceed to analyze the dynamic behavior of a triple-helical gear system simulating a printing unit drive train consisting of a Plate (P), a Blanket (B), and an Impression (I) cylinder gear. A 9-DOF lumped-parameter model is established, considering translational vibrations in the radial (y) and axial (z) directions, as well as torsional vibration (\( \theta \)) for each gear. The generalized coordinate vector is:
$$ \mathbf{q} = \{ y_B, z_B, \theta_B, y_P, z_P, \theta_P, y_I, z_I, \theta_I \}^T $$
The equations of motion are derived using Newton’s second law, incorporating the TVMS \( k_1(t), k_2(t) \) for the P-B and B-I gear pairs, respectively, meshing damping \( c_1, c_2 \), static transmission error \( e_1(t), e_2(t) \), and backlash nonlinearity \( f(\delta) \). The dynamic meshing force \( F_{ni} \) for each pair is:
$$ F_{ni} = c_i \dot{\delta}_i + k_i(t) f(\delta_i), \quad i=1,2 $$
where \( \delta_i \) is the relative displacement along the line of action. For the P-B pair:
$$ \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) $$
Here, \( \beta_b \) is the base circle helix angle. Similar equations govern the B-I pair and the translational motions. The TVMS functions \( k_1(t) \) and \( k_2(t) \) obtained from IPEM are approximated by their third-order Fourier series for efficient integration into the dynamic solver:
$$ k_i(t) = k_{mi} + \sum_{j=1}^{3} [a_{ij} \cos(j\omega_{ni} t) + b_{ij} \sin(j\omega_{ni} t)], \quad i=1,2 $$
where \( k_{mi} \) is the mean stiffness, \( \omega_{ni} \) is the meshing frequency, and \( a_{ij}, b_{ij} \) are Fourier coefficients.
The dynamic model was solved numerically in the MATLAB/Simulink environment. The analysis focused on the dynamic meshing forces and vibration responses. The key findings are as follows:
1. Dynamic Meshing Force: The time-domain responses of the meshing forces for both gear pairs are nearly identical due to symmetric design and loading. The force fluctuates between approximately 1970 N and 3445 N, with a mean value around 2700 N. This mean value has a relative error of only 4.63% compared to the theoretical static force (2832 N), validating the dynamic model’s consistency. The frequency spectrum is dominated by the meshing frequency component at 249.4 Hz.
2. Radial Vibration (y-direction): The radial vibration characteristics differ among the gears. The Blanket (B) gear exhibits the smallest vibration displacement (peak-to-peak ~2.1e-6 mm), velocity, and acceleration. This is because it is subject to opposing meshing forces from both sides. The Plate (P) and Impression (I) gears show larger and comparable radial vibrations but in opposite phases. Their displacement peak-to-peak is about 1.17e-4 mm, with acceleration peaks reaching 134 mm/s² and 89.4 mm/s², respectively.
| Vibration Metric | Plate Gear (P) | Blanket Gear (B) | Impression Gear (I) |
|---|---|---|---|
| Radial Displacement (mm) | 1.13e-4 to 2.30e-4 | -5.43e-6 to -3.32e-6 | -2.49e-4 to -1.11e-4 |
| Radial Velocity (mm/s) | -9.09e-2 to 9.09e-2 | -1.64e-3 to 1.64e-3 | -1.09e-1 to 1.09e-1 |
| Radial Acceleration (mm/s²) | -134 to 134 | -2.26 to 2.26 | -89.4 to 89.4 |
| Axial Displacement (mm) | -1.96e-6 to -1.04e-6 | 3.95e-5 to 7.59e-5 | -1.79e-4 to -1.12e-4 |
| Axial Velocity (mm/s) | -5.03e-4 to 5.03e-4 | -2.93e-2 to 2.93e-2 | -5.08e-4 to 5.08e-4 |
| Axial Acceleration (mm/s²) | -0.93 to 0.93 | -48.6 to 48.6 | -0.81 to 0.81 |
3. Axial Vibration (z-direction): The axial vibration, induced by the helix angle, presents a different pattern. The Plate (P) gear shows minimal axial movement, indicating stable axial positioning. The Blanket (B) gear, surprisingly, experiences the largest axial vibration displacement (peak-to-peak ~3.64e-5 mm) and notably high axial acceleration (peak ~48.6 mm/s²). This is a consequence of the complex coupling of radial forces from two meshing interfaces transforming into axial motion. The Impression (I) gear’s axial vibration is significant in displacement but low in acceleration, similar to the Plate gear.
In conclusion, this study presents a comprehensive framework for analyzing the meshing characteristics of helical gear pairs. The proposed Improved Potential Energy Method provides a robust, efficient, and accurate analytical tool for calculating time-varying meshing stiffness, bridging the gap between the oversimplified ISO standard, the computationally expensive FEM, and the inaccurate traditional PEM. The close agreement with FEA validates its engineering applicability. The subsequent dynamic analysis of a triple-helical gear system using this stiffness model reveals detailed insights into force transmission and vibration patterns. Key findings include the phase opposition of radial vibrations in the centrally driven gears and the potentially significant axial vibration of the intermediate idler gear, which is critical for designing thrust bearings and assessing system stability. This integrated methodology, from stiffness modeling to system dynamics, offers a powerful foundation for the design optimization, fault simulation, and performance prediction of advanced gear transmission systems, particularly in high-precision applications like printing machinery where smooth operation is paramount.