In mechanical transmission systems, straight bevel gears are widely used due to their ability to transmit motion and power between intersecting shafts. However, the presence of nonlinear factors such as backlash and time-varying mesh stiffness can lead to complex dynamic behaviors, including amplitude jumps, multi-valued solutions, and tooth surface impacts. These phenomena significantly affect the reliability and performance of straight bevel gear systems. In this study, we investigate the dynamic characteristics of a straight bevel gear system with backlash using the harmonic balance method. We establish the nonlinear dynamic equations and employ numerical techniques like the Broyden quasi-Newton method and pseudo-arc-length continuation algorithm to analyze the solution domain structure in parameter planes. Our focus is on exploring the sensitivity of amplitude jumps, multi-valued solutions, and tooth surface impacts to key parameters, providing insights for the design and optimization of straight bevel gear systems.
The dynamic model of the straight bevel gear system considers a 7-degree-of-freedom system with elastic supports, as illustrated in the figure below. The system includes vibrations along three coordinate directions and rotations about the axes. The relative torsional displacement along the mesh line is defined, accounting for vibrations and static transmission errors. The mesh force and its components are derived, incorporating time-varying mesh stiffness and a nonlinear backlash function. The equations of motion are formulated and normalized for analysis.
The nonlinear dynamic equations for the straight bevel gear system are given by:
$$ \begin{aligned}
m_1 \ddot{X}_1 + c_{x1} \dot{X}_1 + k_{x1} X_1 &= F_x \\
m_1 \ddot{Y}_1 + c_{y1} \dot{Y}_1 + k_{y1} Y_1 &= F_y \\
m_1 \ddot{Z}_1 + c_{z1} \dot{Z}_1 + k_{z1} Z_1 &= F_z \\
J_1 \ddot{\theta}_1 &= T_1 – F_z r_1 \\
m_2 \ddot{X}_2 + c_{x2} \dot{X}_2 + k_{x2} X_2 &= -F_x \\
m_2 \ddot{Y}_2 + c_{y2} \dot{Y}_2 + k_{y2} Y_2 &= -F_y \\
m_2 \ddot{Z}_2 + c_{z2} \dot{Z}_2 + k_{z2} Z_2 &= -F_z \\
J_2 \ddot{\theta}_2 &= -T_2 + F_z r_2
\end{aligned} $$
where $m_1$ and $m_2$ are the masses, $J_1$ and $J_2$ are the moments of inertia, $c_{ij}$ and $k_{ij}$ are damping and stiffness coefficients, $F_x$, $F_y$, $F_z$ are the mesh force components, and $T_1$ and $T_2$ are the torque inputs. The relative displacement $\Lambda_n$ along the mesh line is defined as:
$$ \Lambda_n = (X_1 – X_2) a_1 – (Y_1 – Y_2) a_2 – (Z_1 – Z_2) + r_1 \theta_1 – r_2 \theta_2 – e_n(t) $$
with $a_1 = \cos \delta_1 \sin \alpha_n$, $a_2 = \cos \delta_1 \cos \alpha_n$, $a_3 = \cos \alpha_n$, where $\delta_1$ is the pitch cone angle and $\alpha_n$ is the normal pressure angle. The static transmission error $e_n(t)$ is expressed as a Fourier series:
$$ e_n(t) = \sum_{l=1}^N A_l \cos(l \Omega_h t + \Phi_l) $$
The mesh force $F_n$ includes the time-varying mesh stiffness $k_h(t)$ and damping $c_h$:
$$ F_n = k_h(t) f(\Lambda_n) + c_h \dot{\Lambda}_n $$
where $f(\Lambda_n)$ is the nonlinear backlash function:
$$ 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 time-varying mesh stiffness $k_h(t)$ is also represented as a Fourier series:
$$ k_h(t) = k_m + \sum_{l=1}^N A_{kl} \cos(l \Omega_h t + \Phi_{kl}) $$
After normalization, the dimensionless equations are derived. The dimensionless relative displacement $\lambda = \Lambda_n / b_h$ is introduced, and the system is analyzed in the frequency domain using the harmonic balance method. The steady-state response is approximated by a first-order harmonic solution:
$$ x_l = x_{ml} + x_{cl} \cos(\Omega \tau) + x_{sl} \sin(\Omega \tau) $$
where $x_{ml}$, $x_{cl}$, and $x_{sl}$ are the mean, cosine, and sine components, respectively. The backlash function is linearized using describing functions:
$$ f(x_l, b_l) = N_{ml} x_{ml} + N_{al} x_{cl} \cos(\Omega \tau) + N_{al} x_{sl} \sin(\Omega \tau) $$
with the describing functions $N_{ml}$ and $N_{al}$ given by:
$$ \begin{aligned}
N_{ml} &= 1 + x_{al} [G(\mu_{l+}) – G(\mu_{l-})] / (2 x_{ml}) \\
N_{al} &= 1 – [H(\mu_{l+}) – H(\mu_{l-})] / 2
\end{aligned} $$
where $\mu_{l\pm} = (\pm b_l – x_{ml}) / x_{al}$, and the functions $G(\mu)$ and $H(\mu)$ are defined as:
$$ G(\mu) = \begin{cases}
\frac{2}{\pi} (\mu \arcsin \mu + \sqrt{1 – \mu^2}) & |\mu| \leq 1 \\
|\mu| & |\mu| > 1
\end{cases} $$
$$ H(\mu) = \begin{cases}
\frac{2}{\pi} (\arcsin \mu + \mu \sqrt{1 – \mu^2}) & |\mu| \leq 1 \\
\text{sign}(\mu) & |\mu| > 1
\end{cases} $$
Applying the harmonic balance method, we obtain a system of nonlinear algebraic equations:
$$ \begin{aligned}
\mathbf{K}_m \mathbf{G}_m – \mathbf{F}_m &= 0 \\
-\Omega^2 \mathbf{M} \mathbf{X}_c + \Omega \mathbf{C} \mathbf{X}_s + \mathbf{K}_m \mathbf{G}_c + \mathbf{K}_c \mathbf{G}_m – \mathbf{F}_c &= 0 \\
-\Omega^2 \mathbf{M} \mathbf{X}_s + \Omega \mathbf{C} \mathbf{X}_c + \mathbf{K}_m \mathbf{G}_s + \mathbf{K}_s \mathbf{G}_m – \mathbf{F}_s &= 0
\end{aligned} $$
where $\mathbf{G}_m$, $\mathbf{G}_c$, $\mathbf{G}_s$ are vectors of the describing functions, $\mathbf{K}_m$, $\mathbf{K}_c$, $\mathbf{K}_s$ are stiffness matrices, $\mathbf{M}$ is the mass matrix, $\mathbf{C}$ is the damping matrix, and $\mathbf{F}_m$, $\mathbf{F}_c$, $\mathbf{F}_s$ are force vectors. This system is solved using the Broyden quasi-Newton method and pseudo-arc-length continuation algorithm to capture the global solution curves and multi-valued responses.
For numerical simulation, we adopt dimensionless parameters: $\xi_{i1} = \xi_{i2} = 0.01$, $\xi_{h1} = \xi_{h2} = 0.0125$, $\xi_h = 0.05$, $k_{i1} = k_{i2} = 1.0$ (for $i = x, y, z$), $k_{h1} = k_{h2} = 0.5$, $k_h = 1 + \alpha \cos(\Omega \tau)$, $f_{pm} = 0.5$, $f_{pv} = 0.0$, $f_e = 0.2 \Omega^2 \cos(\Omega \tau)$, $\alpha = 0.2$, and $b = 1.0$. The frequency range is $\Omega \in [0.05, 2.00]$. We analyze the effects of backlash $b$, time-varying mesh stiffness coefficient $\alpha$, static transmission error coefficient $f_e$, and load coefficient $f_{pm}$ on the dynamic characteristics.
The amplitude-frequency response of the straight bevel gear system exhibits nonlinear jumps and multi-valued solutions in the primary resonance regions near the mesh frequency and shaft frequency. For instance, at $\Omega \in [0.5575, 0.6325]$ and $[1.220, 1.233]$, multiple stable and unstable solutions coexist, indicating amplitude jumps. The backbone curves shift towards lower frequencies, showing stiffness softening behavior. Under small backlash conditions, severe tooth surface impacts occur, while for $b > 0.98$, the jumps and impacts stabilize. The dynamic characteristics are less sensitive to time-varying mesh stiffness, but high-speed light-load conditions and large tooth surface errors exacerbate nonlinear jumps and impacts.
To quantify the tooth surface impact states, we define an impact index $I$ based on the maximum and minimum displacements $\lambda_m \pm \lambda_a$:
$$ I = \begin{cases}
0 & \lambda_m – \lambda_a > b \\
1 & (\lambda_m – \lambda_a \leq b) \cap (\lambda_m + \lambda_a > b) \\
2 & (\lambda_m – \lambda_a < -b) \cap (\lambda_m + \lambda_a > b)
\end{cases} $$
We use the notation “n/I” to represent the number of solutions $n$ and impact state $I$. For example, 2/1 denotes two solutions with single-sided impact, and 2/0-1 indicates two solutions with no impact and single-sided impact coexistence.
The solution domain structure in parameter planes reveals the transition of dynamic behaviors. Below, we summarize the effects of key parameters using tables and formulas.
| Backlash $b$ | Amplitude Jump Region | Tooth Surface Impact | Stability |
|---|---|---|---|
| 0.1 | Near linear, softening | Severe bilateral impact | Unstable |
| 0.2 – 0.6 | Increasing jump region | Transition to unilateral impact | Moderate |
| > 0.98 | Stabilized jumps | Reduced impact | Stable |
The dimensionless mesh stiffness $k_h(\tau)$ is given by:
$$ k_h(\tau) = 1 + \sum_{l=1}^N \frac{A_{kl}}{k_m} \cos(l \Omega \tau + \Phi_{kl}) $$
and the dimensionless force components are:
$$ \begin{aligned}
f_{pm} &= \frac{F_{pm}}{m_e b_h \Omega_n^2} \\
f_{pv} &= \frac{F_{pv}}{m_e b_h \Omega_n^2} \\
f_e &= \sum_{l=1}^N \frac{A_{el}}{b_m} (l \Omega)^2 \cos(l \Omega \tau + \Phi_{el})
\end{aligned} $$
In the $\Omega \times b$ parameter plane, the n/I solution domain structure shows that for non-resonant regions, the system primarily exhibits 1/0 characteristics (single solution, no impact). In primary resonance regions, multiple states like 1/2, 2/1, 2/2, 2/0-1, and 2/0-2 coexist. For shaft frequency resonance, when $b > 0.2$, the n/I characteristics stabilize. For mesh frequency resonance, when $b > 0.6$, the system transitions from 1/2 to 2/0-1, and when $b > 0.98$, 2/0-2 disappears, indicating stabilized jumps and impacts.
| $\alpha$ | Amplitude Jump | Tooth Surface Impact | Resonance Shift |
|---|---|---|---|
| 0.0 – 0.25 | Prominent in mesh frequency | Bilateral impact in shaft frequency | Softening to hardening |
| > 0.4 | Reduced in shaft frequency | Unilateral to no impact | Stable in mesh frequency |
The equivalent mass $m_e$ is defined as:
$$ m_e = \frac{J_1 J_2}{r_1^2 J_2 + r_2^2 J_1} $$
and the dimensionless equations of motion in matrix form are:
$$ \begin{aligned}
\mathbf{M} \ddot{\mathbf{x}} + \mathbf{C} \dot{\mathbf{x}} + \mathbf{K} \mathbf{x} + \mathbf{K}_h \mathbf{f}(\lambda) = \mathbf{F}
\end{aligned} $$
where $\mathbf{x} = [x_1, y_1, z_1, x_2, y_2, z_2, \lambda]^T$ is the displacement vector, $\mathbf{K}_h$ is the nonlinear stiffness matrix due to backlash, and $\mathbf{F}$ is the force vector.
For the $\Omega \times f_e$ parameter plane, as $f_e$ increases, the resonance regions expand and shift to lower frequencies, indicating increased stiffness softening. At $f_e = 0.12$, the shaft frequency region transitions from 1/0 to 1/1, and at $f_e = 0.21$, 2/0-2 appears in the mesh frequency region. At $f_e = 0.34$, 2/2 and 2/1-2 states emerge, accompanied by secondary jumps. High-frequency regions show transitions from 1/0 to 1/1, highlighting the sensitivity to static transmission errors.
| $f_{pm}$ | Amplitude Jump | Tooth Surface Impact | Resonance Behavior |
|---|---|---|---|
| 0.5 – 0.57 | Prominent jumps | Bilateral impact in shaft frequency | Softening |
| 0.75 – 1.0 | Reduced jumps | Unilateral to no impact | Weak linearity |
The harmonic balance equations are solved iteratively. The Jacobian matrix for the Broyden method is updated using:
$$ \mathbf{J}_{k+1} = \mathbf{J}_k + \frac{(\mathbf{y} – \mathbf{J}_k \mathbf{s}) \mathbf{s}^T}{\mathbf{s}^T \mathbf{s}} $$
where $\mathbf{s}$ is the step vector and $\mathbf{y}$ is the residual change. The pseudo-arc-length continuation ensures continuity in the solution curves by parameterizing the arc length $s$:
$$ \left\| \frac{d\mathbf{x}}{ds} \right\|^2 + \left\| \frac{d\Omega}{ds} \right\|^2 = 1 $$
Numerical results validate the n/I characteristics. For example, at $b = 0.1$, the amplitude-frequency response shows multiple solution branches and impact transitions. The mean and amplitude components $\lambda_m$ and $\lambda_a$ are plotted against frequency, confirming the coexistence of 2/2 and 2/1-2 states in jump regions. The solution domain structure maps provide a comprehensive view of parameter sensitivities, aiding in the selection of design parameters to avoid undesirable dynamic behaviors.
In conclusion, the dynamic characteristics of straight bevel gear systems are highly influenced by backlash, time-varying stiffness, static transmission errors, and load conditions. Amplitude jumps and multi-valued solutions occur primarily in the mesh and shaft frequency resonance regions, with the mesh frequency region exhibiting more complex behaviors. Small backlash leads to severe tooth surface impacts, but increasing backlash beyond 0.98 stabilizes the system. Time-varying mesh stiffness has a limited effect, while high-speed light-load conditions and large errors intensify nonlinear phenomena. The solution domain structure in parameter planes offers valuable data for designing straight bevel gear systems with improved stability and reduced impact, ensuring reliable operation in practical applications.
The analysis underscores the importance of considering nonlinear factors in the design process of straight bevel gear systems. Future work could explore additional parameters, such as damping variations and thermal effects, to further enhance the dynamic performance of straight bevel gears.