This analysis focuses on the dynamic behavior of a nutation drive system utilizing a two-stage bilateral meshing configuration with double circular arc spiral bevel gears. The study establishes a comprehensive nonlinear dynamic model to investigate the system’s vibration characteristics under operational conditions, considering key internal excitations inherent to gear meshing.
The nutation drive, inspired by planetary motion principles, offers a compact, high-ratio, and efficient transmission solution. The specific configuration examined here features two opposing meshing stages. The first-stage gear pair consists of an internal spiral bevel gear (Gear 1) and a fixed external spiral bevel gear (Gear 2). The second-stage pair involves another internal spiral bevel gear (Gear 3, rigidly connected to Gear 1 to form the nutating element) meshing with an output external spiral bevel gear (Gear 4). The input shaft axis is inclined at the nutation angle $\epsilon$ relative to the axis of the nutating gear assembly. The primary parameters for the spiral bevel gears in this system are summarized in Table 1.
| Parameter | Gear 2 (Ext.) | Gear 1 (Int.) | Gear 3 (Int.) | Gear 4 (Ext.) |
|---|---|---|---|---|
| Number of Teeth, $z$ | 26 | 28 | 32 | 30 |
| Normal Module, $m_n$ (mm) | 2 | 2 | 2 | 2 |
| Normal Pressure Angle, $\alpha_n$ (°) | 24 | 24 | 24 | 24 |
| Spiral Angle, $\beta$ (°) | 25 | 25 | 25 | 25 |
| Pitch Cone Angle, $\delta$ (°) | 47.19 | 127.81 | 123.96 | 51.04 |
The unique double circular arc profile of the spiral bevel gears contributes to favorable contact conditions but introduces significant geometric complexity, which directly influences the system’s dynamic excitation, particularly the time-varying mesh stiffness.
Dynamic Modeling of the Nutation Drive System
A coupled translational-torsional nonlinear dynamic model with 12 degrees of freedom (DOF) is developed using the lumped parameter method. The following assumptions are made: gear teeth are treated as lumped masses and inertias; the input shaft is rigid; bearing elasticity is modeled via linear springs and dampers; and friction at the gear mesh interfaces is neglected. Damping is considered as viscous.
The two internal spiral bevel gears (Gear 1 and Gear 3) are rigidly connected, sharing the same translational and rotational displacements, thus reducing their independent DOF. The dynamic model accounts for the nutation angle $\epsilon$ in the coordinate transformation for the nutating gear assembly. Two coordinate systems are defined: a fixed spatial mesh coordinate system $O-xyz$ and a rotating coordinate system $O-x_0y_0z_0$ attached to the nutating gears, where the $z_0$-axis is tilted by angle $\epsilon$ relative to the $z$-axis.
The governing equations of motion are derived using Newton’s second law, considering forces from time-varying mesh stiffness, mesh damping, static transmission error, and backlash nonlinearity. The dynamic relative displacement along the line of action for each mesh is crucial. For the first-stage mesh (between the nutating Gear 1 and fixed Gear 2), the displacement is derived as:
$$
\lambda_{n1} = \mathbf{a_1}^T \cdot \mathbf{u_{13}} – \mathbf{a_2}^T \cdot \mathbf{u_{2}} – e_{n1}(t)
$$
where $\mathbf{a_1}$ and $\mathbf{a_2}$ are projection vectors containing geometric parameters (cone angles $\delta$, pressure angle $\alpha_n$, spiral angle $\beta$, and nutation angle $\epsilon$). The displacement vectors are $\mathbf{u_{13}} = [x_{13}, y_{13}+r_{m1}\theta_{13}, z_{13}]^T$ for the nutating assembly and $\mathbf{u_{2}} = [x_{2}, y_{2}+r_{m2}\theta_{2}, z_{2}]^T$ for Gear 2. $e_{n1}(t)$ is the static transmission error. A similar expression is derived for the second-stage mesh displacement $\lambda_{n2}$.
The nonlinear mesh force due to backlash for the $j$-th mesh is defined by a piecewise function:
$$
f_{mj}(\lambda_{nj}) =
\begin{cases}
\lambda_{nj} – b, & \lambda_{nj} > b \\
0, & |\lambda_{nj}| \le b \\
\lambda_{nj} + b, & \lambda_{nj} < -b
\end{cases}
$$
where $b$ is half the total gear backlash $\Delta$, typically defined as $\Delta = 0.06m_n$ for the module range used.
The final system of nonlinear differential equations incorporates translational motions in $x$, $y$, $z$ directions and torsional motion for each independent component (Gear 2, the Nutating Assembly G1/G3, and Gear 4). The equations for Gear 2 in the $x$-direction and for the nutating assembly in the $y$-direction, for example, are:
$$
\begin{aligned}
m_2 \ddot{x}_2 + c_{2x}\dot{x}_2 + k_{2x}x_2 + a_{21}F_{n1} &= 0 \\
m_{13} \ddot{y}_{13} + 2c_{13y}\dot{y}_{13} + 2k_{13y}y_{13} – a_{22}F_{n1} + a_{42}F_{n2} &= L_{01}/r_{m1} + L_{03}/r_{m3}
\end{aligned}
$$
where $F_{nj} = k_{mj}(t) f_{mj}(\lambda_{nj}) + c_{mj}\dot{\lambda}_{nj}$ is the dynamic mesh force for mesh $j$, and $L_{0}$ represents the gyroscopic moment arising from the nutation motion. The torsional equations involve the input torque $T_2$ and output load torque $T_4$.
Analysis of Internal Excitations
The primary internal excitations in this spiral bevel gear system are the time-varying mesh stiffness $k_m(t)$ and the static transmission error $e_n(t)$.
Time-Varying Mesh Stiffness of Double Circular Arc Spiral Bevel Gears: Determining the mesh stiffness for double circular arc spiral bevel gears is complex due to the varying tooth thickness and multiple concurrent contact points along the spatially curved surface. The total mesh stiffness is the sum of the individual contact point stiffnesses: $k_m = \sum_{i=1}^{r} k_{n}^{(i)}$, where $r$ is the number of simultaneous contact points. The stiffness at a single point $k_n$ is calculated from the contact force and the associated elastic deformation (Hertzian contact and bending deflection).
Since an analytical solution is impractical, finite element analysis is employed to obtain the mesh stiffness under load. The resulting periodic stiffness function for each mesh is then approximated by an 8th-order Fourier series to facilitate numerical solution of the dynamics model:
$$
k_{mj}(t) = k_{av,j} + \sum_{n=1}^{8} \left[ a_{n,j} \sin(n\omega_{mj} t) + b_{n,j} \cos(n\omega_{mj} t) \right]
$$
where $k_{av,j}$ is the average mesh stiffness, $\omega_{mj}$ is the mesh frequency, and $a_{n,j}, b_{n,j}$ are harmonic coefficients. The average stiffness values were found to be $k_{av1}=210$ N/µm and $k_{av2}=215.4$ N/µm. The harmonic coefficients for the first mesh are listed in Table 2.
| Harmonic Order, $n$ | $a_{n1}$ (N/µm) | $b_{n1}$ (N/µm) |
|---|---|---|
| 1 | 4.071 | 16.070 |
| 2 | 17.720 | -7.198 |
| 3 | 18.790 | 25.500 |
| 4 | -7.096 | 8.304 |
| 5 | 12.960 | 7.546 |
| 6 | -5.891 | 17.990 |
| 7 | -4.250 | -0.9519 |
| 8 | 1.084 | 6.479 |
Static Transmission Error: This excitation is modeled as a sinusoidal function: $e_{nj}(t) = e_{r,j} \sin(\omega_{mj} t + \phi_j)$. For a gear accuracy grade of 7, an amplitude $e_{r,j} = 25$ µm is used, with an initial phase $\phi_j = 0$.
Dimensionless Equations and Solution Method
To generalize the analysis and improve numerical stability, the equations are non-dimensionalized. The relative displacements $\lambda_1$ and $\lambda_2$ are chosen as key coordinates to eliminate rigid-body motion. The displacement scale is half the backlash $b$, and the time scale is the inverse of the first stage’s natural frequency $\omega_n = \sqrt{k_{av1}/m_{e1}}$, where $m_{e1}$ is the equivalent mass of the first-stage mesh. The dimensionless time is $\tau = \omega_n t$.
The resulting set of coupled, nonlinear, second-order ordinary differential equations is solved numerically using a variable-step fourth-order Runge-Kutta algorithm. The dynamic response is evaluated after the transients have decayed, typically discarding the first 200 mesh cycles and analyzing the subsequent steady-state response over 300 cycles.
A key performance metric, the dynamic load factor $K_j$, is calculated for each mesh as the ratio of the dynamic mesh force $F_{nj}$ to the static load $P_{nj}$: $K_j = F_{nj} / P_{nj}$.
Dynamic Response and Parametric Influence
The nonlinear dynamics of the nutation drive system with double circular arc spiral bevel gears exhibit rich behavior influenced by operating conditions and design parameters.
Influence of Excitation Frequency: The dimensionless excitation frequency $\Omega = \omega_{m1} / \omega_n$ is varied to study the system’s response. As $\Omega$ increases, the system transitions through different dynamic states. For example, at $\Omega=0.9$, the system exhibits a 7-period subharmonic response, evidenced by a phase plot resembling a closed curve and a Poincaré map with seven distinct points. At $\Omega=1.9$, the response becomes quasi-periodic, characterized by a torus-shaped phase plot and a Poincaré map showing points arranged in closed curves. At $\Omega=2.9$, the system enters a chaotic state, indicated by a complex, non-repeating phase trajectory and a Poincaré map with a fractal-like structure of scattered points.
The variation of the maximum dynamic load factor $K_{max}$ with $\Omega$ reveals critical regions. Two significant peaks occur at $\Omega \approx 0.5$ and $\Omega \approx 1.0$, where the system experiences the most severe dynamic loading due to resonance effects. This underscores the importance of avoiding these critical speed ratios during operation or design.
Influence of Support Stiffness: The stiffness of the bearings supporting the gears significantly affects the system’s dynamic load. The analysis varies the support stiffness as a multiple of a baseline value. The results indicate that in the range of 0.5 to 3.5 times the baseline stiffness, increasing the support stiffness effectively reduces the maximum dynamic load factor $K_{max}$, thereby enhancing system stability. However, beyond a factor of about 3.5, further increases in stiffness yield diminishing returns and may even slightly increase the dynamic load due to changes in the system’s modal properties. This finding provides clear guidance for optimizing the support structure to mitigate vibrations in the nutation drive.
Model Validation
To validate the accuracy of the lumped-parameter dynamic model, its predictions were compared against results from a multi-body dynamics simulation conducted in a commercial software (ADAMS). Under identical operating conditions (input speed, load torque), the average dynamic mesh forces calculated by both methods showed close agreement, with a maximum discrepancy of less than 4.4% for either mesh stage. The lumped-parameter model predicted larger force fluctuations, which is attributed to its explicit inclusion of the time-varying mesh stiffness, a factor often simplified in penalty-based contact algorithms. This close correlation validates the fidelity of the established 12-DOF nonlinear dynamic model for the two-stage bilateral nutation drive with double circular arc spiral bevel gears.
Conclusion
This analysis successfully developed and analyzed a comprehensive nonlinear dynamic model for a two-stage bilateral nutation drive incorporating double circular arc spiral bevel gears. The model integrates key nonlinear factors such as time-varying mesh stiffness, transmission error, and backlash. Numerical simulation of the model reveals that the system’s dynamic response can transition from periodic to quasi-periodic and chaotic states as the excitation frequency increases. Critical resonance frequencies were identified at dimensionless frequencies of approximately 0.5 and 1.0, where dynamic load factors peak. Furthermore, it was demonstrated that judiciously increasing the support stiffness within a specific range (0.5 to 3.5 times baseline) can significantly reduce the dynamic load factor and improve operational stability. The results and the validated model provide a essential theoretical foundation for the dynamic design, optimization, and vibration control of high-performance nutation transmission systems utilizing advanced spiral bevel gear geometries.