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 backlash, time-varying stiffness, and static transmission errors often leads to nonlinear dynamic behaviors such as amplitude jumping, multi-valued solutions, and tooth surface impacts. These phenomena can significantly affect the performance, reliability, and lifespan of straight bevel gear systems. This article investigates the dynamic characteristics of a straight bevel gear system with backlash using the harmonic balance method. By analyzing the parameter solution domain structure in two-parameter planes, we explore the sensitivity of amplitude jumping and tooth surface impacts to various parameters, providing data support for the structural design of straight bevel gears.
The dynamic model of a straight bevel gear system considers elastic supports, time-varying meshing stiffness, and backlash. The system consists of two straight bevel gears with orthogonal axes intersecting at point O. A global coordinate system ∑: {Oxyz} is established, with the supports equivalent to stiffness and damping elements at the midpoints of the gear widths. The vibrations of the rigid gear bodies include translations along the three coordinate axes and rotations about the axes, represented by {X1, Y1, Z1, θ1, X2, Y2, Z2, θ2}. The relative torsional displacement along the meshing line, considering vibrations and errors, is given by:
$$ \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) a_3 – e_n(t) $$
where a1 = cosδ1 sinαn, a2 = cosδ1 cosαn, a3 = cosαn, δ1 is the pitch cone angle of the driving gear, αn is the normal pressure angle, r1 and r2 are the base circle radii at the midpoints, and en(t) is the static transmission error expressed as a Fourier series:
$$ e_n(t) = \sum_{l=1}^{N} A_l \cos(l \Omega_h t + \Phi_l) $$
The meshing force and its components along the coordinate directions are:
$$ \begin{cases} F_n = k_h(t) f(\Lambda_n) + c_h \dot{\Lambda}_n \\ F_x = -F_n (\sin \alpha_n \cos \delta_1 + \cos \alpha_n \sin \delta_1) = -a_4 F_n \\ F_y = F_n (\sin \alpha_n \sin \delta_1 – \cos \alpha_n \cos \delta_1) = a_5 F_n \\ F_z = F_n \cos \alpha_n = a_3 F_n \end{cases} $$
Here, kh(t) is the time-varying meshing stiffness, represented as:
$$ k_h(t) = k_m + \sum_{l=1}^{N} A_{kl} \cos(l \Omega_h t + \Phi_{kl}) $$
where km is the average meshing stiffness, Akl is the amplitude of the l-th harmonic, and Φkl is the phase angle. The nonlinear backlash function f(Λn) is defined as:
$$ 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} $$
with 2bh representing the average backlash. The vibration equations for the straight bevel gear system are derived as:
$$ \begin{cases} 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{cases} $$
By introducing the dimensionless relative displacement λ = Λn / bh and other dimensionless parameters, the equations are normalized. The dimensionless vibration equations become:
$$ \begin{cases} \ddot{x}_1 + 2\xi_{x1} \dot{x}_1 + 2a_4 \xi_{h1} \dot{\lambda} + k_{x1} x_1 + k_{h1} f(\lambda) = 0 \\ \ddot{y}_1 + 2\xi_{y1} \dot{y}_1 – 2a_5 \xi_{h1} \dot{\lambda} + k_{y1} y_1 – k_{h1} f(\lambda) = 0 \\ \ddot{z}_1 + 2\xi_{z1} \dot{z}_1 – 2a_3 \xi_{h1} \dot{\lambda} + k_{z1} z_1 – k_{h1} f(\lambda) = 0 \\ \ddot{x}_2 + 2\xi_{x2} \dot{x}_2 – 2a_4 \xi_{h2} \dot{\lambda} + k_{x2} x_2 – k_{h2} f(\lambda) = 0 \\ \ddot{y}_2 + 2\xi_{y2} \dot{y}_2 + 2a_5 \xi_{h2} \dot{\lambda} + k_{y2} y_2 + k_{h2} f(\lambda) = 0 \\ \ddot{z}_2 + 2\xi_{z2} \dot{z}_2 + 2a_3 \xi_{h2} \dot{\lambda} + k_{z2} z_2 + 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 \end{cases} $$
where the dimensionless parameters are defined as: xj = Xj / bh, yj = Yj / bh, zj = Zj / bh, λ = Λn / bh, Ωn = (km / me)^{1/2}, τ = Ωt, Ω = Ωh / Ωn, Ωij = (kij / mj)^{1/2}, ξij = cij / (2mj Ωn), kij = (Ωj / Ωn)^2, ξhj = ch / (2mj Ωn), khj(τ) = kh(τ) / (mj Ωn^2), ξh = ch / (2me Ωn), fpm = Fpm / (me bh Ωn^2), fpv = Fpv / (me bh Ωn^2), and fe = Σ Ael / bm (lΩ)^2 cos(lΩτ + Φel). The time-varying stiffness is kh = 1 + Σ (Akl / km) cos(lΩτ + Φkl).
The harmonic balance method is employed to solve the nonlinear dynamic equations. The steady-state response is approximated by a first-harmonic solution:
$$ x_l = x_{ml} + x_{cl} \cos(\Omega \tau) + x_{sl} \sin(\Omega \tau) $$
where xml, xcl, xsl, and xal = (xcl^2 + xsl^2)^{1/2} are the mean, cosine, sine, and amplitude components of the l-th dimension vector. The nonlinear backlash function f(xl, bl) is expanded using Fourier series as:
$$ 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 Nml and Nal given by:
$$ \begin{cases} 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{cases} $$
where μl± = (±bl – xml) / xal. The functions G(μ) and H(μ) are defined as:
$$ G(\mu) = \begin{cases} 2 (\mu \arcsin \mu + \sqrt{1 – \mu^2}) / \pi, & |\mu| \leq 1 \\ |\mu|, & |\mu| > 1 \end{cases} $$
$$ H(\mu) = \begin{cases} 2 (\arcsin \mu + \mu \sqrt{1 – \mu^2}) / \pi, & |\mu| \leq 1 \\ \text{sign}(\mu), & |\mu| > 1 \end{cases} $$
Substituting the harmonic solutions into the dimensionless equations and equating coefficients of like harmonics yields a system of nonlinear algebraic equations:
$$ \begin{cases} K_m G_m – F_m = 0 \\ -\Omega^2 M X_c + \Omega C X_s + K_m G_c + K_c G_m – F_c = 0 \\ -\Omega^2 M X_s – \Omega C X_c + K_m G_s + K_s G_m – F_s = 0 \end{cases} $$
where Gm = {Nmi xmi}7×1, Gc = {Nai xci}7×1, Gs = {Nai xsi}7×1, and the matrices Km, Kc, Ks, C, M, Fm, Fc, Fs are constructed from the system parameters. This system is solved using the Broyden quasi-Newton method and the pseudo-arc length continuation algorithm to obtain the amplitude-frequency response and multi-valued solutions.
Numerical simulations are conducted with dimensionless parameters: ξi1 = ξi2 = 0.01, ξh1 = ξh2 = 0.0125, ξh = 0.05, ki1 = ki2 = 1.0 (i = x, y, z), kh1 = kh2 = 0.5, kh = 1 + α cos(Ωτ), fpm = 0.5, fpv = 0.0, fe = 0.2 Ω^2 cos(Ωτ), α = 0.2, and b = 1.0. The frequency range Ω ∈ [0.05, 2.00] is considered, with parameters b, α, fe, and fpm varied to analyze their effects on dynamic characteristics.
The amplitude-frequency response of the straight bevel gear system exhibits nonlinear jumping and multi-valued solutions in the main resonance regions near the meshing frequency (Ω ≈ 0.68) and the shaft frequency (Ω ≈ 1.240). For instance, with fpm = 0.5, amplitude jumping occurs in intervals [0.5575, 0.6325] and [1.220, 1.233], where three solutions (two stable and one unstable) coexist. The backbone curves shift toward lower frequencies, indicating stiffness softening. When fe = 0.3, double jumps appear near the meshing frequency, and the backbone shows both hardening and softening behaviors.
The sensitivity of amplitude jumping to parameters is analyzed. For backlash b, smaller values (e.g., b = 0.1) result in nearly linear responses with mild stiffness softening. As b increases, softening intensifies, and the main resonance shifts to lower frequencies. The jumping region near the meshing frequency expands, while the shaft frequency region remains relatively insensitive for b > 0.3. For time-varying stiffness amplitude α, increasing α causes secondary jumps near the meshing frequency and transitions from softening to hardening. In the shaft frequency region, jumping weakens, suggesting that reduced contact ratio exacerbates nonlinearities. For static error fe, larger values intensify jumping in both resonance regions, with stiffness softening transitioning to hardening at the meshing frequency. Improved machining accuracy can mitigate these effects. For load fpm, higher loads reduce nonlinear jumping and softening, leading to weaker linear behavior, consistent with reduced tooth surface impacts.
Tooth surface impacts are classified into no impact (I=0), single-sided impact (I=1), and double-sided impact (I=2) based on the maximum and minimum displacements λm ± λ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} $$
The multi-valued solution and impact state are denoted as “n/I”, where n is the number of solutions (1 or 2), and I is the impact state. For example, 2/1 indicates two solutions with single-sided impact, and 2/1-0 represents two solutions with mixed no-impact and single-sided impact. In the main resonance regions, multiple states coexist, such as 2/0-1, and transitions between states can lead to complex dynamics, including chaos.
The parameter solution domain structure is analyzed in two-parameter planes: Ω × b, Ω × α, Ω × fe, and Ω × fpm. Each plane uses 501 × 801 parameter points, with colors indicating different n/I states. The distribution reveals the evolution of dynamic characteristics with parameters.
| Parameter | Effect on Meshing Frequency Region | Effect on Shaft Frequency Region | Critical Values |
|---|---|---|---|
| Backlash (b) | Jumping intensifies, softening increases, transitions to 2/0-1 at b>0.6, stable at b>0.98 | Stable for b>0.2, n/I states weaken | b=0.98 for stability |
| Stiffness Amplitude (α) | Region expands, complex n/I states | Region shrinks, jumping weakens, 2/1-2 disappears at α>0.25, 1/1 to 1/0 at α>0.4 | α=0.4 for linearity |
| Static Error (fe) | Region expands and shifts left, stiffness softening, 2/0-2 at fe=0.21, 2/2 and 2/1-2 at fe=0.34 | Region expands, 1/0 to 1/1 at fe=0.12, 2/0-1 at fe=0.19 | fe=0.34 for severe jumping |
| Load (fpm) | Region shifts right, jumping weakens, 2/0-1 disappears | Region shifts right, 2/1-2 disappears at fpm=0.57, 1/1 to 1/0 at fpm=0.75 | fpm=0.75 for reduced impact |
In the Ω × b plane, non-resonant regions are dominated by 1/0 states. In main resonance regions, 1/2, 2/1, 2/2, 2/0-1, and 2/0-2 states coexist. For the shaft frequency, stability is achieved at b > 0.2. For the meshing frequency, transitions occur at b > 0.6 (to 2/0-1) and b > 0.98 (2/0-2消失), where jumping and impacts stabilize. Increasing backlash suppresses tooth surface impacts.
In the Ω × α plane, the shaft frequency region shrinks with increasing α, and n/I states become less sensitive. At α > 0.25, 2/1-2 disappears, and at α > 0.4, 1/1 transitions to 1/0. The meshing frequency region expands with stable n/I characteristics.
In the Ω × fe plane, both resonance regions expand and shift leftward with increasing fe, indicating stiffness softening. At the shaft frequency, transitions occur at fe = 0.12 (1/0 to 1/1) and fe = 0.19 (2/0-1). At the meshing frequency, 2/0-2 appears at fe = 0.21, and 2/2 and 2/1-2 at fe = 0.34, accompanied by secondary jumps. High-frequency regions show increased impacts due to error excitation amplification.
In the Ω × fpm plane, resonance regions shift rightward with increasing fpm, transitioning from softening to weak linearity. At the shaft frequency, 2/1-2 disappears at fpm = 0.57, and 1/1 at fpm = 0.75. At the meshing frequency, 2/0-1 vanishes. High-speed, light-load conditions exacerbate nonlinearities and impacts.
To validate the n/I characteristics, a section at b = 0.1 is analyzed. The mean λm, amplitude λa, and λm ± λa distributions clearly show amplitude jumping, multi-valued solutions, and impact state transitions. In double-solution regions, multiple impact states (e.g., 2/2 and 2/1-2) coexist, consistent with the solution domain structure.
In conclusion, the dynamic characteristics of straight bevel gear systems are highly sensitive to parameters such as backlash, time-varying stiffness, static error, and load. Nonlinear jumping and tooth surface impacts are prominent in main resonance regions, especially at the meshing frequency. Small backlash leads to severe impacts, while larger backlash (b > 0.98) stabilizes the system. Time-varying stiffness has limited sensitivity, but static error and light loads intensify nonlinearities. The parameter solution domain structure provides a comprehensive view of n/I state transitions, aiding in the design of straight bevel gears to avoid undesirable dynamic behaviors. For reliable operation, parameters should be selected from regions with single-valued, no-impact states (1/0). This analysis underscores the importance of considering nonlinear dynamics in the design and optimization of straight bevel gear systems.