In the field of mechanical engineering, straight bevel gear systems are widely employed in various applications such as automotive transmissions, aerospace mechanisms, and industrial machinery due to their ability to transmit power between intersecting shafts. However, the presence of backlash and time-varying parameters in these systems introduces nonlinear vibrations that can lead to dynamic instabilities, including bifurcations, tooth surface impacts, and non-meshing phenomena. Understanding the coupling between periodic motion, impact characteristics, and dynamic loads is crucial for designing reliable straight bevel gear systems. In this study, we investigate the parametric solution domain structures of a straight bevel gear system with backlash, focusing on the interplay between bifurcation behaviors, tooth surface impacts, non-meshing dynamics, and dynamic load coefficients. We employ an improved Continuous-Poincaré-Newton-Floquet (CPNF) method based on cell mapping principles to analyze the system’s response across a two-parameter plane defined by time-varying meshing stiffness and frequency ratios. Our findings reveal complex bifurcation patterns, including saddle-node, Hopf, period-doubling, and period-3 bifurcations, coexisting with multiple impact types. The results demonstrate that increasing the time-varying stiffness coefficient exacerbates impact and chaos phenomena, while non-meshing and dynamic load characteristics undergo abrupt changes influenced by bifurcations and impacts. This comprehensive analysis provides insights into parameter selection for optimizing straight bevel gear performance and mitigating adverse dynamic effects.
The dynamic behavior of straight bevel gear systems is governed by nonlinear factors such as time-varying meshing stiffness, backlash, and external excitations. To model this, we consider a 7-degree-of-freedom (DOF) system that accounts for translational and rotational motions along the coordinate axes. The equations of motion are derived using Newton’s second law, incorporating stiffness, damping, and gap functions. The dimensionless form of the dynamic equations facilitates numerical analysis and parametric studies. Key parameters include the time-varying meshing stiffness coefficient $\alpha$, frequency ratio $\Omega$, and backlash $b$. The meshing stiffness $k_h(t)$ is expressed as a Fourier series: $$k_h(t) = k_m + \sum_{l=1}^{N} k_{kl} \cos(l\Omega_h t + \Phi_{kl}),$$ where $k_m$ is the mean stiffness, $k_{kl}$ are harmonic amplitudes, and $\Phi_{kl}$ are phase angles. Similarly, the static transmission error $e_n(t)$ is given by: $$e_n(t) = \sum_{l=1}^{N} A_l \cos(l\Omega_h t + \Phi_l).$$ The gap function $f(\Lambda_n)$ defines the contact conditions: $$f(\Lambda_n) = \begin{cases} \Lambda_n – b_h, & \Lambda_n > b_h, \\ 0, & |\Lambda_n| \leq b_h, \\ \Lambda_n + b_h, & \Lambda_n < -b_h. \end{cases}$$ The dynamic transmission error $\Lambda_n$ along the line of action is: $$\Lambda_n = (X_1 – X_2)a_1 – (Y_1 – Y_2)a_2 – (Z_1 + r_1\theta_1 – Z_2 – r_2\theta_2)a_3 – e_n(t),$$ where $a_1 = \cos\delta_1 \sin\alpha_n$, $a_2 = \cos\delta_1 \cos\alpha_n$, $a_3 = \cos\alpha_n$, $\delta_1$ is the pinion pitch angle, and $\alpha_n$ is the normal pressure angle. The governing equations in dimensionless form are: $$\begin{aligned} &\ddot{x}_1 + 2\xi_{x1}\dot{x}_1 + 2a_4\xi_{h1}\dot{\lambda} + k_{x1}x_1 + a_4 k_{h1} f(\lambda) = 0, \\ &\ddot{y}_1 + 2\xi_{y1}\dot{y}_1 – 2a_5\xi_{h1}\dot{\lambda} + k_{y1}y_1 – a_5 k_{h1} f(\lambda) = 0, \\ &\ddot{z}_1 + 2\xi_{z1}\dot{z}_1 – 2a_3\xi_{h1}\dot{\lambda} + k_{z1}z_1 – a_3 k_{h1} f(\lambda) = 0, \\ &\ddot{x}_2 + 2\xi_{x2}\dot{x}_2 – 2a_4\xi_{h2}\dot{\lambda} + k_{x2}x_2 – a_4 k_{h2} f(\lambda) = 0, \\ &\ddot{y}_2 + 2\xi_{y2}\dot{y}_2 + 2a_5\xi_{h2}\dot{\lambda} + k_{y2}y_2 + a_5 k_{h2} f(\lambda) = 0, \\ &\ddot{z}_2 + 2\xi_{z2}\dot{z}_2 + 2a_3\xi_{h2}\dot{\lambda} + k_{z2}z_2 + a_3 k_{h2} f(\lambda) = 0, \\ &-a_1\ddot{x}_1 + a_2\ddot{y}_1 + a_3\ddot{z}_1 + a_1\ddot{x}_2 – a_2\ddot{y}_2 – a_3\ddot{z}_2 + \ddot{\lambda} + 2a_3\xi_h\dot{\lambda} + a_3 k_h f(\lambda) = f_{pm} + f_{pv} + f_e \Omega^2 \cos(\Omega\tau), \end{aligned}$$ where $\xi_{ij}$ are damping ratios, $k_{ij}$ are stiffness coefficients, and $f_{pm}$, $f_{pv}$ are dimensionless force components.
To quantify the dynamic characteristics of the straight bevel gear system, we define several performance indicators. The non-meshing duty coefficient ($\delta_{NMDC}$) and back-meshing duty coefficient ($\delta_{BMDC}$) measure the severity of tooth separation and back-side contact, respectively. These are calculated based on the minimum dynamic transmission error $\lambda_{min}$ relative to the backlash $b$: $$\delta_{NMDC} = \frac{t_1}{T}, \quad \delta_{BMDC} = \frac{t_2}{T},$$ where $t_1$ and $t_2$ are the times of non-meshing and back-meshing within one period $T$. The impact state $I$ is categorized as: $$I = \begin{cases} 0, & \lambda_{min} > b \text{ (no impact, meshing)}, \\ 0, & \lambda_{min} = b \text{ (tooth surface friction slide, meshing)}, \\ 1, & \lambda_{min} < b \text{ (unilateral impact, non-meshing)}, \\ 1, & \lambda_{min} = -b \text{ (tooth back friction slide, non-meshing)}, \\ 2, & \lambda_{min} < -b \text{ (bilateral impact, tooth back meshing)}. \end{cases}$$ The dynamic load coefficient ($\delta_{DLC}$) represents the amplification of dynamic meshing forces: $$\delta_{DLC} = \frac{|f_n(\tau)|_{max}}{f_{pm}} = \frac{|k(\tau) f(\lambda, b) + 2\xi_h \dot{\lambda}|_{max}}{f_{pm}}.$$ These indicators help assess the system’s dynamic performance and identify critical parameter regions.
We analyze the parametric solution domain structures using cell mapping theory and the CPNF method. The two-parameter plane is defined by $\Omega \times \alpha \in [1.2, 1.7] \times [0.05, 0.55]$, discretized into $501 \times 501$ cells. The CPNF method efficiently tracks periodic solutions and determines their stability via Floquet multipliers. The Jacobian matrix $\mathbf{f_X}(\mathbf{X}, \tau)$ is computed using finite differences near discontinuity points. The continuation process involves predicting new fixed points using Euler integration: $$\mathbf{X}_{k+1,0} = \mathbf{X}_k – [\mathbf{G_X}(\Omega_k, \mathbf{X}_k)]^{-1} \mathbf{G_\Omega}(\Omega_k, \mathbf{X}_k) \Delta\Omega,$$ where $\mathbf{G_X} = \mathbf{I} – \mathbf{DP}(\Omega_k, \mathbf{X}_k)$ and $\mathbf{G_\Omega} = \mathbf{P_\Omega}(\Omega_k, \mathbf{X}_k)$. This approach allows for rapid exploration of the parameter space and identification of bifurcation boundaries.
The results reveal rich dynamic behaviors in the straight bevel gear system. The $I/P$ solution domain structure (where $I$ is impact type and $P$ is period) shows predominant period-1 motion with regions of period-2, period-3, period-9, period-10, period-18, period-32, quasi-periodic, and chaotic motions. Three impact types coexist, with unilateral impact being most common. Key bifurcations include saddle-node (SN), Hopf (HP), period-doubling (PD), period-3 (3T), and catastrophe (CIC) bifurcations. For low $\alpha$ ($\alpha < 0.15$), increasing $\Omega$ induces SN and 3T bifurcations, leading to $0/1$, $1/1$, and $1/3$ states. In the range $0.15 < \alpha < 0.3$, PD, 3T, and HP bifurcations drive the system into chaos. For $0.3 < \alpha < 0.55$, grazing (G) bifurcations, 3T chaos, and CIC bifurcations occur, with complex dynamics near $\Omega = 1.6$. The non-meshing and back-meshing characteristics exhibit abrupt changes at bifurcation boundaries, with $\delta_{NMDC}$ reaching up to 0.55 in $1/1$ regions and $\delta_{BMDC}$ peaking at 0.38 in $2/1$ regions. The dynamic load coefficient $\delta_{DLC}$ shows mutations, with values up to 2.8 in bilateral impact zones. Overall, increasing $\alpha$ intensifies impact and chaos, while higher $\Omega$ reduces non-meshing and back-meshing in the same domain.
| Parameter Region | Predominant Motion | Bifurcation Types | Impact States | $\delta_{NMDC}$ Range | $\delta_{BMDC}$ Range | $\delta_{DLC}$ Range |
|---|---|---|---|---|---|---|
| $\alpha < 0.15$ | Period-1 | SN, 3T | 0/1, 1/1 | 0.1–0.3 | 0.0–0.1 | 1.2–1.8 |
| $0.15 < \alpha < 0.3$ | Period-1, Chaos | PD, HP, 3T | 1/1, 1/N | 0.2–0.5 | 0.1–0.3 | 1.5–2.5 |
| $0.3 < \alpha < 0.55$ | Period-1, Period-3 | G, 3T, CIC | 2/1, 2/N | 0.3–0.55 | 0.2–0.38 | 2.0–2.8 |
The bifurcation diagrams and Poincaré maps validate the CPNF results. For $\alpha = 0.20$, the system transitions from $1/1$ to $2/3$ via 3T bifurcation, then to quasi-periodic and chaotic motions through Hopf bifurcations. For $\alpha = 0.40$, period-3 bifurcations lead to period-9, period-18, and chaos. The phase plots show distinct attractors, confirming the period-doubling route to chaos. The straight bevel gear system’s dynamics are highly sensitive to parameter variations, emphasizing the need for careful design to avoid unstable regions.
In conclusion, the parametric solution domain analysis of straight bevel gear systems reveals intricate couplings between bifurcations, impacts, non-meshing, and dynamic loads. The CPNF method proves effective for global dynamics exploration. Designers should avoid parameter combinations that induce chaos or severe impacts, preferentially selecting regions with period-1 motion and minimal backlash effects. Future work could extend this approach to include friction, lubrication, and multi-mesh interactions in straight bevel gear systems.