Nonlinear Frequency Response Analysis of Herringbone Gears Using Incremental Harmonic Balance Method

In this study, we investigate the nonlinear frequency response characteristics of herringbone gear pairs, which are widely used in heavy machinery such as vehicles and ships due to their high load capacity, compact structure, and balanced axial forces. The dynamic behavior of herringbone gears is critical for system stability and reliability, especially under nonlinear effects like time-varying mesh stiffness, backlash, and excitation forces. Previous research has extensively explored spur gears and planetary systems, but detailed analyses of herringbone gears’ nonlinear dynamics remain limited. Here, we develop a torsional dynamic model incorporating key nonlinear factors and apply the Incremental Harmonic Balance Method (IHBM) to analyze frequency responses. We validate results with numerical methods and examine parameter influences to provide insights for vibration control and optimization. Throughout this article, the term ‘herringbone gears’ is emphasized to highlight our focus on this specific gear type.

We begin by establishing a nonlinear torsional dynamic model for a herringbone gear pair. The model considers time-varying mesh stiffness, constant and dynamic backlash, static transmission error, and external dynamic excitations. The system is represented as a four-degree-of-freedom model, as shown in the following schematic, which illustrates the torsional dynamics involved.

The equations of motion are derived using Newton’s second law. Let \(I_{pL}\), \(I_{pR}\), \(I_{gL}\), and \(I_{gR}\) represent the polar moments of inertia at the equivalent nodes for the left and right helical gears of the pinion and gear, respectively. The angular displacements are denoted by \(\theta_{pL}\), \(\theta_{pR}\), \(\theta_{gL}\), and \(\theta_{gR}\). The base circle radii are \(r_{pL}\), \(r_{pR}\), \(r_{gL}\), and \(r_{gR}\), and the helix angles are \(\beta_L\) and \(\beta_R\). The mesh stiffness and damping for the left and right gear pairs are \(k_{mL}\), \(k_{mR}\), \(c_{mL}\), and \(c_{mR}\), while the torsional stiffness and damping of the intermediate shaft segments are \(k_{pt}\), \(k_{gt}\), \(c_{pt}\), and \(c_{gt}\). The relative mesh displacements \(\delta_L\) and \(\delta_R\) are defined as:

$$ \delta_L = (r_{pL}\theta_{pL} – r_{gL}\theta_{gL})\cos\beta_L – e_L(t) $$
$$ \delta_R = (r_{pR}\theta_{pR} – r_{gR}\theta_{gR})\cos\beta_R – e_R(t) $$

where \(e_L(t)\) and \(e_R(t)\) are static transmission errors. The relative torsional displacements between coaxial gears are:

$$ \gamma_P = \theta_{pL} – \theta_{pR}, \quad \gamma_G = \theta_{gL} – \theta_{gG} $$

The backlash function \(f(\delta_j)\) for \(j = L, R\) is piecewise-linear:

$$ f(\delta_j) =
\begin{cases}
\delta_j – b, & \delta_j > b \\
0, & -b \leq \delta_j \leq b \\
\delta_j + b, & \delta_j < -b
\end{cases} $$

Here, \(b\) is half the backlash. We consider two backlash models: constant backlash \(b_1 = b_0(1 + W)\) and dynamic backlash \(b_2 = b_0\left(1 + \sum_{l=1}^L b_l \cos(l\omega t)\right)\), where \(b_0\) is the mean backlash, \(W\) accounts for wear, and \(b_l\) are dynamic backlash amplitudes. The equations of motion in matrix form are:

$$ \tilde{M}\tilde{q}”(\bar{t}) + \tilde{C}\tilde{q}'(\bar{t}) + \tilde{K}\tilde{q}(\bar{t}) = \tilde{F}(\bar{t}) $$

where \(\tilde{M}\), \(\tilde{C}\), and \(\tilde{K}\) are dimensionless mass, damping, and stiffness matrices, \(\tilde{q}\) is the displacement vector, and \(\tilde{F}\) is the excitation vector. The dimensionless parameters are defined using a nominal displacement scale \(b_c = b_0\), natural frequency \(\omega_n = \sqrt{k_{\text{mean}}/m_c}\), and normalized time \(\bar{t} = \omega_n t\). The dimensionless stiffness matrix \(\tilde{K}\) is:

$$ \tilde{K} = \begin{bmatrix}
k_{11} & 0 & k_{13} & k_{14} \\
0 & k_{22} & k_{23} & k_{24} \\
k_{31} & k_{32} & k_{33} & 0 \\
k_{41} & k_{42} & 0 & k_{44}
\end{bmatrix} $$

with elements such as \(k_{11} = \frac{\cos^2\beta_L k_{mL}}{m_c \omega_n^2}\), \(k_{13} = \frac{r_{pL}^2 k_{pt} \cos\beta_L}{I_{pL} \omega_n^2}\), and others derived from gear geometry. Time-varying mesh stiffness and static transmission error are expressed as Fourier series:

$$ k_{mj}(t) = k_{\text{mean}}\left[1 – \sum_{l=1}^L \epsilon_l \cos(l\omega t)\right], \quad j = L, R $$
$$ e_j(t) = e_m + \sum_{l=1}^L e_l \sin(l\omega t), \quad j = L, R $$

where \(\epsilon_l\) and \(e_l\) are fluctuation coefficients. Damping is modeled as \(c_{mj} = 2\xi \sqrt{m_c k_{\text{mean}}}\) for mesh damping and \(c_{it} = 2\xi \sqrt{m_c k_{it}}\) for torsional damping, with \(\xi\) as the damping ratio.

To analyze the frequency response, we apply the Incremental Harmonic Balance Method (IHBM). This method is effective for solving nonlinear differential equations with periodic coefficients. We set \(\tau = \Omega \bar{t}\), where \(\Omega = \omega/\omega_n\) is the normalized mesh frequency. The equation becomes:

$$ \Omega^2 \tilde{M} q”(\tau) + \Omega \tilde{C} q'(\tau) + \tilde{K} q(\tau) = \tilde{F}(\tau) $$

We seek a periodic solution of the form \(q(\tau) = q_0(\tau) + \Delta q(\tau)\), where \(q_0(\tau)\) is an initial approximation and \(\Delta q(\tau)\) is an increment. The solution is expanded in Fourier series:

$$ q_{j0} = \sum_{n=0}^N \left[a_{jn} \cos(n\tau) + b_{jn} \sin(n\tau)\right] $$
$$ \Delta q_j = \sum_{n=0}^N \left[\Delta a_{jn} \cos(n\tau) + \Delta b_{jn} \sin(n\tau)\right] $$

for \(j = 1, 2, 3, 4\), representing the four degrees of freedom. The backlash function is linearized as \(f(\Delta_j) = f(\Delta_{j0}) + f'(\Delta_{j0})\Delta(\Delta_j)\). Substituting into the equation and applying the Galerkin procedure, we obtain a linear system for the Fourier coefficients:

$$ C \Delta A = R $$

where \(\Delta A\) is the vector of incremental coefficients, \(C\) is a matrix derived from the system parameters, and \(R\) is the residual vector. This system is solved iteratively using the Newton-Raphson method until convergence. The IHBM allows us to capture harmonic responses accurately, including super-harmonic components due to nonlinearities.

We now present the frequency response characteristics of herringbone gears under various parameters. The base system parameters are: \(L = 3\), \(N = 11\), \(F_1 = 1.3534\), \(\xi = 0.045\), \(\epsilon_1 = 0.16\), \(\epsilon_2 = 0.07\), \(\epsilon_3 = 0.03\), \(e_1 = 0.03\), \(e_2 = 0.01\), \(e_3 = 0.01\), \(f_1 = 0.06\), \(f_2 = 0.03\), \(f_3 = 0.01\), \(W = 0\), \(b_1 = 0.05\), \(b_2 = 0.02\), \(b_3 = 0.01\). We compute amplitude-frequency curves for both constant and dynamic backlash models and validate with Runge-Kutta numerical integration. The results show that herringbone gears exhibit nonlinear phenomena such as jump discontinuities and multi-valued solutions, influenced by parameters like stiffness, damping, and excitations.

The following table summarizes the effects of key parameters on the amplitude-frequency response of herringbone gears, based on our analysis using IHBM:

Parameter Effect on Amplitude-Frequency Response Influence on Nonlinear Phenomena
Time-Varying Mesh Stiffness Increases amplitude; reduces with lower fluctuation coefficients Promotes jumps and multi-valued solutions; reduces with decreased fluctuations
Damping Ratio Decreases amplitude; suppresses vibrations Reduces jumps and multi-valued intervals; stabilizes response
Static Transmission Error Increases amplitude; larger errors amplify response Enhances jumps; not a primary source of nonlinearity
External Excitation Amplitude Increases amplitude slightly; minimal effect on state changes Little impact on jump frequencies; multi-valued regions may shrink
Constant Backlash Increases amplitude with larger backlash Widens multi-valued intervals; jumps persist
Dynamic Backlash Slightly increases amplitude; reduces multi-valued regions Narrows or eliminates multi-valued intervals; controls impacts

To quantify the frequency response, we define the normalized amplitude \(A = \max(|\Delta_j|)\) over a period. For a given \(\Omega\), the response can include primary resonance at \(\Omega = 1\) and super-harmonic resonances at \(\Omega = 1/2, 1/3\), etc., due to higher harmonics in excitations. The equation for the dimensionless amplitude in terms of system parameters can be expressed as:

$$ A(\Omega) = \frac{\sqrt{\left(\sum_{l} F_l \cos(l\Omega\tau)\right)^2 + \left(\sum_{l} (l\Omega)^2 e_l \sin(l\Omega\tau)\right)^2}}{\sqrt{\left(\Omega^2 \tilde{M} – \tilde{K}_{\text{eff}}\right)^2 + \left(\Omega \tilde{C}\right)^2}} $$

where \(\tilde{K}_{\text{eff}}\) is the effective stiffness including nonlinear terms. This formula highlights how excitations and stiffness variations drive the response. For herringbone gears, the coupling between left and right gear pairs adds complexity, but the IHBM effectively decomposes it into harmonic components.

We further analyze the impact of backlash models. For constant backlash, variations in \(W\) from -0.6 to 0.6 show that increased backlash raises vibration amplitudes and expands multi-valued regions. For dynamic backlash, increasing the amplitudes \(b_l\) from zero to 0.15 leads to a slight amplitude rise but shrinks multi-valued intervals, indicating better control over nonlinear vibrations. This is summarized in the table below for specific cases:

Backlash Model Parameter Values Amplitude Trend Multi-Valued Interval (\(\Omega\) range)
Constant (\(b_1\)) \(W = -0.6, -0.3, 0, 0.3, 0.6\) Increases with \(W\) Widens from 0.78-0.86 to broader ranges
Dynamic (\(b_2\)) \(b_1 = 0, 0.05, 0.10, 0.15; b_2, b_3\) proportional Slight increase Narrows from 0.83-0.84 to disappearance

The role of time-varying mesh stiffness is critical in herringbone gears. We test three cases: (i) \(\epsilon_1 = 0.16, \epsilon_2 = 0.07, \epsilon_3 = 0.03\), (ii) \(\epsilon_1 = 0.13, \epsilon_2 = 0.04, \epsilon_3 = 0.01\), and (iii) constant stiffness (\(\epsilon_l = 0\)). The amplitude-frequency curves demonstrate that reducing stiffness fluctuations decreases amplitudes, especially at lower frequencies, and mitigates jumps. In dynamic backlash models, this reduction is more effective for vibration control. The mathematical representation of stiffness influence is embedded in the \(\tilde{K}\) matrix, where terms like \(k_{11}\) depend on \(k_{mL}\), which varies with \(\epsilon_l\).

Damping effects are examined for \(\xi = 0.03, 0.045, 0.06\). Higher damping ratios reduce amplitudes, suppress jump phenomena, and narrow multi-valued regions. This aligns with the damping term \(\Omega \tilde{C} q'(\tau)\) in the equation, which dissipates energy. For herringbone gears, increased damping is particularly beneficial in dynamic backlash scenarios to stabilize responses.

Static transmission error variations are studied with \(e_1 = 0.1, 0.03, 0\) (and proportional \(e_2, e_3\)). Larger errors amplify amplitudes and enhance jumps, but even zero error still yields nonlinear behavior, confirming that other factors like backlash dominate nonlinearity in herringbone gears. The error term appears in \(\tilde{F}(\tau)\) as \(-e”_j(\tau)\), contributing to excitation forces.

External excitation amplitude, represented by \(F_1\) values of 0.8534, 1.3534, and 3.3534, shows that higher excitations increase amplitudes but have minimal effect on jump frequencies. Multi-valued intervals may decrease with higher \(F_1\), suggesting that load changes do not alter the fundamental nonlinear state of herringbone gears significantly. This is captured in the forcing vector \(\tilde{F}(\tau)\), which includes terms like \(F_1 \cos \beta_L\).

To elaborate on the IHBM solution process, we derive the linearized equation for increments. After Taylor expansion and Galerkin projection, the coefficients matrix \(C\) has elements computed from integrals over \([0, 2\pi]\). For example, for each harmonic \(i\), we have:

$$ C_{ij} = \int_0^{2\pi} \left[\Omega^2 \tilde{M} \phi_i” \phi_j” + \Omega \tilde{C} \phi_i’ \phi_j’ + \tilde{K}_{\text{eff}} \phi_i \phi_j\right] d\tau $$

where \(\phi_i\) are basis functions (\(\cos(i\tau)\) or \(\sin(i\tau)\)). The residual vector \(R\) involves the mismatch between the initial guess and the equation. Iterative updates \(\Delta A = C^{-1} R\) converge to a precise solution. This method efficiently handles the piecewise nonlinearity of herringbone gears, avoiding numerical instability seen in pure integration methods at high frequencies (e.g., \(\Omega \approx 1.8\)).

In terms of super-harmonic responses, we observe that for \(\Omega = 1/2\) and \(\Omega = 1/3\), the system exhibits periodic motions with multiple harmonics, evident in phase plots that are non-circular closed curves. These are driven by higher-order terms in the stiffness and excitation series. The general solution form includes contributions from all harmonics:

$$ q(\tau) = \sum_{n=0}^N \left[A_n \cos(n\Omega\tau) + B_n \sin(n\Omega\tau)\right] $$

where \(A_n\) and \(B_n\) are determined by IHBM. This reveals that herringbone gears can experience complex vibrations beyond primary resonance, impacting noise and wear.

For design implications, our analysis suggests that dynamic backlash models offer better vibration control for herringbone gears compared to constant backlash, as they reduce multi-valued regions. Adjusting stiffness fluctuations and damping can further optimize performance. The equations and tables provided here serve as a reference for engineers working with herringbone gear systems.

In conclusion, our study using the Incremental Harmonic Balance Method elucidates the nonlinear frequency response of herringbone gears. Key findings include the presence of primary and super-harmonic resonances, the influence of parameters on amplitude and jumps, and the advantage of dynamic backlash for vibration control. This work enhances understanding of herringbone gear dynamics and aids in developing quieter, more reliable transmission systems. Future research could extend to multi-mesh herringbone gear trains or experimental validation.

Scroll to Top