# Effect of gear topology discontinuities on the nonlinear dynamic response of a multi-degree-of-freedom gear train

Adrien Mélot, Youness Benaïcha, Emmanuel Rigaud, Joël Perret-Liaudet, Fabrice Thouverez

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

Adrien Mélot, Youness Benaïcha, Emmanuel Rigaud, Joël Perret-Liaudet, Fabrice Thouverez. Effect of gear topology discontinuities on the nonlinear dynamic response of a multi-degree-of-freedom gear train. Journal of Sound and Vibration, 2021, 10.1016/j.jsv.2021.116495. hal-03233115

HAL Id: hal-03233115

https://hal.science/hal-03233115

Submitted on 23 May 2021

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>
# Effect of gear topology discontinuities on the nonlinear dynamic response of a multi-degree-of-freedom gear train

Adrien Mélot<sup>a,*</sup>, Youness Benaïcha<sup>a,b</sup>, Emmanuel Rigaud<sup>a</sup>, Joël Perret-Liaudet<sup>a</sup> and Fabrice Thouverez<sup>a</sup>

## ARTICLE INFO

### Keywords:

Computational nonlinear dynamics

Vibro-impacts

Harmonic Balance Method

Bifurcation analysis

Gear backlash

Static transmission error

## ABSTRACT

Weight reduction is a recurring concern in the design of modern mechanical systems. This search may lead engineers to resort to using gears with holes in order to meet their requirements. This paper presents a methodology to carry out nonlinear dynamic analyses of a gear transmission with holes in the gear blanks subjected to a multiharmonic internal excitation. This work investigates the influence of these holes on the vibration levels, occurrence of contact loss and possible bifurcations. The numerical model features two flexible shafts coupled by a spur gear pair with holes. The gear model consists in two lumped masses and inertias and includes the gear backlash as well as the internal excitation sources that are the time-varying stiffness and the static transmission error (STE). The resulting mechanical system is solved in the frequency domain by the Harmonic Balance Method (HBM) coupled with an arc-length continuation algorithm and its stability is evaluated with Hill's method. Results show that adding holes not only impacts the STE but also the mesh stiffness. The interactions of these two quantities has a substantial influence on the bifurcation structure along the main solution branch and leads to a decrease of the span of vibro-impact regions. As the applied static torque is increased, it is found that using holed gear blanks can effectively prevent contact loss and lead to a linear response.

## 1. Introduction

Geared systems are known to be prone to unwanted noise and vibrations stemming from the meshing process and possible contact loss between gear teeth. These vibrations have motivated many studies over the last few decades. Such systems are characterised by an internal excitation source known as the static transmission error (STE) whose origin lies in the teeth deflection under load, manufacturing defects and tooth profile modifications. The STE is defined as the difference between the actual position of the output gear and the position it would occupy if the gear pair were perfectly conjugate [59] and the mesh stiffness is expressed as the derivative of the transmitted load relative to the STE [8, 13].

It is well-established that the STE, together with mesh stiffness variations, are responsible for parametric resonances. The induced dynamic mesh forces are transmitted to the housing through the shafts and the bearings. The resulting vibrations are at the origin of the so-called whining noise [8, 37]. However, under particular loading conditions, tooth separation can occur due to the gear backlash necessary to allow for assembly and operation [1, 2, 26, 33]. These events induce nonlinear vibrations which generate a noise that can be classified into two categories depending on the loading scenario: rattling noise in case of lightly loaded gears [27, 46, 25] and hammering noise in case of a high static torque [43].

Most studies concerning the nonlinear behaviour of gears are limited to single degree of freedom (DoF) models [11, 50]. A few consider the bearings flexibility [38, 45] but nonlinear analyses of transmission models including the shafts flexibility are not always carried out [32, 56]. This is particularly restrictive since the number of considered degrees-of-freedom directly affects the stiffness distribution [47] and therefore the strain energy at the mesh which can result in an inaccurate prediction of contact loss thresholds.

A wide variety of excitations have been considered in the literature for the analysis of the nonlinear dynamic response. The response to the STE has naturally been studied [63, 61]. The response to torque fluctuations was

<figure><img src="image_3.png" /></figure>

<footnote><sup>a</sup>Laboratoire de Tribologie et Dynamique des Systèmes, UMR CNRS 5513, Ecole Centrale de Lyon, 36 avenue Guy de Collongue, 69134 Ecully Cedex, France</footnote>

<footnote><sup>b</sup>ANSYS SAS, 35-37 rue Louis Guérin, 69100 Villeurbanne Cedex, France</footnote>

<footnote>*Corresponding author</footnote>

<footnote>adrien.melot@ec-lyon.fr (A. Mélot)</footnote>

<footnote>ORCID(s):</footnote>
also investigated using both harmonic [25, 27, 57] and multiharmonic excitations [4, 24]. Several studies were also carried out considering a stochastic excitation [29, 58], variable backlash [41, 48] and variable tooth profile deviations [39]. Nonetheless, very few studies deal with the impact of holes in the gear webs. Guilbert [17, 18] et al. developed a mixed lumped-parameter/finite element model to investigate how holes influence the static and transient linear dynamic response of a gear pair. This model, although able to deal with highly flexible gears, is to this day limited to linear dynamic analyses and requires a large computational effort. Shweiki [53] et al. used a FE simulation to compute the STE and a single degree-of-freedom torsional model to study the impact of two gear blank topologies (circular and oblong holes) on the dynamics solved by time integration. However, the parametric study they proposed relies on an analytical approximation of the STE by arbitrarily changing its mean and oscillating components. This could represent a sensitivity analysis to the applied torque but it lacks a fully fledged sensitivity analysis to the gear blank topology. The aim of this study is thus to investigate the effect of holes on the nonlinear response and bifurcation structure of a simple, yet representative, transmission model including bearing and shaft flexibility.

The paper is structured as follows: first, we introduce the dynamic model of the whole transmission and the computational procedure used to compute the static transmission error and mesh stiffness. In section 3, we give a brief theoretical background on the harmonic balance method and continuation techniques employed to build the nonlinear forced response curves and carry out the stability analyses. Finally, a thorough analysis of the presence of holes and key design parameters and their influence on the nonlinear dynamic response is presented and conclusions are drawn.

## 2. Flexible transmission model

### 2.1. Dynamic model

Figure 1 depicts the model of a kinematic chain equipped with a motor which drives two counter-rotating shafts using a reserve gear pair. The nonlinear dynamic behaviour of this system is the main focus of this work. The shafts are modelled by torsion springs of stiffness $K_{1,2}$ and are supported by two bearings with translational stiffness $K_b$. On one end of both shafts, the input and the output are modelled by two lumped inertias, respectively denoted $I_{in}$ for the input, corresponding to the motor inertia, and $I_{out}$ for the output. The shafts are coupled by a reverse spur gear pair whose characteristics are reported in table 2. The gears are modelled as lumped masses and inertias denoted respectively $m_1$ and $I_1$ for the input gear and $m_2$ and $I_2$ for the output gear.

The behaviour of a gear is modelled with a single node and two degrees-of-freedom at the centre of each gear: one in the torsional direction and one in the translational direction along the line of action. This assumption is reasonable since the tooth width is supposed much smaller than the length of the shafts so that any deformation that would induce a twist of the gear pair is negligible. In the general case the nonlinear force reads

$$ \mathbf{f}_{nl}(\mathbf{q}) = \sum_{e} \mathbf{G}_{e} F_{nl,e} (\mathbf{G}_{e}^{T} \mathbf{q}) \quad (1) $$

where $\mathbf{q} = (\theta_{in}, \theta_{out}, \theta_1, u_1, \theta_2, u_2)^T$ is the vector of generalised coordinates, $e$ denotes the nonlinear element number and $\mathbf{G}_e$ is the force direction associated to the element. In our case, the gears are connected by a nonlinear element consisting of a time-varying piecewise linear stiffness. It includes the backlash $2b$ and static transmission error as a gap function $g(t)$. The nonlinear mesh force $f_{nl}$ can be written as follows:

$$ \mathbf{f}_{nl}(\mathbf{q}) = \mathbf{G} k_m(t) (\mathbf{G}^T \mathbf{q} - g(t)) \mathcal{H} (\mathbf{G}^T \mathbf{q} - g(t)) + \mathbf{G} k_m(t) (\mathbf{G}^T \mathbf{q} + g(t)) \mathcal{H} (-\mathbf{G}^T \mathbf{q} - g(t)) \quad (2) $$

Here, $\mathcal{H}$ is the Heaviside step function and $\mathbf{G}$ is a $(6 \times 1)$ column vector allowing for the projection of the displacements in the global reference frame on the line of action (see Fig. 2). Only the restriction to the nonlinear degrees-of-freedom has non-zero values. It is expressed as

$$ \mathbf{G} = (0 \ 0 \ r_{b1} \tan(\beta) \ 1 \ r_{b2} \tan(\beta) \ -1)^T. \quad (3) $$

The expression of the gap function, as depicted in Fig. 3, can be deduced from the knowledge of the loaded STE. Indeed, the STE $q_s(t)$ is defined as the static deformation to the input torque $T$. By considering a piecewise linear stiffness, contact between the active tooth flanks (positive gap) necessarily occurs at

$$ g(t) = b + q_s(t) - \frac{F_s}{k_m(t)} \quad (4) $$

where $F_s = T/r_b$ corresponds to the transmitted load.
<figure><img src="image_4.png" /><figcaption>**Figure 1:** Dynamic model of the flexible transmission considered in this study.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 2:** Local reference frame of the gear pair.</figcaption></figure>

## 2.2. Computation of the static transmission error with a nonlinear quasi-static analysis

The static transmission error and the mesh stiffness are computed by means of the decomposition approach introduced by Benaïcha et al. in [5], which is based on substructuring and a nonlinear quasi-static simulation. The goal is to obtain these important quantities without heavy computational effort in comparison with a full nonlinear multibody analysis while retaining a satisfying accuracy. To this end, the gear pair with holes is separated in two substructures: one consisting of the gear pair without holes denoted 'wh' and a second one, labelled 'gb', comprising solely the gear body with holes.

More precisely, substructure 'wh' is a gear pair where the gear body is replaced by a rigid disk in order to compute the static deflection of the gear teeth only. A nonlinear quasi-static analysis is performed over a mesh period to evaluate the static deflection along the line of action. The contact is described by the Signorini conditions [60] and an augmented lagrangian formalism [30] is used for its numerical treatment.

The hole-induced deflection is computed with substructure 'gb'. Indeed, the torsional displacement of the holed gear blank is computed off-line of the multibody analysis by enforcing instantaneous tangential displacements on the DoFs located on the edge of the gear blank. These displacements are determined by discretizing the angle $2\pi/N_h$
<figure><img src="image_6.png" /><figcaption>**Figure 3:** Nonlinear mesh force (dashed blue line) and contact model of the gear pair excited by the STE with a time-varying stiffness (solid red line).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 4:** Flowchart of the decomposition procedure</figcaption></figure>

into 16 intermediate positions and performing a static simulation of the full multibody model at each angular position. This allows for a significant reduction of the computation time since the step size is much larger than that needed with classical multibody approaches [6].

The contributions of the different substructures are then summed along the line of action in order to reconstruct the static transmission error of the holed gear pair (see Fig.4). The mesh stiffness is determined by numerical differentiation of the transmitted load relative to the STE [8, 13].
## 3. Numerical procedure for the nonlinear dynamic analysis

The computational strategy developed hereafter is based on a frequency domain technique, namely the Harmonic Balance Method (HBM). The method has been used by the nonlinear dynamics community for several decades [55] in a wide variety of fields such as bolted structures [23], turbomachinery [49], bladed disks dynamics [12, 44], vibro-impact systems [3] and gear dynamics [1, 2, 62].

### 3.1. Harmonic Balance Method

The equation of motion of the previously described system takes the general form

$$[\mathbf{M}]\ddot{\mathbf{q}} + [\mathbf{C}]\dot{\mathbf{q}} + [\mathbf{K}]\mathbf{q} + \mathbf{f}_{nl}(\mathbf{q}) = \mathbf{f}_{ex} \quad (5)$$

where $\mathbf{q}$ contains the generalised displacement of each DoF and $[\mathbf{M}]$, $[\mathbf{C}]$, $[\mathbf{K}]$ are respectively the mass, damping and stiffness matrices. $\mathbf{f}_{ex}$ is the vector of external periodic forcing and $\mathbf{f}_{nl}$ the vector of nonlinear forces, i.e. the mesh force caused by the contact, or lack thereof, between gear teeth. The following study is interested in periodic solutions of fundamental frequency $\Omega$ which, in our case, corresponds to the angular frequency of both shafts. Note that the proposed methodology can be used with all sorts of gear ratios, the only difference would be the fundamental frequency. The terms in equation (5) are thus expanded as truncated Fourier series. More specifically, the generalised displacements are expressed as

$$\mathbf{q} \approx \sum_{k=0}^{H} \mathbf{a}_k \cos(k\Omega t) + \mathbf{b}_k \sin(k\Omega t) = [\mathbf{T} \otimes [I_n]]\tilde{\mathbf{q}} \quad (6)$$

where $\mathbf{T}$ contains the harmonic base functions up to the truncation order $H$, $\otimes$ is the Kronecker product, $[I_n]$ the identity matrix of size $n$ and $\tilde{\mathbf{q}}$ is the vector in which the Fourier coefficients are stored. The velocities and accelerations are expressed by means of the frequency domain differential operator $[\nabla]$ [10]

$$\dot{\mathbf{q}} = [\mathbf{T}[\nabla] \otimes [I_n]]\tilde{\mathbf{q}} \quad (7)$$

$$\ddot{\mathbf{q}} = [\mathbf{T}[\nabla]]^2 \otimes [I_n] \tilde{\mathbf{q}} \quad (8)$$

Substituting these expressions into the equation of motion (5) and applying a Galerkin procedure with the same harmonic base functions $\mathbf{T}$ yields the nonlinear algebraic problem to be solved

$$\mathbf{R}(\tilde{\mathbf{q}}, \Omega) = [\mathbf{Z}(\Omega)] \tilde{\mathbf{q}} + \tilde{\mathbf{f}}_{nl}(\tilde{\mathbf{q}}) - \tilde{\mathbf{f}}_{ex} = \mathbf{0} \quad (9)$$

where $\tilde{\mathbf{f}}_{nl}$ and $\tilde{\mathbf{f}}_{ex}$ contain the Fourier coefficients of the nonlinear forces and external forcing, respectively, and $[\mathbf{Z}(\Omega)]$ is the dynamic stiffness matrix duplicated on each considered harmonic:

$$[\mathbf{Z}(\Omega)] = [[\nabla]]^2 \otimes [\mathbf{M}] + [\nabla] \otimes [\mathbf{C}] + [\mathbf{I}_{2H+1}] \otimes [\mathbf{K}] \quad (10)$$

### 3.2. Numerical strategy

The nonlinear algebraic system (9) has to be solved in an iterative way. The expression of the nonlinear forces and their derivatives with respect to the unknowns (displacement and frequency in case of a nonlinear forced response analysis) are thus required at each iteration of the solver (detailed hereafter). Two possibilities are available at this stage. The first option is to derive the expression of the Fourier coefficients analytically [28, 42], although it can prove tedious for strong, non-polynomial nonlinearities. A more straightforward approach is to use the alternating frequency/time (AFT) procedure [7]. The main underlying idea is to take advantage of the existence of a closed-form expression of the nonlinear forces in the time-domain. To this end, the displacements in the time-domain are computed by an IDFT (inverse discrete Fourier transform) or IFFT (inverse fast Fourier transform) of the Fourier coefficients stored in vector $\tilde{\mathbf{q}}$. The resulting time-domain signals are used to compute the nonlinear force $f_{nl}$ over one period and are then transferred back into the frequency domain by a DFT or FFT as illustrated in Fig. 5.

A converged solution of equation (9) is an approximation of the real solution at a single frequency. To build the nonlinear response over the whole frequency range of interest and to avoid difficulties induced by a vertical tangent at turning points, the parameter $\Omega$ has to be treated as an unknown. Among the available methods, this paper employs one of the predictor-corrector type [52], namely the arclength continuation with a tangent predictor. This method relies on parameterizing the unknowns $(\tilde{\mathbf{q}}, \Omega)$ by the curvilinear abscissa $s$. An iteration is done in two steps:
<figure><img src="image_8.png" /><figcaption>Figure 5: Illustration of the AFT procedure.</figcaption></figure>

**Prediction:** an estimation of the solution is first computed in the direction defined by the tangent of the solution branch at the previous iteration

$$ \left( \begin{array}{c} \tilde{q} \\ \Omega \end{array} \right)_{(p)} = \left( \begin{array}{c} \tilde{q} \\ \Omega \end{array} \right)_{(k)} + \Delta s \left( \begin{array}{c} \Delta \tilde{q} \\ \Delta \Omega \end{array} \right) \qquad (11) $$

where the subscripts (p) and (k) denote the predicted solution and the iteration number, respectively. The tangent vector is obtained by solving the following linear algebraic system

$$ \left[ \begin{array}{cc} \left[ \begin{array}{c} \partial_{\tilde{q}} \mathbf{R} \\ \Delta \tilde{q}^T \end{array} \right] & \partial_{\Omega} \mathbf{R} \\ \Delta \tilde{q}^T & \Delta \Omega \end{array} \right]_{(k)} \left( \begin{array}{c} \Delta \tilde{q} \\ \Delta \Omega \end{array} \right) = \left( \begin{array}{c} \mathbf{0} \\ 1 \end{array} \right) \qquad (12) $$

where the last row of the system constrains the length of the tangent vector.

**Correction:** a constraint equation is then added to the system to account for the additional unknown $\Omega$. This equation is typically the equation of a hypersphere centered on the previous converged solution point:

$$ P(\tilde{q}, \Omega, s) = (\Delta \tilde{q})^T (\Delta \tilde{q}) + \Delta \Omega^2 - \Delta s^2 = 0 \qquad (13) $$

Geometrically speaking, adding this constraint equation is tantamount to searching for the intersection between the solution curve and a hypersphere of radius $\Delta s$. The complete system of equations which forms an extended residual $\mathbf{R}_{ex} = (\mathbf{R}(\tilde{q}, \Omega) P(\tilde{q}, \Omega, s))^T$ then has to be solved with a Newton-like procedure. Linearizing the extended residual $\mathbf{R}_{ex}$ yields

$$ \left( \begin{array}{c} \tilde{q} \\ \Omega \end{array} \right)_{(k+1)} = \left( \begin{array}{c} \tilde{q} \\ \Omega \end{array} \right)_{(k)} - \left[ \begin{array}{cc} \partial_{\tilde{q}} \mathbf{R} & \partial_{\Omega} \mathbf{R} \\ \partial_{\tilde{q}} P & \partial_{\Omega} P \end{array} \right]_{(k)}^{-1} \mathbf{R}_{ex}^{(k)} \qquad (14) $$

In this study we used the MATLAB `fsolve` subroutine which is based on a trust-region algorithm [9].

### 3.3. Stability analysis

Since the harmonic balance method does not simulate the transient response of the system, a stability analysis must be conducted to assess whether the computed solution is stable or unstable. Stability can be determined by the well-known Floquet theory in the time domain [40]. However, an attractive alternative when using the harmonic balance is to use Hill's method which is a frequency domain variant of Floquet theory. The main idea is to introduce a periodic solution $\mathbf{q}$ of eq. (5) perturbed by an exponential term:

$$ \mathbf{x}_0(t) = \mathbf{q}(t) + \mathbf{p}(t)e^{\lambda t} \qquad (15) $$

Substituting this expression into the equation of motion, applying the harmonic balance formalism and a first-order Taylor expansion of the nonlinear term yields a quadratic eigenvalue problem of size $n(2H + 1)$:

$$ (\Lambda^2[\tilde{M}] + \Lambda[\tilde{C}] + [\partial_{\tilde{q}} \mathbf{R}]) \tilde{p} = 0 \qquad (16) $$
with

$$ 
\begin{aligned}
[\tilde{\mathbf{M}}] &= [\mathbf{I}_{2H+1}] \otimes [\mathbf{M}] \\
[\tilde{\mathbf{C}}] &= [\mathbf{V}] \otimes 2[\mathbf{M}] + [I_{2H+1}] \otimes [\mathbf{C}]
\end{aligned}
\qquad (17) $$

Equation (16) can be linearized [54, 10] to avoid the expensive computation of a quadratic eigenvalue problem at the cost of doubling the problem size. There exist several linearized forms. The first companion form was adopted in this study:

$$ \begin{bmatrix} [\tilde{\mathbf{C}}] & [\partial_{\tilde{\mathbf{q}}} \mathbf{R}] \\ -[\mathbf{I}_{n(2H+1)}] & [\mathbf{0}] \end{bmatrix} + \lambda \begin{bmatrix} [\tilde{\mathbf{M}}] & [\mathbf{0}] \\ [\mathbf{0}] & [\mathbf{I}_{n(2H+1)}] \end{bmatrix} = 0 \qquad (18) $$

Solving this eigenvalue problem yields $2n(2H + 1)$ solutions which exceed the number of theoretically expected eigenvalues ($2n$) by a factor of $(2H + 1)$. It is thus required to discriminate the redundant eigenvalues, stemming from the presence of multiple harmonics, from the ones with a physical signification. There exist two methods for identifying the relevant eigenvalues:

* sorting the eigenvalues consists in only keeping the $2n$ eigenvalues with smallest imaginary parts. This method has been proved to converge as the number of retained harmonics is increased [34] and is used in this study.

* sorting the eigenvalues consists in keeping the $2n$ eigenvalues associated with the most symmetric eigenvectors [19, 31]. This approach appears to be more robust, although there is no formal mathematical proof regarding its convergence.

A converged solution point is eventually deemed stable if the real part of all retained eigenvalues is negative, and unstable otherwise.

### 3.4. Detection of bifurcations

A bifurcation corresponds to a topological change of the vector field describing the dynamic behaviour of the system. They are usually detected by dedicated test functions, whose zeros correspond to a specific bifurcation. Branch points (BP) and saddle-node (SN) bifurcations are detected when the eigenvalue matrix $[V] = diag(\lambda)$ obtained with Hill's Method is singular, that is:

$$ \phi_{BP,SN} = \det([V]) = |[V]| = 0 \qquad (19) $$

Both bifurcations are detected with this test function. A saddle-node bifurcation is associated with a change of direction of the solution curve relative to the continuation parameter. SN bifurcations are thus discriminated from BP bifurcations by monitoring the sign of the frequency-component of the tangent vector at each iteration:

$$ \phi_{SN} = \Delta\Omega \qquad (20) $$

A third test function is defined to detect Neimark-Sacker (NS) bifurcations, corresponding to the creation of a branch of quasi-periodic solutions, and period-doubling (PD) bifurcations. Both bifurcations are characterized by a pair of complex conjugate eigenvalues crossing the imaginary axis. Their detection rely on the bialternate matrix product [16]. For $[A], [B]$ two square matrices of size $n$, the bialternate product $A \odot B$ is defined as

$$ [[A] \odot [B]]_{(p,q),(r,s)} = \frac{1}{2} \left( \begin{vmatrix} a_{pr} & a_{ps} \\ b_{qr} & b_{qs} \end{vmatrix} + \begin{vmatrix} b_{pr} & b_{ps} \\ a_{qr} & a_{qs} \end{vmatrix} \right) \qquad (21) $$

with $(1 \le q < p \le n)$ and $(1 \le s < r \le n)$. The product $[2[A] \odot [I_n]]$ is a diagonal matrix composed of the sums $\lambda_i + \lambda_j$ of the eigenvalues of matrix $[A]$ [15] and is singular when a pair of complex conjugate eigenvalues crosses the imaginary axis. A suitable test function is thus:

$$ \phi_{NS,PD} = |[2[V] \odot [I]]| \qquad (22) $$

In order to differentiate between NS and PD bifurcations, one can monitor the imaginary part of the crossing eigenvalues. In case of a PD bifurcation, the imaginary part equals $\Omega/2$. A NS bifurcation is detected otherwise. In the special case of eigenvalues that are real conjugates, a neutral saddle point is detected and can therefore be ignored.
Due to the piecewise linear nature of the nonlinearity, another type of bifurcation can arise in this vibro-impact system, a so-called discontinuity-induced bifurcation. At the frontier between impacting and non-impacting orbits there exists grazing orbits characterized by impacts with zero normal velocity. This transition where, under a change of parameter, an initially linear response loses regularity until it exhibits vibro-impacts is termed a grazing bifurcation [22]. Contrary to their smooth counterparts, non-smooth bifurcations are not necessarily associated with eigenvalues smoothly crossing the imaginary axis. It is thus impossible to employ classical bifurcation analysis tools. In the present work, grazing bifurcations are detected by monitoring the mean value of the Heaviside function of the nonlinear force (2) over a rotation period at each frequency

$$ \phi_{GRA} = \mathcal{H} \left( \left\langle \mathcal{H} (\mathbf{G}^T \mathbf{q} - g(t)) \right\rangle - 1 \right) - \mathcal{H} \left( 1 - \left\langle \mathcal{H} (\mathbf{G}^T \mathbf{q} - g(t)) \right\rangle \right) \quad (23) $$

where $\langle \cdot \rangle$ denotes the mean value over one fundamental period. Such function equals 1 when the gear teeth stay in contact and -1 in case of contact loss, ensuring zero-crossings at the onset of impacting orbits.

# 4. Results and discussion

The numerical procedure introduced in the previous section enables the analysis of the aforementioned transmission model. First, the computation of a reference configuration, namely, a transmission without holes in the gear webs, provides an understanding of the underlying dynamics. Five lightweight gear configurations are then investigated and their effect on the nonlinear response are discussed. The analysis is carried out for two different input torques: $T = 20$ N·m and $T = 115$ N·m, the latter corresponding to the optimal torque of the reference configuration for which the STE fluctuation is minimised. The effect of holes is expected to be more pronounced as the applied torque and consequently the static deflection is increased.

## 4.1. Reference configuration

A modal damping ratio of $\xi = 3\%$ is applied on all modes. The physical damping matrix is defined as

$$ [C] = [M]^T [\Psi] [\text{diag}(2\xi_j\omega_j)] [\Psi]^T [M] \quad (24) $$

where $[\Psi]$ is the modal basis (with mass normalized modes) of the whole transmission computed with the average mesh stiffness linking the two gears. No excitation sources other than the STE and mesh stiffness fluctuations are considered. The system's response is investigated between $\Omega = 10$ rad·s⁻¹ and $\Omega = 1466$ rad·s⁻¹, corresponding to a maximum rotational speed of 14000 rpm. To avoid any numerical issue, the continuation procedure is initiated in the linear regime by applying the HBM to system (5) with a dedicated linear element that enforces contact be kept at all time. It is derived by linearizing the mesh force around the static equilibrium:

$$ \mathbf{f}_{nl}(\mathbf{q}) \approx \mathbf{f}_{nl}(\mathbf{q}_s) + \left. \frac{\partial \mathbf{f}_{nl}}{\partial \mathbf{q}} \right|_{\mathbf{q}_s} (\mathbf{q} - \mathbf{q}_s) \quad (25) $$

Since only a constant input torque is applied, it comes:

$$ [M]\ddot{\mathbf{q}} + [C]\dot{\mathbf{q}} + [K](\mathbf{q} - \mathbf{q}_s) + \left. \frac{\partial \mathbf{f}_{nl}}{\partial \mathbf{q}} \right|_{\mathbf{q}_s} (\mathbf{q} - \mathbf{q}_s) = \mathbf{0} \quad (26) $$

### 4.1.1. Input data and selection of harmonics

The STE and mesh stiffness are computed with the finite element analysis introduced in section 2.2. FFTs are performed on the resulting time domain signals and 7 harmonics are used to reconstruct the STE and mesh stiffness for the nonlinear analysis (see Fig. 6). The response being approximated by means of truncated Fourier series, the harmonic order ought to be as high as possible for a satisfying approximation of the response while remaining as low as possible for a fast computation. In our case, the STE and the mesh stiffness fluctuations are characterized by harmonics of the mesh frequency $H_m = kZ_1$ with $k \in \mathbb{N}$ and $Z = Z_1 = Z_2 = 50$. In addition, the convolution between mesh stiffness fluctuations and the STE is expected to generate high frequency harmonics. In addition, the nonlinearity is responsible for the creation of harmonic orders that are different from the mesh harmonics $H_m = kZ_1$, $k \in \mathbb{N}$.
<figure><img src="image_9.png" /><figcaption>Figure 6: Static transmission error (a) and mesh stiffness (b) computed at $T = 20$ N·m and Static transmission error (c) and mesh stiffness (d) computed at $T = 115$ N·m. Results from the FE computation (---) and signals reconstructed with 7 harmonics (—)</figcaption></figure>

With $Z_1 = 50$ teeth, keeping 7 harmonics of the internal excitation results in an excitation up to harmonic $H_{350}$ of the rotation. Experience shows that keeping harmonics $kZ_1, k \in \mathbb{N}$, up to 2 times the maximum excitation order ensures sufficient accuracy for the computation of the primary resonance and a good representation of the superharmonic resonances. Additional harmonics around each meshing harmonics are included for the convergence of the detection of grazing bifurcations although their contribution to the displacement remains small.

### 4.1.2. Results for an input torque $T = 20$ N·m

The following results correspond to the nonlinear forced response analysis of the reference configuration. Figure 7 depicts the standard deviation (SD) of the dynamic transmission error $DTE(t)$, defined as

$$DTE(t) = G^T q(t) \qquad (27)$$

The standard deviation, denoted $DTE_{SD}$ in the following, is plotted versus the angular frequency of the input shaft.

There are 10 bifurcations in the frequency range of interest. The primary resonance, near $\Omega = 650$ rad·s⁻¹, exhibits two saddle-node (SN) bifurcations responsible for the well-known amplitude-jump phenomenon. A small interval of stable vibro-impact regimes is observed as the grazing bifurcation occurs before the loss of stability due to the SN bifurcation. This happens because we consider damped dynamics. In the conservative case, the transition between impacting and non-impacting orbits is expected to take place through a SN bifurcation [35, 36].

At lower frequencies, two resonances arise at around $\Omega = 230$ rad·s⁻¹ and $\Omega = 370$ rad·s⁻¹, respectively. The latter leads to the creation of a stable vibro-impact region that spans $\Delta\Omega = 8$ rad·s⁻¹, while the former exhibits another unstable vibro-impact region.

The response near the low frequency resonances (see Fig. 9) is characterized by a harmonic content with higher order harmonics than that of the behaviour around the primary resonance (see Fig. 8), which is excited only by the first harmonic of the mesh frequency $H_{50}$. Indeed, both resonances are excited by several harmonics of the internal excitation, as evidenced by the Campbell diagram displayed in Fig. A.1. This leads to shorter free-flight periods as highlighted by the time series and phase portraits shown in Fig. 9. One can note that, despite the shorter free-flight events, the impact velocity remains of the same order of magnitude.
<figure><img src="image_10.png" /><figcaption>**Figure 7:** Forced response curve of the reference configuration for an input torque $T = 20 \text{ N} \cdot \text{m}$ (a) and effect of harmonic truncation order (b). Dashed lines (---) indicate unstable regions and solid lines (—) indicate stable regions. Saddle-node and grazing bifurcations are represented with red circle (●) and magenta diamond (♦) markers, respectively.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 8:** Dynamic transmission error and phase portraits between the grazing and saddle-node bifurcations (a), (c) and beyond (b), (d) the second saddle-node bifurcation on the primary resonance</figcaption></figure>

### 4.1.3. Results for an input torque $T = 115 \text{ N} \cdot \text{m}$

Results for an input torque $T = 115 \text{ N} \cdot \text{m}$ (displayed in Fig. 10) show a much more linear behaviour of the gear transmission over the whole frequency range. Indeed, the only vibro-impact region lies on the primary resonance, around $\Omega = 800 \text{ rad} \cdot \text{s}^{-1}$. More importantly, the vibro-impact region only spans $\Delta\Omega = 30 \text{ rad} \cdot \text{s}^{-1}$ instead of more than $\Delta\Omega = 700 \text{ rad} \cdot \text{s}^{-1}$ with an input torque of 20 N·m.

One can see that the primary resonance peak remains linear until the vibration amplitude becomes sufficiently large
<figure><img src="image_12.png" /><figcaption>**Figure 9:** Dynamic transmission error and phase portraits before (a), (c) and past (b), (d) the first grazing bifurcation near $\Omega = 230$ rad·s$^{-1}$</figcaption></figure>

(about 5.2 µm in SD value) for contact to be lost. This is expected as the increased static load causes a larger static deflection of the gear teeth (mean value of the STE of 11.2 µm at $T = 115$ N·m compared to a mean value of 2.9 µm at $T = 20$ N·m). Furthermore, tooth profile modifications (see Table 2) are such that the internal excitation is minimum at $T = 115$ N·m. As a result, the system only exhibits a mild softening behaviour. At low frequencies, we observe a number of resonances generated by the excitation of modes by high order harmonics of the mesh frequency.

<figure><img src="image_13.png" /><figcaption>**Figure 10:** Forced response curve of the reference configuration for an input torque $T = 115$ N·m. Dashed lines (- -) indicate unstable regions and solid lines (—) indicate stable regions. Saddle-node and grazing bifurcations are represented with red circle (●) and magenta diamond (♦) markers, respectively.</figcaption></figure>
## 4.2. Parametric study of the influence of gear blank topology

In the following, the influence of gear blank topology on the dynamic response is investigated with several configurations. The hypothesis of constant mass is made so that each configuration is defined, and referred to, by the number of holes and their radial position in the gear blank (see Fig. 11). For instance, a configuration with 8 holes at a radial position of $R = 30$ mm from the centre of the gear is labelled 8R30.

<figure><img src="image_14.png" /><figcaption>**Figure 11:** Schematic representation of the holed configurations</figcaption></figure>

The STE and mesh stiffness fluctuations are computed by means of the novel numerical method introduced in [5] and presented in section 2. The effect of holes on the static transmission error and mesh stiffness are discussed in citeYouness. For the sake of completeness, we briefly recall the main conclusions.

The addition of holes in the gear blanks is responsible for the appearance of a low frequency component in the STE generated by the passage of holes at the angular position of the mesh zone. The created harmonic order thus corresponds to the number of holes $N_h$ in the analysed configuration. This also impacts the mesh stiffness, as sidebands $H_{kZ\pm l/N_h}$, $(k,l) \in \mathbb{N}^2$ are created around each harmonic of the mesh frequency, resulting in a modulation visible in Fig. 12.

<figure><img src="image_15.png" /><figcaption>**Figure 12:** Exemple of static transmission error (a) and mesh stiffness (b) computed at $T = 20$ N·m for a gear configuration with 8 holes. Static transmission error (c) and mesh stiffness (d) computed at $T = 115$ N·m. Results from the FE computation (---) and signals reconstructed with 7 harmonics (—)</figcaption></figure>
The time series of the input data used in this study are shown in Appendix C and D. The influence of the radial position of the holes can be assessed by looking at the STE and mesh stiffness of design cases 8R27, 8R30 and 8R34. Figures C.1-(c),(d),(e),(f),(g),(h) make it apparent that the harmonic order $H_{N_h}$ is amplified as radius R is increased, i.e, when the holes are placed close to the gear teeth. The effect of the number of holes can be apprehended by looking at Fig. D.1-(g),(h),(i),(j). It clearly appears that the higher the number of holes, the smaller the STE fluctuations are. Note that the same trends can be observed on the time series of the mesh stiffness.

### 4.2.1. Effect of holes at $T = 20$ N·m

As with the reference computation, we are interested in the evolution of $DTE_{SD}$ with the rotational speed. The results for all tested configurations are summarised in Fig. 13 and comparisons between the reference and individual configurations are reported in Fig. 14 for more clarity.

<figure><img src="image_16.png" /><figcaption>**Figure 13:** Forced response of the transmission model to the internal excitation (a) with close-up on the primary resonance (b) and main superharmonic resonance (c). Dashed lines indicate unstable regions and solid lines indicate stable regions. Saddle-node, Neimark-Sacker and grazing bifurcations are represented with circle (●), triangle (▲) and diamond (♦) markers, respectively. Each color (see online version) corresponds to a test case</figcaption></figure>

These figures clearly show that the addition of holes in the gear blanks give rise to a more complex and rich bifurcation structure along the main solution branch. Indeed, Neimark-Sacker bifurcations are detected for all holed configurations. This indicates that a quasiperiodic (QP) attractor is created in the vicinity of the frequencies at which the NS bifurcations happen.

This should motivate future research, as the 1-T periodic solution loses its stability in favour of a QP solution (with two or more incommensurate frequencies) which cannot be investigated using classical harmonic balance algorithms. A few workarounds were proposed in the literature, e.g., by implementing a multidimensional harmonic balance [21, 51] or considering the greatest common divisor of the incommensurate frequencies as an approximate (monodimensional)
<figure><img src="image_17.png" /><figcaption>**Figure 14:** Nonlinear forced response curves of all tested configurations against the reference for a static torque *T* = 20 N·m: (a) 6R30, (b) 8R27, (c) 8R30, (d) 8R34, (e) 10R34. Dashed lines indicate unstable regions and solid lines indicate stable regions. Saddle-node, Neimark-Sacker and grazing bifurcations are represented with circle (●), triangle (▲) and diamond (♦) markers, respectively.</figcaption></figure>

fundamental frequency [20].

All configurations but 8R27 exhibit two NS bifurcations on the upper branch of the primary resonance. The 8R27 configuration does not show the smaller NS bifurcation region ranging from about $\Omega$ = 810 rad·s⁻¹ to $\Omega$ = 900 rad·s⁻¹. Note that the term *smaller* refers to the frequency interval between two NS bifurcations on the main curve as it is impossible to presuppose the extent of the QP branch without computing it. The suppression of this NS region results from the decrease of the frequency of the last grazing bifurcation to the frequency interval of the second region of NS bifurcations. Interestingly, the extent of the frequency range affected by these bifurcations seems robust to changes in gear blank topology.
Regarding the superharmonic resonances, one can see that all holed designs successfully cancel out the vibro-impact region arising between $\Omega = 368\,\text{rad}\cdot\text{s}^{-1}$ and $\Omega = 375\,\text{rad}\cdot\text{s}^{-1}$ in the reference configuration. Furthermore, on the main superharmonic resonance, all tested configurations suppress the SN bifurcation although they are not able to prevent contact loss.

One of the more interesting results of this study is the variability of the range of vibro-impact regions located around the primary resonance. Due to the softening nonlinearity, the vibro-impact region around the primary resonance is bounded by the SN bifurcation located at the top of the branch and a grazing bifurcation. Figure 15 summarises the location of detected grazing bifurcations and ensuing contact loss events. In the reference configuration, the vibro-impact region spans $\Delta\Omega = 890\,\text{rad}\cdot\text{s}^{-1}$ while holed configurations exhibit contact loss in frequency intervals whose length ranges from $\Delta\Omega = 769\,\text{rad}\cdot\text{s}^{-1}$ to $\Delta\Omega = 303\,\text{rad}\cdot\text{s}^{-1}$.

<figure><img src="image_18.png" /><figcaption>**Figure 15:** Bifurcation map for a static torque $T = 20\,\text{N}\cdot\text{m}$. Grazing bifurcations and the SN bifurcations located at the top of the primary resonances are indicated by colored diamond (◇) and circle (○) markers whose colors correspond to the tested configuration. Vibro-impact regions are represented by double-sided arrows.</figcaption></figure>

The design parameter governing this phenomenon is not easily identified. Results suggest that placing holes closer to the shaft lead to a shrinkage of the frequency interval exhibiting contact loss (see configurations 8R34, 8R30 and 8R27). This is associated to the increase of the mean value of the STE as holes are placed closer to the shaft, hence increasing the static deflection and the distance to cover before losing contact. Besides, for a fixed radial position, larger holes do not necessarily lead to a narrower vibro-impact region. Indeed, switching from configuration 8R30 to 6R30 does provide an enhancement, however configuration 10R34 shows better results than configuration 8R34. This could indicate the presence of a saddle point in the parameter space.

### 4.2.2. Effect of holes at $T = 115\,\text{N}\cdot\text{m}$

As the load is increased, so is the hole-induced deflection. This results in a strong variation of the mean value of the mesh stiffness and a significant frequency shift of the primary resonance of the system. This shift is clearly observed in figure 16 where the primary resonance of the reference lies at $\Omega = 800\,\text{rad}\cdot\text{s}^{-1}$ and shifts to a little above $\Omega = 600\,\text{rad}\cdot\text{s}^{-1}$ for the most flexible case 8R27. As expected with piecewise-smooth systems, the system behaviour is linear until contact is lost, which does not happen with all holed configurations. This is a particularly interesting result, since this completely suppresses the vibro-impact-induced hammering noise. The linear response can be a source of whining noise. However, classical strategies such as adding tooth profile modifications [14] to reduce the peak-to-peak STE can then be employed for mitigating the resulting noise.

## 5. Conclusion

Gears are a central part of many mechanical systems and engineering applications. They are often designed with holes in order to meet the ever-increasing demand for lightweight structures. It is of paramount importance to have a
<figure><img src="image_19.png" /><figcaption>**Figure 16:** Forced response curve of the reference configuration for a static torque *T* = 115 N·m. Dashed lines (---) indicate unstable regions and solid lines (—) indicate stable regions. Saddle-node and grazing bifurcations are represented with circle (●) and diamond (♦) markers, respectively.</figcaption></figure>

profound understanding of the response of nonlinear systems as they behave quite differently from their linear counterparts. With that in mind, this work presents a parametric study of lightweight gear designs in terms of nonlinear forced response and bifurcation analysis of periodic orbits.

The study shows that adding holes in the gear blank has a real influence for both low and high input torques. For high input torques, the response of the reference configuration was found to be mostly linear, except around the primary resonance, where a small vibro-impact region and saddle-node bifurcation was detected. All lightweight configurations lead to the cancellation of the contact loss events. These findings are of particular interest since most studies and research towards vibrations reduction of geared systems aim at reducing the peak-to-peak fluctuations of the STE. This is a perfectly sound endeavour, provided the transmission stays in a linear regime. However, this study demonstrates this could prove insufficient as the input torque *T* = 115 N·m is the optimal torque of the reference configuration, i.e, the torque at which the STE fluctuations are minimised, and contact loss still occured. Moreover, although the investigated designs are not able to prevent vibro-impacts from happening at low torques, they clearly have a beneficial effect in curtailing the frequency range affected by intermittent contact. These results thus offer new insights and research perspectives when dealing with nonlinear gear dynamics.

**Acknowledgements**

This work was performed within the framework of the LabCom LADAGE (LAboratoire de Dynamique des en-grenAGEs), created by the LTDS and the Vibratec Company and operated by the French National Research Agency (ANR-14-LAB6-0003). It was also performed within the framework of the LABEX CeLyA (ANR-10-LABX-0060) of Université de Lyon, within the program "Investissements d'Avenir" (ANR-16-IDEX-0005) operated by the French National Research Agency (ANR). Author A. Mélot acknowledges Camille Grégoire from the LTDS for the precious feedback regarding the manuscript.

**Conflict of interest**

The authors declare that they have no conflict of interest.

**CRediT authorship contribution statement**

**Adrien Mélot:** Conceptualization, Methodology, Software, Validation, Investigation, Data Curation, Writing - Original Draft, Writing - Review and Editing, Visualization. **Youness Benaïcha:** Conceptualization, Validation, Resources, Data Curation, Writing - Original Draft, Writing - Review and Editing. **Emmanuel Rigaud:** Supervision, Funding acquisition, Writing - Review and Editing. **Joël Perret-Liaudet:** Supervision, Funding acquisition, Writing - Review and Editing. **Fabrice Thouverez:** Funding acquisition, Writing - Review and Editing.
A. Campbell diagram of the reference configuration

<figure><img src="image_20.png" /><figcaption>**Figure A.1:** Campbell diagram for the reference configuration showing the frequency coincidences between the linear modes of the system and the harmonics of the internal excitation.</figcaption></figure>

B. Harmonic spectra of the internal excitation

<figure><img src="image_21.png" /><figcaption>**Figure B.1:** Harmonic content of the fluctuations of the static transmission error (a) and mesh stiffness (b) of the reference configuration at *T* = 115 N·m and fluctuations of the static transmission error (c) and mesh stiffness (d) of a 10-hole configuration computed at *T* = 115 N·m.</figcaption></figure>
### C. Input data at $T = 20$ N·m

<figure><img src="image_22.png" /><figcaption>**Figure C.1:** Static transmission error and mesh stiffness of the holed gear configurations computed at $T = 20$ N·m. 6R30 (a,b), 8R27 (c,d), 8R30 (e,f), 8R34 (g,h), 10R34 (i,j)</figcaption></figure>
D. Input data at $T = 115$ N·m

<figure><img src="image_23.png" /><figcaption>Figure D.1: Static transmission error and mesh stiffness of the holed gear configurations computed at $T = 115$ N·m. 6R30 (a,b), 8R27 (c,d), 8R30 (e,f), 8R34 (g,h), 10R34 (i,j)</figcaption></figure>
References

[1] Al-shyyab, A., Kahraman, A., 2005a. Non-linear dynamic analysis of a multi-mesh gear train using multi-term harmonic balance method: period-one motions. Journal of Sound and Vibration 284, 151–172. doi:10.1016/j.jsv.2004.06.010.

[2] Al-shyyab, A., Kahraman, A., 2005b. Non-linear dynamic analysis of a multi-mesh gear train using multi-term harmonic balance method: sub-harmonic motions. Journal of Sound and Vibration 279, 417–451. doi:10.1016/j.jsv.2003.11.029.

[3] Alcorta, R., Baguet, S., Prabel, B., Piteau, P., Jacquet-Richardet, G., 2019. Period doubling bifurcation analysis and isolated sub-harmonic resonances in an oscillator with asymmetric clearances. Nonlinear Dynamics 98, 2939–2960. doi:10.1007/s11071-019-05245-6.

[4] Barthod, M., Hayne, B., Tébec, J.L., Pin, J.C., 2007. Experimental study of dynamic and noise produced by a gearing excited by a multi-harmonic excitation. Applied Acoustics 68, 982–1002. doi:https://doi.org/10.1016/j.apacoust.2006.04.012.

[5] Benaicha, Y., Mélot, A., Rigaud, E., Beley, J.D., Thouverez, F., Perret-Liaudet, J., 2021. A decomposition method for the fast computation of the transmission error of gears with holes. Working paper or preprint.

[6] Benaicha, Y., Perret-Liaudet, J., Beley, J.D., Rigaud, E., Thouverez, F., 2021. On a flexible multibody modelling approach using FE-based contact formulation for describing gear transmission error. Working paper or preprint.

[7] Cameron, T.M., Griffin, J., 1989. An Alternating Frequency/Time Domain Method for Calculating the Steady-State Response of Nonlinear Dynamic Systems. Journal of Applied Mechanics. doi:10.1115/1.3176036.

[8] Carbonelli, A., Rigaud, E., Perret-Liaudet, J., 2016. Vibro-Acoustic Analysis of Geared Systems—Predicting and Controlling the Whining Noise. Springer International Publishing, Cham. pp. 63–79. doi:10.1007/978-3-319-24055-8_5.

[9] Conn, A.R., Gould, N.I.M., Toint, P.L., 2000. Trust Region Methods. Society for Industrial and Applied Mathematics. doi:10.1137/1.9780898719857.

[10] Detroux, T., Renson, L., Masset, L., Kerschen, G., 2015. The harmonic balance method for bifurcation analysis of large-scale nonlinear mechanical systems. Computer Methods in Applied Mechanics and Engineering 296, 18–38. doi:10.1016/j.cma.2015.07.017.

[11] Dion, J.L., Le Moyne, S., Chevallier, G., Sebbah, H., 2009. Gear impacts and idle gear noise: Experimental study and non-linear dynamic model. Mechanical Systems and Signal Processing 23, 2608–2628. doi:https://doi.org/10.1016/j.ymssp.2009.05.007.

[12] Fontanela, F., Grolet, A., Salles, L., Hoffmann, N., 2019. Computation of quasi-periodic localised vibrations in nonlinear cyclic and symmetric structures using harmonic balance methods. Journal of Sound and Vibration 438, 54–65. doi:https://doi.org/10.1016/j.jsv.2018.09.002.

[13] Garambois, P., Donnard, G., Rigaud, E., Perret-Liaudet, J., 2017a. Multiphysics coupling between periodic gear mesh excitation and input/output fluctuating torques: Application to a roots vacuum pump. Journal of Sound and Vibration 405, 158–174. doi:10.1016/j.jsv.2017.05.043.

[14] Garambois, P., Perret-Liaudet, J., Rigaud, E., 2017b. NVH robust optimization of gear macro and microgeometries using an efficient tooth contact model. Mechanism and Machine Theory 117, 78–95. doi:https://doi.org/10.1016/j.mechmachtheory.2017.07.008.

[15] Govaerts, W., 2000. Numerical bifurcation analysis for odes. Journal of Computational and Applied Mathematics 125, 57–68. doi:https://doi.org/10.1016/S0377-0427(00)00458-1.

[16] Govaerts, W., Sijnave, B., 1999. Matrix manifolds and the jordan structure of the bialternate matrix product. Linear Algebra and its Applications 292, 245–266. doi:https://doi.org/10.1016/S0024-3795(99)00039-7.

[17] Guilbert, B., Philippe, C., Velex, P., 2017. Static and dynamic analyses of thin-rimmed gears with holes, in: VDI, International Conference on Gears, Munich, Germany.

[18] Guilbert, B., Velex, P., Dureisseix, D., Cutuli, P., 2016. A Mortar-Based Mesh Interface for Hybrid Finite-Element/Lumped-Parameter Gear Dynamic Models—Applications to Thin-Rimmed Geared Systems. Journal of Mechanical Design 138. doi:10.1115/1.4034220.

[19] Guillot, L., Lazarus, A., Thomas, O., Vergez, C., Cochelin, B., 2020. A purely frequency based floquet-hill formulation for the efficient stability computation of periodic solutions of ordinary differential systems. Journal of Computational Physics 416, 109477. doi:https://doi.org/10.1016/j.jcp.2020.109477.

[20] Guskov, M., Thouverez, F., 2012. Harmonic Balance-Based Approach for Quasi-Periodic Motions and Stability Analysis. Journal of Vibration and Acoustics 134. doi:10.1115/1.4005823.

[21] Haslam, A.H., Schwingshackl, C.W., Rix, A.I.J., 2020. A parametric study of an unbalanced Jeffcott rotor supported by a rolling-element bearing. Nonlinear Dynamics 99, 2571–2604. doi:10.1007/s11071-020-05470-4.

[22] Ibrahim, R.A., 2009. Vibro-Impact Dynamics. Springer-Verlag Berlin Heidelberg. doi:https://doi.org/10.1007/978-3-642-00275-5.

[23] Jaumouillé, V., Sinou, J.J., Petitjean, B., 2010. An adaptive harmonic balance method for predicting the nonlinear dynamic responses of mechanical systems—application to bolted structures. Journal of Sound and Vibration 329, 4048–4067. doi:https://doi.org/10.1016/j.jsv.2010.04.008.

[24] Kadmiri, Y., Perret-Liaudet, J., Rigaud, E., Le Bot, A., Vary, L., 2011. Influence of Multiharmonics Excitation on Rattle Noise in Automotive Gearboxes. Advances in Acoustics and Vibration 2011, 659797. doi:10.1155/2011/659797. publisher: Hindawi Publishing Corporation.

[25] Kadmiri, Y., Rigaud, E., Perret-Liaudet, J., Vary, L., 2012. Experimental and numerical analysis of automotive gearbox rattle noise. Journal of Sound and Vibration 331, 3144–3157. doi:https://doi.org/10.1016/j.jsv.2012.02.009.

[26] Kahraman, A., Singh, R., 1990. Non-linear dynamics of a spur gear pair. Journal of Sound and Vibration 142, 49–75. doi:10.1016/0022-460X(90)90582-K.

[27] Karagiannis, K., Pfeiffer, F., 1991. Theoretical and experimental investigations of gear-rattling. Nonlinear Dynamics 2, 367–387. doi:10.1007/BF00045670.

[28] Krack, M., Panning-von Scheidt, L., Wallaschek, J., 2013. A high-order harmonic balance method for systems with distinct states. Journal of Sound and Vibration 332, 5476–5488. doi:10.1016/j.jsv.2013.04.048.

[29] Kunert, A., Pfeiffer, F., 1990. Stochastic model for rattling in gear-boxes, in: Schiehlen, W. (Ed.), Nonlinear Dynamics in Engineering Systems, Springer Berlin Heidelberg, Berlin, Heidelberg. pp. 173–180.
[30] Laursen, T.A., Maker, B.N., 1995. An augmented lagrangian quasi-newton solver for constrained nonlinear finite element applications. International Journal for Numerical Methods in Engineering 38, 3571-3590. doi:https://doi.org/10.1002/nme.1620382103.

[31] Lazarus, A., Thomas, O., 2010. A harmonic-based method for computing the stability of periodic solutions of dynamical systems. Comptes Rendus Mécanique 338, 510-517. doi:https://doi.org/10.1016/j.cme.2010.07.020.

[32] Liu, C., Qin, D., Wei, J., Liao, Y., 2018. Investigation of nonlinear characteristics of the motor-gear transmission system by trajectory-based stability preserving dimension reduction methodology. Nonlinear Dynamics 94, 1835-1850. URL: https://doi.org/10.1007/s11071-018-4460-2, doi:10.1007/s11071-018-4460-2.

[33] Liu, G., Parker, R.G., 2008. Nonlinear dynamics of idler gear systems. Nonlinear Dynamics 53, 345-367. doi:10.1007/s11071-007-9317-z.

[34] Moore, G., 2005. Floquet theory as a computational tool. SIAM Journal on Numerical Analysis 42, 2522-2568.

[35] Mora, K., Budd, C., Glendinning, P., Keogh, P., 2014. Non-smooth hopf-type bifurcations arising from impact-friction contact events in rotating machinery. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 470, 20140490. doi:10.1098/rspa.2014.0490.

[36] Mora, K., Champneys, A.R., Shaw, A.D., Friswell, M.I., 2020. Explanation of the onset of bouncing cycles in isotropic rotor dynamics; a grazing bifurcation analysis. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 476, 20190549. doi:10.1098/rspa.2019.0549.

[37] Neufond, J., Denimal, E., Rigaud, E., Perret-Liaudet, J., Carbonelli, A., 2019. Whining noise computation of a planetary gear set induced by the multi-mesh excitations. Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science 233, 7236-7245.

[38] Nevzat Özgüven, H., Houser, D., 1988. Mathematical models used in gear dynamics—a review. Journal of Sound and Vibration 121, 383-411. doi:https://doi.org/10.1016/S0022-460X(88)80365-1.

[39] Ottewill, J., Neild, S., Wilson, R., 2010. An investigation into the effect of tooth profile errors on gear rattle. Journal of Sound and Vibration 329, 3495-3506. doi:https://doi.org/10.1016/j.jsv.2010.03.014.

[40] Peletan, L., Baguet, S., Torkhani, M., Jacquet-Richardet, G., 2013. A comparison of stability computational methods for periodic solution of nonlinear problems with application to rotordynamics. Nonlinear Dynamics 72, 671-682. doi:10.1007/s11071-012-0744-0.

[41] Perret-Liaudet, J., Rigaud, E., 2007. Some Effects of Gear Eccentricities on Automotive Rattle Noise, in: IDETC-CIE2007, Volume 7: 10th International Power Transmission and Gearing Conference. pp. 561-568. doi:10.1115/DETC2007-34794.

[42] Petrov, E.P., Ewins, D.J., 2003. Analytical Formulation of Friction Interface Elements for Analysis of Nonlinear Multi-Harmonic Vibrations of Bladed Disks. Journal of Turbomachinery 125, 364-371. doi:10.1115/1.1539868.

[43] Pfeiffer, F., Prestl, W., 1994. Hammering in diesel-engine driveline systems. Nonlinear Dynamics 5, 477-492. doi:10.1007/BF00052455.

[44] Quaegbeur, S., Chouvion, B., Thouverez, F., Berthe, L., 2020. Energy transfer between nodal diameters of cyclic symmetric structures exhibiting polynomial nonlinearities: Cyclic condition and analysis. Mechanical Systems and Signal Processing 139, 106604. doi:https://doi.org/10.1016/j.ymssp.2019.106604.

[45] Raghothama, A., Narayanan, S., 1999. Bifurcation and chaos in geared rotor bearing system by incremental harmonic balance method. Journal of Sound and Vibration 226, 469-492. doi:https://doi.org/10.1006/jsvi.1999.2264.

[46] Rigaud, E., Perret-Liaudet, J., 2020. Investigation of gear rattle noise including visualization of vibro-impact regimes. Journal of Sound and Vibration 467, 115026. doi:https://doi.org/10.1016/j.jsv.2019.115026.

[47] Rigaud, E., Sabot, J., 1996. Effect of Elasticity of Shafts, Bearings, Casing and Couplings on the Critical Rotational Speeds of a Gearbox, in: International Conference on Gears, Dresde, Germany. pp. 833-845.

[48] Rocca, E., Russo, R., 2011. Theoretical and experimental investigation into the influence of the periodic backlash fluctuations on the gear rattle. Journal of Sound and Vibration 330, 4738-4752. doi:https://doi.org/10.1016/j.jsv.2011.04.008.

[49] Salles, L., Staples, B., Hoffmann, N., Schwingshackl, C., 2016. Continuation techniques for analysis of whole aeroengine dynamics with imperfect bifurcations and isolated solutions. Nonlinear Dynamics 86, 1897-1911. doi:10.1007/s11071-016-3003-y.

[50] Sato, K., Yamamoto, S., kawakami, T., 1991. Bifurcation sets and chaotic states of a gear system subjected to harmonic excitation. Computational Mechanics 7, 173-182. doi:10.1007/BF00369977.

[51] Schilder, F., Vogt, W., Schreiber, S., Osinga, H.M., 2006. Fourier methods for quasi-periodic oscillations. International Journal for Numerical Methods in Engineering 67, 629-671. doi:10.1002/nme.1632.

[52] Seydel, R., 2010. Practical Bifurcation and Stability Analysis. Springer-Verlag New York. doi:10.1007/978-1-4419-1740-9.

[53] Shweiki, S., Palermo, A., Mundo, D., 2017. A Study on the Dynamic Behaviour of Lightweight Gears. Shock and Vibration, 7982170 doi:10.1155/2017/7982170.

[54] Tisseur, F., Meerbergen, K., 2001. The quadratic eigenvalue problem. SIAM Review 43, 235-286.

[55] Urabe, M., 1965. Galerkin's procedure for nonlinear periodic systems. Archive for Rational Mechanics and Analysis 20, 120-152. URL: https://doi.org/10.1007/BF00284614, doi:10.1007/BF00284614.

[56] Wang, J., Zhang, J., Yao, Z., Yang, X., Sun, R., Zhao, Y., 2019. Nonlinear characteristics of a multi-degree-of-freedom spur gear system with bending-torsional coupling vibration. Mechanical Systems and Signal Processing 121, 810-827. URL: https://www.sciencedirect.com/science/article/pii/S0888327018307799, doi:https://doi.org/10.1016/j.ymssp.2018.12.002.

[57] Wang, M., Zhao, W., Manoj, R., 2002. Numerical modelling and analysis of automotive transmission rattle. Journal of Vibration and Control 8, 921-943. doi:https://doi.org/10.1177/10775402029594.

[58] Wei, S., Han, Q.K., Dong, X.J., Peng, Z.K., Chu, F.L., 2017. Dynamic response of a single-mesh gear system with periodic mesh stiffness and backlash nonlinearity under uncertainty. Nonlinear Dynamics 89, 49-60. doi:10.1007/s11071-017-3435-z.

[59] Welbourn, D., 1979. Fundamental knowledge of gear noise: A survey, in: Proceedings of conf. on Noise and Vibrations of Engines and Transmissions. C177/79, pp. 9-29.

[60] Wriggers, P., 2006. Computational Contact Mechanics. Springer, Berlin, Heidelberg. doi:10.1007/978-3-540-32609-0.
[61] Yang, Y., Cao, L., Li, H., Dai, Y., 2019. Nonlinear dynamic response of a spur gear pair based on the modeling of periodic mesh stiffness and static transmission error. *Applied Mathematical Modelling* 72, 444–469. doi:https://doi.org/10.1016/j.apm.2019.03.026.

[62] Yavuz, S.D., Saribay, Z.B., Cigeroglu, E., 2018. Nonlinear time-varying dynamic analysis of a spiral bevel geared system. *Nonlinear Dynamics* 92, 1901–1919. doi:10.1007/s11071-018-4170-9.

[63] Özgüven, H., Houser, D., 1988. Dynamic analysis of high speed gears by using loaded static transmission error. *Journal of Sound and Vibration* 125, 71–83. doi:https://doi.org/10.1016/0022-460X(88)90416-6.
Nomenclature

Matrices and vectors

[∇] Frequency domain differential operator

[Ψ] Modal basis matrix

$f_{ex}$ External load vector

$f_{nl}$ Vector of nonlinear force

$q$ Vector of generalized displacements

$\bar{q}$ Vector of Fourier coefficients

$G$ Force direction vector

$R$ HBM residual

$R_{ex}$ Extended residual with arclength equation

$T$ Vector of harmonic base functions

[C] Damping matrix

$[I_n]$ Identity matrix of size n

[K] Stiffness matrix

[M] Mass matrix

[Z] HBM dynamic stiffness matrix

Scalars

$\alpha$ Pressure angle

$\phi$ Test function for bifurcation detection

$\xi_j$ Modal damping ratio of mode j

$\Omega$ Fundamental frequency

$b$ Half backlash

$b_f$ Face width

$g(t)$ Gap function

$h_a$ Addendum coefficient

$h_d$ Dedendum coefficient

$k_m(t)$ Time-varying mesh stiffness

$m$ Gear module

$m_1$ Mass of gear 1 (driving wheel)

$m_2$ Mass of gear 2 (driven wheel)

$q_s$ Static transmission error

$r_b$ Base radius

$s$ Curvilinear abscissa

$x$ Profile shift coefficient

$H$ Heaviside step function

$F_s$ Transmitted load

$F_{nl}$ Nonlinear force

$H$ harmonic truncation order

$H_m$ mesh harmonics

$I_1$ Inertia of gear 1 (driving wheel)

$I_2$ Inertia of gear 2 (driven wheel)

$I_I$ Input inertia

$I_O$ Output inertia

$K_1$ Torsional stiffness of the driving shaft

$K_2$ Torsional stiffness of the driven shaft

$K_b$ Bearing stiffness

$N_h$ Number of holes

$T$ Torque

$Z$ Number of teeth

Subscripts

$e$ Nonlinear element number

$gb$ Gear blank

$wh$ Without holes

$SD$ Standard deviation

Operators

$\otimes$ Kronecker product

$\partial_x$ Partial derivative with respect to x

Abbreviations

$DoF$ Degree of freedom

$DTE$ Dynamic transmission error

$STE$ Static transmission error
<table><caption>**Table 1**
Parameter values for the dynamic model</caption><thead><tr><th>Designation</th><th>Without holes</th><th>With holes</th><th>Unit</th></tr></thead><tbody><tr><td><i>I</i><sub>l</sub></td><td>3e-3</td><td>3e-3</td><td>kg.m<sup>2</sup></td></tr><tr><td><i>I</i><sub>o</sub></td><td>5e-2</td><td>5e-2</td><td>kg.m<sup>2</sup></td></tr><tr><td><i>K</i><sub>1</sub></td><td>4e5</td><td>4e5</td><td>N·m.rad<sup>-1</sup></td></tr><tr><td><i>K</i><sub>2</sub></td><td>4e5</td><td>4e5</td><td>N·m.rad<sup>-1</sup></td></tr><tr><td><i>K</i><sub>b</sub></td><td>1e7</td><td>1e7</td><td>N·m<sup>-1</sup></td></tr><tr><td><i>I</i><sub>1</sub></td><td>1.52e-3</td><td>1.24e-3</td><td>kg.m<sup>2</sup></td></tr><tr><td><i>I</i><sub>2</sub></td><td>1.52e-3</td><td>1.24e-3</td><td>kg.m<sup>2</sup></td></tr><tr><td><i>m</i><sub>1</sub></td><td>1.11</td><td>0.81</td><td>kg</td></tr><tr><td><i>m</i><sub>2</sub></td><td>1.11</td><td>0.81</td><td>kg</td></tr></tbody></table>

<table><caption>**Table 2**
Characteristics of the gear pair</caption><thead><tr><th>Name</th><th>Designation</th><th>Gear 1</th><th>Gear 2</th><th>Unit</th></tr></thead><tbody><tr><td>Module</td><td><i>m</i></td><td>2</td><td></td><td>mm</td></tr><tr><td>Number of teeth</td><td><i>Z</i></td><td>50</td><td>50</td><td>-</td></tr><tr><td>Pressure angle</td><td><i>α</i></td><td>20</td><td></td><td>deg</td></tr><tr><td>Base radius</td><td><i>r</i><sub>b</sub></td><td>46.984</td><td>46.984</td><td>mm</td></tr><tr><td>Profile shift coefficient</td><td><i>x</i></td><td>0</td><td>0</td><td>-</td></tr><tr><td>Addendum coefficient</td><td><i>h</i><sub>a</sub></td><td>1</td><td>1</td><td>-</td></tr><tr><td>Dedendum coefficient</td><td><i>h</i><sub>d</sub></td><td>1.25</td><td>1.25</td><td>-</td></tr><tr><td>Face width</td><td><i>b</i><sub>f</sub></td><td>20</td><td>20</td><td>mm</td></tr><tr><td colspan="5">Tip relief modification</td></tr><tr><td>Length</td><td><i>l</i></td><td>1.75</td><td>1.75</td><td>mm</td></tr><tr><td>Amount</td><td><i>a</i></td><td>5</td><td>5</td><td>μm</td></tr></tbody></table>