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

HAL
open science

# Vibro-impact dynamics of large-scale geared systems

Adrien Mélot, Joël Perret-Liaudet, Emmanuel Rigaud

► To cite this version:

Adrien Mélot, Joël Perret-Liaudet, Emmanuel Rigaud. Vibro-impact dynamics of large-scale geared systems. Nonlinear Dynamics, 2022, 10.1007/s11071-022-08144-5. http://hal-03740450v2

HAL Id: hal-03740450

https://hal.science/hal-03740450v2

Submitted on 15 Jan 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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>

Distributed under a Creative Commons Attribution 4.0 International License
# Vibro-impact dynamics of large-scale geared systems

Adrien Mélot · Joël Perret-Liaudet · Emmanuel Rigaud

Received: date / Accepted: date

**Abstract** This work is concerned with the analysis of vibro-impact responses observed in large-scale nonlinear geared systems. Emphasis is laid on the interactions between the high frequency internal excitation generated by the meshing process, i.e. the static transmission error (STE) and time-varying mesh stiffness (TVMS), and low frequency external excitations. To this end, a three-dimensional finite element (FE) model of a pump equipped with a reverse spur gear pair (gear ratio 1:1) is built. The model takes into account the flexibility of the kinematic chain, the bearings and the housing and the gear backlash nonlinearity. A reduced-order model is solved with the Harmonic Balance Method (HBM) coupled to an arc-length continuation algorithm which allows one to compute the periodic solutions of the system. The onset and disappearance of vibro-impact responses is studied through the computation of grazing bifurcations. Results show that the coupling between the external excitation and the time-varying mesh stiffness term greatly modify the characteristics of the responses in terms of number and periodicity of impacts and contact loss duration.

**Keywords** Grazing bifurcations · Harmonic balance method · Reduced-order modelling · Numerical continuation · Vibro-impacts · Gear backlash

## 1 Introduction

Gears are key components of numerous mechanical systems used to adapt and transmit speed and torques.

They can be found in automotive transmissions, wind turbines [1], geared turbofan engines [2] and pumps [3, 4], to name a few. It is well known that geared systems are prone to vibration and noise issues. The meshing process, associated to the contact between gear teeth, generate an internal excitation characterized by a transmission error and time-varying mesh stiffness. The transmission error 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 [5]. Under standard operating conditions, gear teeth are in permanent contact. However, a gap, or gear backlash, is necessary to allow for assembly and operation. Some specific operating conditions may therefore lead to contact loss and repeated impacts between gear teeth [6]. In both cases, the resulting dynamic mesh load is transmitted through the shafts to the housing. The vibratory state of the latter generates a high frequency, multi-harmonic noise, known as whining noise [7] in case of permanent contact, or rattle [8]/hammering noise [9] when impacts occur.

Most nonlinear gear dynamics studies in the literature resort to using lumped parameter, torsional models of the gear pair without considering the flexibility of other components [10, 11, 12, 13], although some consider the shafts with a lumped torsional stiffness [14, 15] and the bearings [16, 17]. Such models are limited to a single, up to a dozen of degrees of freedom (DoF) and are quite restrictive in their assumptions. More recently, a few studies used FE models with beam elements to model the shafts [18, 19]. Both types of models suffer from a limited predictability as the effects of the housing on the driveline dynamic behaviour is neglected [20]. Furthermore, these models cannot be used to retrieve the displacement field of the housing which is essential to compute the vibro-impact-induced noise.

<footnote>A. Mélot ☒ · J. Perret-Liaudet · E. Rigaud
Laboratoire de Tribologie et Dynamique des Systèmes, UMR CNRS 5513, Ecole Centrale de Lyon, 36 avenue Guy de Collongue, 69134 Ecully Cedex, France
E-mail: adrien.melot@ec-lyon.fr</footnote>
The most common numerical methods used in nonlinear gear dynamics studies are time integration techniques such as Runge-Kutta (RK) schemes [11, 21, 22] and the Harmonic Balance [23, 24, 25, 18, 15, 26]. RK-based integrators are found in numerous studies as they provide the means to investigate both periodic and non-periodic responses. However, gear dynamics are characterized by stiff equations due to the contact between gear teeth and high frequency excitations generated by the meshing process. Small time steps are thus required to ensure sufficient accuracy. Besides, taking into account the whole system usually involves considering other sources of excitation. These other sources are generally low frequency excitations, so that time integration has to be carried out over long time intervals. All these reasons, coupled to the need to compute the transient response, lead to hefty computational efforts when employing direct time integration methods.

The HBM, on the other hand, relies on solving the equations of motion in the frequency domain by approximating the solutions with Fourier series. In contrast to time integration methods, the transient response does not need to be computed and problems associated to considering both low and high frequency excitations are avoided. The accuracy of the solution and computational effort directly depend on the truncation order of the Fourier series. However, the HBM framework offers ways to efficiently reduce the number of equations in the frequency domain, e.g. either by using only a subset of harmonics, as is the case with systems with dry friction where only odd harmonics contribute to the response [27], or by dynamically adapting the truncation order [28, 29, 30].

All excitations acting on the system are also seldom taken into account [31, 4]. Garambois et al. [4] used a FE model of a kinematic chain with beam elements and a linear model of the gear mesh, i.e. without considering the backlash nonlinearity, to study the coupling between the internal excitation and the fluctuating torques. They showed that such coupling induces a spectral enrichment as well as an increased RMS value of the dynamic response resulting in a marked influence on the perceived radiated noise.

The main contributions of the present can be summarized as follows:

- To propose a general numerical strategy to compute the nonlinear dynamic response of a geared system with backlash nonlinearity. The strategy can be applied to systems of varying complexity, from single degree of freedom models to large-scale FE models including, for the first time, a detailed, high fidelity 3D model of the housing.

- To take into account complex loading scenarios, consisting of the various multiharmonic excitations stemming from different physics (contact between gear teeth, aerodynamic forces, rotor unbalance, etc.).

- To carry out an in-depth study of the coupling between these excitations, more specifically between the low frequency aerodynamic and unbalance forces and the parametric high frequency gear excitation.

The structure of the paper is as follows: Section 2 describes the FE model and reduced-order model of the studied system. In Section 3, the numerical methods used to solve the equations of motion are introduced. The results are discussed in Section 4 and Section 5 draws the main conclusions and suggests directions for future research.

## 2 Model description

### 2.1 Mathematical model

The system under investigation in this study is an industrial pump consisting of two counter-rotating shafts coupled by a spur gear pair without tooth profile modifications whose gear ratio is 1:1 with 76 teeth. The input and output shafts thus have the same rotational speed $\Omega$. The design characteristics of the gear pair are reported in Table. 1. Since their first mode is far above the frequency range of interest, each gear blank is considered as a rigid disk with lumped mass and inertia. Each gear thus possesses six degrees of freedom (three translations and three rotations) located at its centre, as depicted in Fig. 1. More details on the contact modelling are given in Sect. 2.3.

<table><caption>Table 1: Characteristics of the gear pair</caption><thead><tr><th>Name</th><th></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>1</td><td></td><td>mm</td></tr><tr><td>Number of teeth</td><td><i>Z</i></td><td>76</td><td>76</td><td>-</td></tr><tr><td>Helix angle</td><td><i>β</i></td><td>0</td><td></td><td>deg</td></tr><tr><td>Normal pressure angle</td><td><i>α</i></td><td>20</td><td></td><td>deg</td></tr><tr><td>Base radius</td><td><i>r<sub>b</sub></i></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<sub>a</sub></i></td><td>1</td><td>1</td><td>-</td></tr><tr><td>Dedendum coefficient</td><td><i>h<sub>d</sub></i></td><td>1.25</td><td>1.25</td><td>-</td></tr><tr><td>Face width</td><td><i>b<sub>f</sub></i></td><td>20</td><td></td><td>mm</td></tr></tbody></table>

Each shaft is supported by two ball bearings modelled by linear axial and radial stiffness elements. 3D quadratic tetrahedral elements are used to generate the
<figure><img src="image_3.png" /><figcaption>Fig. 1: Local frame of the gear pair.</figcaption></figure>

FE mesh of the shafts and the housing. They are described, respectively, by around 400,000 and 2.4 million degrees of freedom. Figure 3a depicts a schematic representation of the system and Figure 3b shows the geometry and mesh of the housing.

Several excitations are considered in the analysis. The external aerodynamic loading stemming from the pumping process is a narrowband random process approximated as a multi-harmonic periodic excitation. Due to the pump geometry, only even harmonics are significant. Both forces and torques are applied on twelve centreline nodes (six on each shaft) that are rigidly linked to the FE mesh of the shafts, as depicted in Fig. 3b.

$$ \mathbf{f}_{ex}^{aero} = \sum_{k=1}^{\infty} \hat{\mathbf{f}}_{ex,2k}^{aero} e^{i2k\Omega t} \quad (1) $$

In the following, two harmonic components are considered, corresponding to up to harmonic $H_4$ of the shafts rotation $\Omega$. Figure 2 shows an example of the applied aerodynamic forcing, computed in the time domain with an in-house CFD software. Its amplitude is equal to approximately 15% of the applied torque $T = 2.5$ N·m. One can also see that the forcing exhibits a phase shift equal to $\pi/2$ between the input and output shafts and are of opposite signs.

An additional excitation $f_{unb}$ in the form of a mass unbalance:

$$ f_{unb,x} = -m_{unb}\Omega^2 \sin \Omega t, \quad (2) $$

$$ f_{unb,y} = -m_{unb}\Omega^2 \cos \Omega t, \quad (3) $$

is uniformly distributed along the shafts with a phase shift equal to $\pi/2$ between the input and output shafts. The values considered is the following are $m_{unb} = 30 \cdot 10^{-6}$ kg·m for the input shaft and $m_{unb} = 20 \cdot 10^{-6}$ kg·m for the output shaft.

Considering the shaft geometry and operating speed, rotation and gyroscopic effects are neglected. The equation of motion of the system can therefore be written in the following form:

$$ \mathbf{M}\ddot{\mathbf{u}} + \mathbf{K}\mathbf{u} + \mathbf{f}_{nl}(\mathbf{u}) = \mathbf{f}_{ex}, \quad (4) $$

<figure><img src="image_4.png" /><figcaption>Fig. 2: Example of aerodynamic forcing. The blue solid line and red solid line correspond to the forcing on the input and output shafts, respectively.</figcaption></figure>

where $\mathbf{M}$ and $\mathbf{K}$ are respectively the mass and stiffness matrices and $\mathbf{u}$ is the vector of displacements of each DoF. $\mathbf{f}_{ex}$ is the vector of external periodic forcing and $\mathbf{f}_{nl}$ is the vector of nonlinear forces (cf. Section 2.3). Note that obvious time dependence is omitted in the notations for the sake of clarity.

## 2.2 Reduced-order model

The above-described FE model consists of about 3 million DoF. A direct calculation of the system's response with classical numerical strategies would lead to intractable computations. As a first step, the FE model is thus approximated by building a representative reduced-order model.

In this work, we rely on a component mode synthesis (CMS) technique [32, 33, 34]. CMS methods rely on defining a number of substructures with $n_i$ internal DoF $q_i$ and $n_c$ contact DoF $q_c$ at the interface between substructures. The internal DoF are then reduced using a set of $n_{CM}$ component modes. Since only internal DoF are reduced, CMS-based reduction methods are particularly well suited when the contact interfaces comprise a small number of DoF. In our case, the contact interface only consists of the twelve DoF used to model the gear pair (see Fig. 1) which lead to a significantly smaller model. The Craig-Bampton (CB) [35] method is employed hereafter, which consists in computing a reduction basis spanned by fixed-interface modes and static constraint modes. First, the mass and stiffness matrices are partitioned into internal and contact DoF:

$$ \mathbf{K} = \begin{bmatrix} \mathbf{K}_{ii} & \mathbf{K}_{ic} \\ \mathbf{K}_{ci} & \mathbf{K}_{cc} \end{bmatrix} \quad \text{and} \quad \mathbf{M} = \begin{bmatrix} \mathbf{M}_{ii} & \mathbf{M}_{ic} \\ \mathbf{M}_{ci} & \mathbf{M}_{cc} \end{bmatrix}. \quad (5) $$
<figure><img src="image_5.png" /><figcaption>Fig. 3: Schematic representation of the system (a) and view of the FE model of the housing (b)</figcaption></figure>

Fixed-interface modes are computed by solving the following eigenvalue problem:

$$ (\mathbf{K}_{ii} - \omega^2 \mathbf{M}_{ii}) \boldsymbol{\phi} = \mathbf{0}. \qquad (6) $$

Static constraint modes are computed by applying a unit displacement on each contact DoF, sequentially:

$$ \boldsymbol{\Psi} = -\mathbf{K}_{ii}^{-1} \mathbf{K}_{ic} \qquad (7) $$

The reduction matrix **T** is then formed by concatenating the fixed-interface and static constraint modes:

$$ \mathbf{T} = \begin{bmatrix} \boldsymbol{\Phi} & \boldsymbol{\Psi} \\ \mathbf{0} & \mathbf{I} \end{bmatrix} \qquad (8) $$

A Galerkin projection then yields $n_r = n_c + n_{CM}$ equations of motion in the reduced subspace written in matrix form:

$$ \underbrace{\mathbf{T}^T \mathbf{M} \mathbf{T} \ddot{\mathbf{q}}}_{\mathbf{M}_r} + \underbrace{\mathbf{T}^T \mathbf{K} \mathbf{T} \mathbf{q}}_{\mathbf{K}_r} + \underbrace{\mathbf{T}^T \mathbf{f}_{nl} (\mathbf{T} \mathbf{q})}_{f_{nl}^r} = \underbrace{\mathbf{T}^T \mathbf{f}_{ex}}_{f_{ex}^r} \qquad (9) $$

where $\mathbf{M}_r$, $\mathbf{K}_r$ are the mass and stiffness matrices of the reduced-order model. $f_{nl}^r$ and $f_{ex}^r$ are the vectors of nonlinear forces and external forces, respectively.

## 2.3 Contact modelling

The contact between gear teeth acts in the normal direction of the tooth profile. The displacement along the line of action describes the relative deformation and contact gap between gear teeth. This relative deformation can be retrieved from the generalized coordinates

$\mathbf{q}$ of the reduced-order model by using a projection vector $\mathbf{G}$ of size $n_r$ defined in the local frame of the gear pair (see. Fig. 1):

$$ \Delta = \mathbf{G}^T \mathbf{q}. \qquad (10) $$

Only the restriction of $\mathbf{G}$ to the gear DoF $\mathbf{q}_c = (u_1, u_2, u_3, \theta_1, \theta_2, \theta_3, u_4, u_5, u_6, \theta_4, \theta_5, \theta_6)$ has non-zero elements, i.e.:

$$ \mathbf{G}|_{q_c} = (0, 1, \tan \beta, r_{b,1} \tan \alpha \tan \beta, -r_{b,1} \tan \beta, r_{b,1}, \qquad (11) \\ 0, -1, -\tan \beta, r_{b,2} \tan \alpha \tan \beta, -r_{b,2} \tan \beta, r_{b,2})^T, $$

where $\alpha$ and $\beta$ the pressure and helix angles and $r_{b,1}$ and $r_{b,2}$ are the base radii of the pinion and the gear, respectively.

Depending on the relative displacement on the line of action, the following behaviours can be observed:

- a linear behaviour, corresponding to permanent contact between gear teeth, when the vibration amplitude is smaller than the static deflection,

- single-sided impacts when oscillations are larger than the static deflection, yet small enough not to cross the gear backlash,

- double-sided impacts when the amplitude of the dynamic response is sufficient to cross the gear backlash and lead to impacts on the reverse flank of the adjacent tooth.

The transmitted torque induces a static deflection $q_s(\theta)$ corresponding to the static transmission error (STE), which depends on the angular position of the gear pair. When the gear pair rotates, the STE generates a displacement excitation which is periodic under
stationary operating conditions. The nonlinear contact force $f_{nl}$ acting between gear teeth can be linearized around the static equilibrium, i.e. $q_s(t)$, to define a time-varying mesh stiffness $k_m(t)$. Contact loss can therefore occur at a threshold $g(t)$ (see Fig. 4):

$$g(t) = b + q_s(t) - \frac{F_s}{k_m(t)}, \qquad (12)$$

where $b$ is the constant half backlash and $F_s$ is the transmitted load.

<figure><img src="image_6.png" /><figcaption>Fig. 4: Nonlinear force model.</figcaption></figure>

The nonlinear mesh force acting along the line of action is expressed as:

$$f_{nl}^{LoA}(q, t) = k_m(t) (\mathbf{G}^T \mathbf{q} - g) \mathcal{H}(\mathbf{G}^T \mathbf{q} - g) \\ + k_m(t) (\mathbf{G}^T \mathbf{q} + g) \mathcal{H}(-\mathbf{G}^T \mathbf{q} - g), \quad (13)$$

where $\mathcal{H}$ is the Heaviside step function. The generalized nonlinear force vector is thus written as:

$$\mathbf{f}_{nl}^r(q, t) = \mathbf{G} \mathbf{f}_{nl}^{LoA}(q, t) \qquad (14)$$

The STE and mesh stiffness are computed with an in-house code based on the Reissner-Mindlin thick plate theory [36]. Due to the symmetry of the gears, an angular sector corresponding to one tooth (one mesh period) is discretized and the contact equations are solved in static conditions for each angular position.

Under stationary operating conditions, the signals are periodic and can be expressed as truncated Fourier series:

$$q_s(t) = \sum_{k=0}^{H_{q_s}} q_s^{c,k} \cos(k\Omega t) + q_s^{s,k} \sin(k\Omega t) \qquad (15)$$

$$k_m(t) = \sum_{k=0}^{H_{k_m}} k_m^{c,k} \cos(k\Omega t) + k_m^{s,k} \sin(k\Omega t) \qquad (16)$$

In the following, three harmonics are used both for the STE and mesh stiffness. Figure 5 shows the STE and mesh stiffness over two mesh period computed with an input torque $T = 2.5$ N·m. For more realistic results, tooth profile and helix errors corresponding to a quality class 6 (ISO 1328-1:2013) are taken into account in the computation.

<figure><img src="image_7.png" /><figcaption>Fig. 5: Static transmission error (a) and time-varying mesh stiffness (b) reconstructed with three harmonics at an input torque $T = 2.5$ N·m.</figcaption></figure>

## 2.4 Damping modelling

There exist numerous sources of energy dissipation and damping mechanisms in geared system, which makes them quite difficult to model accurately. Equivalent viscous damping terms acting at the gear mesh are commonly used in single degree of freedom gear models [11, 37]. In order to generalize this idea, a constant modal damping is defined with the modal basis of the underlying reduced linear system where the gears are coupled
with the average mesh stiffness $\overline{k_m}$, i.e.:

$$ \left( \left[ \mathbf{K}_r + \overline{k_m} \mathbf{G} \mathbf{G}^T \right] - \omega_{r,j}^2 \mathbf{M}_r \right) \boldsymbol{\phi}_r^{\text{contact}} = \mathbf{0}. \quad (17) $$

It is thus possible to write the modal damping matrix $\mathbf{C}_r^m$ in the reduced modal subspace:

$$ \mathbf{C}_r^m = \mathbf{diag}(2\xi_j\omega_{r,j}), \quad (18) $$

where $\omega_{r,j}$ is the eigenfrequency of the j-th mode of the reduced-order model and $\xi_j = 5\%$ in the following. The damping matrix in the reduced subspace finally reads:

$$ \mathbf{C}_r = \mathbf{M}_r^T \boldsymbol{\Phi}_r \mathbf{C}_r^{\text{mod}} \boldsymbol{\Phi}_r^T \mathbf{M}_r, \quad (19) $$

where matrix $\boldsymbol{\Phi}_r$ is formed with the mass-normalized eigenvectors of Eq. (17). Note that more sophisticated damping models, i.e. impact damping can easily be implemented since the physical gear DoFs are retained in the reduced-order model.

# 3 Computational strategy

## 3.1 Harmonic balance method

The general idea behind the HBM consists in approximating a periodic solution $\mathbf{q}$ of Eq. (4) with truncated Fourier series:

$$ \mathbf{q} = \text{Re} \left( \sum_{k=0}^{\infty} \tilde{\mathbf{q}}_k e^{ik\Omega t} \right) \approx \text{Re} \left( \sum_{k=0}^{H} \tilde{\mathbf{q}}_k e^{ik\Omega t} \right), \quad (20) $$

where $\tilde{\mathbf{q}}$ contains the coefficients of the one-sided Fourier transform. Substituting Eq. (20) into Eq. (4) and applying a Galerkin procedure yields a residual which consists of a set of $n(H + 1)$ nonlinear algebraic equations:

$$ \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 (21) $$

where $\mathbf{Z}(\Omega)$ is the dynamic stiffness matrix duplicated on each considered harmonic:

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

where $\otimes$ is the Kronecker product and $\nabla$ the frequency domain differential operator:

$$ \nabla = \mathbf{diag}(0, i, 2i, \dots, Hi). \quad (23) $$

It is worth noting that the present study only investigates period-one solutions. Subharmonic responses can be computed by modifying the Fourier base used in the HBM formulation [24,38].

The Fourier coefficients of the nonlinear forces are evaluated using Eq. (14) and the alternating frequency/time (AFT) scheme [39] whose principle is summarized in Fig. 6.

<figure><img src="image_8.png" /><figcaption>Fig. 6: Illustration of the AFT procedure.</figcaption></figure>

Eq. (21) can be solved using a Newton-Raphson-based iterative solver to yield the solution at a given frequency. In order to compute the response in a given frequency range, the frequency $\Omega$ is treated as an unknown and the response curve is parameterized by its curvilinear abscissa to avoid difficulties with vertical tangents at turning points. An arc-length continuation procedure is also implemented. The arc-length continuation is a predictor/corrector scheme (see Fig. 7) which relies on finding an approximate solution $\tilde{\mathbf{Q}}_{(p)} = (\tilde{\mathbf{q}} \ \Omega)_{(p)}^T$ based on the previously computed solution point $\tilde{\mathbf{Q}}_{(k)}$ and iteratively performing corrections on the predicted solution until convergence is reached. A thorough description of the numerical methods employed in this work can be found in [15,26].

<figure><img src="image_9.png" /><figcaption>Fig. 7: Illustration of the arc-length continuation.</figcaption></figure>
## 3.2 Condensation on the nonlinear degrees of freedom

The nonlinear force stemming from the contact between gear teeth and backlash nonlinearity only act on a small set of DoF, i.e. the twelve DoF representing the gear pair. It is possible to exploit this sparsity to only retain the nonlinear DoF as unknowns of the problem and greatly reduce the computational time [40]. The equations of motion are first partitioned into nonlinear (subscript 'nl') and linear (subscript 'ln') DoF:

$$ \begin{bmatrix} \mathbf{Z}_{nl,nl} & \mathbf{Z}_{nl,nl} \\ \mathbf{Z}_{nl,nl} & \mathbf{Z}_{nl,nl} \end{bmatrix} \begin{pmatrix} \tilde{\mathbf{q}}_{nl} \\ \tilde{\mathbf{q}}_{nl} \end{pmatrix} + \begin{pmatrix} \mathbf{0} \\ \tilde{\mathbf{F}}_{nl} \end{pmatrix} = \begin{pmatrix} \tilde{\mathbf{F}}_{ex,nl} \\ \tilde{\mathbf{F}}_{ex,nl} \end{pmatrix}, \quad (24) $$

which yields:

$$ \mathbf{Z}_r \tilde{\mathbf{q}}_{nl} + \tilde{\mathbf{F}}_{nl} = \tilde{\mathbf{F}}_r, \quad (25) $$

where $\mathbf{Z}_r$ is the Schur complement of the linear part of the dynamic stiffness matrix and $\tilde{\mathbf{F}}_r$ is the condensed external forcing vector:

$$ \begin{aligned} \mathbf{Z}_r &= \mathbf{Z}_{nl,nl} - \mathbf{Z}_{nl,nl} \mathbf{Z}_{ln,ln}^{-1} \mathbf{Z}_{ln,nl}, \\ \tilde{\mathbf{F}}_r &= \tilde{\mathbf{F}}_{ex,nl} - \mathbf{Z}_{nl,nl} \mathbf{Z}_{ln,ln}^{-1} \tilde{\mathbf{F}}_{ex,nl}. \end{aligned} \quad (26) $$

The Jacobian matrix can be expressed by direct differentiation of Eq. (26) with respect to the unknowns:

$$ \begin{aligned} \partial_{\Omega} \mathbf{Z}_r = & \partial_{\Omega} \mathbf{Z}_{nl,nl} \\ & - \left( \partial_{\Omega} \mathbf{Z}_{nl,nl} \mathbf{Z}_{ln,ln}^{-1} \mathbf{Z}_{ln,nl} + \mathbf{Z}_{nl,nl} \partial_{\Omega} \mathbf{Z}_{ln,ln}^{-1} \mathbf{Z}_{ln,nl} + \mathbf{Z}_{nl,nl} \mathbf{Z}_{ln,ln}^{-1} \partial_{\Omega} \mathbf{Z}_{ln,nl} \right), \end{aligned} \quad (27) $$

$$ \begin{aligned} \partial_{\Omega} \tilde{\mathbf{F}}_r = & \partial_{\Omega} \tilde{\mathbf{F}}_{ex,nl} \\ & - \left( \partial_{\Omega} \mathbf{Z}_{nl,nl} \mathbf{Z}_{ln,ln}^{-1} \tilde{\mathbf{F}}_{ex,nl} + \mathbf{Z}_{nl,nl} \partial_{\Omega} \mathbf{Z}_{ln,ln}^{-1} \tilde{\mathbf{F}}_{ex,nl} + \mathbf{Z}_{nl,nl} \mathbf{Z}_{ln,ln}^{-1} \partial_{\Omega} \tilde{\mathbf{F}}_{ex,nl} \right). \end{aligned} \quad (28) $$

The above expressions involve the inverse of the linear part of the dynamic stiffness matrix $\mathbf{Z}_{ln,ln}^{-1}$. This can lead to heavy computational effort as the number of linear degrees of freedom increases, even with the preliminary CMS reduction described in Sect. 2.2. To keep the computational time to a minimum, the inverse $\mathbf{Z}_{ln,ln}^{-1}$ is not computed explicitly. Instead, the terms $\partial_{\Omega} \mathbf{Z}_{nl,nl} \mathbf{Z}_{ln,ln}^{-1}$, $\mathbf{Z}_{nl,nl} \partial_{\Omega} \mathbf{Z}_{ln,ln}^{-1}$ and $\mathbf{Z}_{nl,nl} \mathbf{Z}_{ln,ln}^{-1}$ are evaluated by solving three linear systems [41]:

$$ \begin{aligned} \mathbf{Z}_{nl,nl}^T \mathbf{Z}_{nl,nl} \mathbf{Z}_{ln,ln}^{-T} &= \mathbf{Z}_{nl,nl}^T, \\ \mathbf{Z}_{ln,nl}^T \mathbf{Z}_{nl,nl} \partial_{\Omega} \mathbf{Z}_{ln,ln}^{-T} &= - \partial_{\Omega} \mathbf{Z}_{ln,nl}^T \mathbf{Z}_{nl,nl} \mathbf{Z}_{ln,ln}^{-T}, \\ \mathbf{Z}_{ln,nl}^T \partial_{\Omega} \mathbf{Z}_{nl,nl} \mathbf{Z}_{ln,ln}^{-T} &= \partial_{\Omega} \mathbf{Z}_{nl,nl}^T. \end{aligned} $$

## 3.3 Detection of grazing bifurcations

The non-smoothness of the nonlinear contact force can give rise to discontinuity-induced bifurcations [42]. Of particular importance in this paper are grazing bifurcations, characterized by impacts with zero normal velocity, which mark the onset and disappearance of vibro-impacting orbits. It is usual in bifurcation analysis to define test functions whose zeros indicate bifurcation points. In this work, grazing bifurcations are detected by monitoring the changes of sign of $\phi_{GR}$ defined as:

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

where $\langle \cdot \rangle$ denotes the mean value over one fundamental period. Before a grazing bifurcation, the gear teeth are in permanent contact and, assuming $\mathcal{H}(0) = 0$, $\phi_{GR}$ equals 1. As soon as the teeth lose contact, the value of $\phi_{GR}$ switches to -1 until contact is re-established over the whole rotation period. Thus, a grazing bifurcation corresponding to the onset (resp. disappearance) of vibro-impacting orbits is detected when $\phi_{GR}$ switches from 1 to -1 (resp. from -1 to 1). The accuracy of the detection depends directly on the number of harmonics retained in the HBM approximation and on the number of harmonics used to describe the mesh stiffness and static transmission error.

# 4 Results and discussion

## 4.1 Preliminary analysis

The accuracy of the reduced-order model is dependent on dimension of the reduced subspace, i.e. the number of fixed-interface modes used to approximate the dynamics of the linear part of the system. This number ought to be large enough to describe the vibration behaviour resulting from all excitations. The frequency range in which the reduced-order model is valid is dictated by the maximum frequency of the internal excitation, since its acts at a significantly higher frequency than the external excitation and all fixed-interface modes are considered up to truncation order $n_{CM}$.

In the following, the validity of the reduced-order model is evaluated using the concept of mesh energy $\rho$ [20] defined for each mode of the underlying linear system by Eq. (30):

$$ \rho_j = \frac{\phi_j^T [\overline{k_m \mathbf{G}} \mathbf{G}^T] \phi_j}{\phi_j^T [\mathbf{K} + \overline{k_m \mathbf{G}} \mathbf{G}^T] \phi_j}, \quad (30) $$

with:
- $\phi_j$ the eigenvector of the j-th mode,

- **G** the projection vector defined by Eq. (11),

- $[K]$ the stiffness matrix of the system without coupling between gear teeth.

Modes with a high mesh energy can be excited by the high frequency internal excitation. It is therefore necessary to ensure those modes are correctly described by the reduced-order model. In the following, we consider 3 harmonics of the internal mesh excitation, corresponding to harmonics $H_Z$, $H_{2Z}$ and $H_{3Z}$ of the shaft rotation. Figure 8 shows the Campbell diagram of the underlying linear system, i.e. with the gears coupled by the average mesh stiffness $\bar{k}_m$. Excitation orders are shown in solid black lines and the modes are plotted as horizontal lines with energy-dependent colours. For readability reasons, only modes with a mesh energy greater than 3% are plotted. One can see that several modes have a significant mesh energy. The mode with most mesh energy ($\rho = 32\%$) lies at 4 kHz. The highest frequency of a mode with mesh energy is 5 kHz. By keeping $n_{CM} = 250$ component modes in the reduction, Fig. 9 shows that the reduced-order model is valid until more than $\Omega = 9$ kHz.

<figure><img src="image_10.png" /><figcaption>Fig. 8: Campbell diagram of the underlying linear system. Excitation orders corresponding to the internal excitation are plotted in solid black lines. The colour of the horizontal lines correspond to the mesh energy of the corresponding mode. (Color figure online)</figcaption></figure>

The response being approximated by Fourier series, the accuracy directly depends on the number of retained harmonics. As such, harmonics present in the excitation must be included. The TVMS is a multi-harmonic parametric excitation which is expected to induce a coupling between the external excitation, leading to additional harmonics in the response. In our case,

the frequency of the external excitation is significantly lower (up to harmonic $H_4$ of the shaft rotation) than that of the internal excitation (harmonics $kZ$, $k \in \mathbb{N}$ of the shaft rotation). Sidebands around each considered mesh harmonic should therefore be retained in the Fourier series for an accurate description of the response.

<figure><img src="image_11.png" /><figcaption>Fig. 9: Relative error on the frequency of the modes of the linearized reduced-order model with respect to those of the full FE model.</figcaption></figure>

With large-scale systems each additional harmonic can lead to a significant increase of the computational time. Thus, a convergence study with all excitation sources is first carried out to determine the smallest set of harmonic to consider for the analysis. Figure 10 shows the evolution of the standard deviation (SD) of the mesh load, computed with Eq. (13), with respect to the rotational frequency of the shafts for several number of harmonics. The complete list can be found in table 2. One can see that convergence in term of mesh load is reached with 20 harmonics. Most features of the forced response curve are well captured once all harmonics of the internal excitation ($H_{kZ}$, $k \in \mathbb{N}$ of the shaft rotation) are included in the computation. However, it is clear that the coupling between internal and external excitations severely affects the primary resonance peak and more harmonics need to be taken into account for a satisfying accuracy.

Figure 11 depicts the frequency of all detected grazing bifurcations with respect to the number of harmonics. The vertical dashed lines correspond to a ±1% error centred on the most converged computation. It appears that carrying out the computation with an insufficient number of harmonics leads to the inaccurate detection of vibro-impact responses. For instance, the computa-
tion with 5 harmonics cannot detect the grazing bifurcations occurring at low frequencies. Furthermore, one can see that satisfying accuracy is reached with 14 harmonics, which is smaller than the number of harmonics required for a good description of the forced response curve. A set of 20 harmonics is therefore used throughout the rest of this paper. All computations were performed on a PC equipped with an AMD Ryzen 9 5950X CPU (16 cores, 32 threads) with a base clock speed equal to 3.4 GHz and 64 Gb RAM. The total computation time depends on the continuation step size, number of harmonics and sampling in the AFT procedure. With the parameters used in this paper, the average computation time is 11 s per solution point.

## 4.2 Coupling between internal and external excitations

In order to highlight the coupling between internal and external excitations in the nonlinear regime, we carry out three analyses with an input torque $T = 2.5$ N·m corresponding to a static mesh load $F_s = 70$ N:

- the first one corresponds to a forced response computation considering only the internal excitation,

- the second one considers only the external excitation (aerodynamic loading and mass unbalance) and STE with a constant mesh stiffness (the average value over one vibration cycle),

- the third one includes all excitations.

<figure><img src="image_12.png" /><figcaption>Fig. 10: Convergence of the standard deviation of the dynamic mesh load with respect to the number of retained harmonics. Diamond markers correspond to grazing bifurcation points. (Color figure online)</figcaption></figure>

Figure 12 depicts the evolution of the dynamic mesh load with respect to the rotational frequency for all

<figure><img src="image_13.png" /><figcaption>Fig. 11: Convergence of grazing bifurcations with respect to the number of retained harmonics. The vertical dashed lines correspond to a ±1% error centered on the most converged computation with 24 harmonics.</figcaption></figure>

cases. The regions where the response exhibits vibro-impacts are delimited by coloured diamond markers denoting grazing bifurcations. It appears that including low frequency external excitation lead to a higher amplitude of the dynamic response away from resonances, both at low and high frequencies. Furthermore, contrary to the resonances at $\Omega = 1085$ rpm and $\Omega = 1580$ rpm which are almost unaffected, the primary resonance peak is modified. Its amplitude decreases from 93 N when only the internal excitation is considered to 88 N with all excitations. Besides, its shape is also significantly impacted by the coupling between internal and external excitations, leading to smaller amplitude-jump instabilities.

Figure 13 shows the frequency at which grazing bifurcations occur for the considered excitation scenarios. One can see that the external excitation leads to an increase of the frequency range exhibiting vibro-impacts, as grazing bifurcations occurring at $\Omega = 1000$ rpm and $\Omega = 1300$ rpm appear more spread. Overall the span of vibro-impact regions increases by 14% when all excitations are considered.

## 4.3 Influence of the excitations on the vibro-impact characteristics

Figure 14 depicts the evolution of the standard deviation of the mesh load with respect to the rotational speed. One can see three main resonance peaks located at $\Omega = 1070$ rpm, $\Omega = 1570$ rpm and $\Omega = 3000$ rpm, corresponding to the excitation of the mode with high-
<figure><img src="image_14.png" /><figcaption>Fig. 12: Evolution of the standard deviation of the mesh load with respect to the rotational speed. Computation with all excitations (red), all excitations except the TVMS and external excitations (blue) and internal excitation only (purple). Coloured diamond markers represent grazing bifurcations.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 13: Frequency location of grazing bifurcations for all considered cases. Grazing bifurcations, represented by diamond markers, associated to both ends of a vibro-impact region are plotted in the same color.</figcaption></figure>

The mesh energy by the three harmonics of the internal excitation $H_{3Z}$, $H_{2Z}$ and $H_Z$, respectively.

Five grazing bifurcations are observed on the forced response curve. In the following, they are referred to with their number of appearance with respect to the rotational speed. For clarity, regions of the forced response curve exhibiting vibro-impacts are highlighted in blue. The first grazing bifurcation around $\Omega = 994$ rpm marks the onset of a vibro-impact region spanning $\Delta\Omega = 260$ rpm until the second grazing bifurcation at $\Omega = 1255$ rpm. A small interval of responses

with permanent contact conditions can be observed following this bifurcation until the appearance of a third grazing bifurcation at $\Omega = 1295$ rpm. The responses exhibit vibro-impacts until a speed $\Omega = 2202$ rpm where the fourth grazing bifurcation occurs, leading to another small interval without vibro-impacts. At $\Omega = 2378$ rpm, the forced response curve exhibits a grazing bifurcation responsible for the appearance of vibro-impacting responses over the rest of the frequency range of interest.

<figure><img src="image_16.png" /><figcaption>Fig. 14: Forced response curve with both internal and external excitations. Blue diamond markers denote grazing bifurcations and blue highlights denote vibro-impacting responses.</figcaption></figure>
Time responses can be evaluated from the Fourier coefficients computed with the HBM. Figure 15 depicts the time series of the dynamic mesh load computed with Eq. (13) and the dynamic transmission error computed with Eq. (10) at $\Omega = 2970$ rpm, where the most energetic mesh mode is excited by the first harmonic of the internal excitation $H_Z$. Fig. 15a and 15b show that the dynamic mesh load goes to zero, indicating a contact loss, repeatedly during the rotation period. This induces impacts on the teeth flanks. Note that the mesh load is always positive or null. This indicates that impacts always occur on the active tooth flanks and that the vibration amplitude is never high enough to cross the gear backlash. This is confirmed by Fig. 15d which shows the dynamic transmission error over two fundamental periods. The response exhibits a low frequency modulation due to the aerodynamic loads (harmonic $H_2$ of the shaft rotation, also visible on Fig. 15c) with a peak-to-peak amplitude equal to $4 \ \mu\text{m}$ which is much smaller than the gear backlash $b = 40 \ \mu\text{m}$.

Figure 15c shows the harmonic content of the dynamic mesh load at $\Omega = 2970$ rpm. One can see that the response is strongly multi-harmonic with harmonics of the mesh frequency ($H_Z$, $H_{2Z}$ and $H_{3Z}$), low frequency harmonics stemming from the aerodynamic forcing ($H_2$ and $H_4$). Sidebands generated by interactions between the aerodynamic forcing and the time-varying mesh stiffness can be observed only around mesh harmonics $H_Z$ and $H_{2Z}$. Sidebands around harmonics $H_{3Z}$ and $H_{4Z}$ do not appear for they were not considered in the computation. Note that the amplitude of the sidebands is of the same order of magnitude as the $H_2$ harmonic.

Figure 15e is a close-up of Fig. 15d. The black dashed line represents the gap limit corresponding to the teeth flanks with the STE. It appears clearly that impacts occur once every fundamental period of the STE, i.e. 76 times per shaft rotation. Note that due to the low frequency modulation, the impact velocity varies periodically. This results in significant variations of the dynamic mesh load, as evidenced in Fig. 15b. The smallest impacts generate a mesh load equal to 180 N while the largest ones generate a mesh load equal to 260 N. This corresponds to a variation of 44% of the dynamic mesh load occurring twice per shaft rotation.

Figure 16 shows the time series and harmonic content of the dynamic mesh load computed with Eq. (13) and the time series of the dynamic transmission error computed with Eq. (10) at $\Omega = 6875$ rpm. At this rotational speed, the frequency of the internal excitation is too high to excite the mesh modes. The small dynamic amplification visible on Fig. 14 is induced by the excitation of a shaft mode by the external excitation. One can see from Fig. 16a that the gear response is governed

by the second harmonic of the shaft rotation induced by the aerodynamic forcing, with the amplitude of the mesh harmonics relatively small. The peak-to-peak amplitude of the mesh load is significantly smaller than the one at $\Omega = 2970$ rpm (160 N instead of 275 N).

This is confirmed by Fig. 16c which shows the harmonic content of the dynamic mesh load at $\Omega = 6875$ rpm. Compared to the results at $\Omega = 2970$ rpm, harmonics $H_2$ and $H_4$ appear much more prominent, with harmonic $H_2$ of the same order of magnitude as the mesh harmonics. Furthermore, the amplitude sidebands around each mesh harmonic are quite noticeably smaller.

Fig. 16e shows that vibro-impacts still occur once per mesh period with 76 impacts per shaft rotation. However, Fig. 16a clearly shows intermittent phases of vibro-impacts induced by the low frequency modulation. Phases of permanent contact lasting several mesh periods ($\approx 22\%$ of the shaft rotation) alternate with phases of vibro-impacts, lasting approximately 28% of the shaft rotation, twice per shaft rotation.

The same analysis can be carried out over the whole forced response curve. Figures 17a and 17b shows two close-ups of the time series of the dynamic transmission error at $\Omega = 1700$ rpm. At this rotational speed, the external excitation-induced modulation does not lead to alternating phases of permanent contact and phases with 1 impact per period of the internal excitation (1-IPP). Instead, the modulation induces a change in the number of impacts per period during the rotation of the shafts. The dynamic response exhibits 1-IPP during around 20% of the shaft rotation and 2-IPP for the rest of the rotation. Note that the change in the number of IPP is governed by the external excitation so that it occurs twice per shaft rotation in the following sequence:

- 2-IPP (40% of the rotation),

- 1-IPP (10% of the rotation),

- 2-IPP (40% of the rotation),

- 1-IPP (10% of the rotation).

The response at $\Omega = 2970$ rpm exhibits an even more complex vibro-impacting behaviour. Figure 18 shows close-ups of the time series of the dynamic transmission error at $\Omega = 1070$ rpm, corresponding to the resonance of the most energetic mesh mode excited by the third harmonic of the internal excitation $H_{3Z}$.

One can see that, since the response is governed by the third harmonic of the internal excitation, the response exhibits sequentially permanent contact, 1-IPP, 2-IPP and 3-IPP with continuous, smooth transitions between each type of behaviour in the following sequence:
<figure><img src="image_17.png" /><figcaption>Fig. 15: Time series of the dynamic mesh load (a, b) and dynamic transmission error (d, e) and harmonic content of the dynamic mesh load (c) at $\Omega = 2970$ rpm. Close-ups of the dynamic mesh load (b) and dynamic transmission error (e). Black dashed lines correspond to the gap limit, i.e. teeth flanks.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>Fig. 16: Time series of the dynamic mesh load (a, b) and dynamic transmission error (d, e) and harmonic content of the dynamic mesh load (c) at $\Omega = 6875$ rpm. Close-ups of the dynamic mesh load (b) and dynamic transmission error (e). Black dashed lines correspond to the gap limit, i.e. teeth flanks.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>Fig. 17: Time series of the dynamic transmission error at $\Omega = 1700$ rpm where the response exhibits one impact per period (a) and two impacts per period (b). The black dashed line corresponds to the gap limit.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Fig. 18: Time series of the dynamic transmission error at $\Omega = 1070$ rpm where the response exhibits permanent contact (a), one impact per period (b), two impacts per period (c) and three impacts per period (d). The black dashed line corresponds to the gap limit.</figcaption></figure>
- Permanent contact (14.5% of the rotation),

- 1-IPP (4% of the rotation),

- 2-IPP (5.25% of the rotation),

- 3-IPP (17% of the rotation),

- 2-IPP (5.25% of the rotation),

- 1-IPP (4% of the rotation),

- Permanent contact (14.5% of the rotation),

- 1-IPP (4% of the rotation),

- 2-IPP (5.25% of the rotation),

- 3-IPP (17% of the rotation),

- 2-IPP (5.25% of the rotation),

- 1-IPP (4% of the rotation).

## 5 Concluding remarks

Gears are an essential component of a plethora of mechanical systems present in almost all industrial sectors. A thorough understanding of their dynamic behaviour is thus of crucial importance to achieve efficient designs to meet the ever more stringent demands of the industry. Such demands often require that all components be included in the mechanical analyses to obtain representative models. Besides, mechanical systems are often a significant source of noise that ought to be reduced as much as possible. Thus, it is of vital importance to model the flexibility of the housing and to employ numerical methods able to deal with the considerable increase of the dimension of the system.

In this regard, this paper proposed a computational strategy to study the periodic solutions of large-scale geared system. It is based on a spectral technique, the harmonic balance method, coupled to an arc-length continuation algorithm. The industrial-scale model is reduced in two steps. Firstly, the Craig-Bampton method allows one to a few hundreds of degrees of freedom. Secondly, the number of equations in the spectral domain is reduced to twelve via an exact explicit condensation on the gear degrees of freedom. The proposed methodology is able to compute the periodic solutions of high-fidelity models of geared systems composed of several millions of degrees of freedom. The dynamic response at the gear can then be used to reconstruct the response of all degrees of freedom, allowing one to retrieve the dynamic response of the housing which can be used to evaluate the radiated sound.

Besides, retaining the gear degrees of freedom as physical coordinates with the Craig-Bampton reduction allows one to modify the gear and carry out an optimization of the gear design parameters to minimize the dynamic response without having to compute an updated reduced order model.

The main findings of this paper show that the parametric mesh stiffness induces a coupling with the low

frequency external excitation generated by the pumping process of the pump under study. This coupling induces an enrichment of the harmonic content of the dynamic transmission error with sidebands around each gear mesh harmonic. As a result, depending on the rotational speed and dynamic amplifications due to mesh modes and/or shaft modes, once can observe either alternating phases of permanent contact and vibro-impacts or transitions between responses with 1 impact per period and 2 impacts per period. Future work will focus on computing the radiated noise and evaluating the effect of these modulations on the acoustic perception.

**Acknowledgements** This work was performed within the framework of the LabCom LADAGE (LAboratoire de Dynamique des engrenAGEs), 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).

## Conflict of interest

The authors declare that they have no conflict of interest.

## Data availability

The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request.

## References

1. Nejad, A.R., Keller, J., Guo, Y., Sheng, S., Polinder, H., Watson, S., Dong, J., Qin, Z., Ebrahimi, A., Schelenz, R., Gutiérrez Guzmán, F., Cornel, D., Golafshan, R., Jacobs, G., Blockmans, B., Bosmans, J., Pluymers, B., Carroll, J., Koukoura, S., Hart, E., McDonald, A., Natarajan, A., Torsvik, J., Moghadam, F.K., Daems, P.J., Verstraeten, T., Peeters, C., Helsen, J.: Wind turbine drivetrains: state-of-the-art technologies and future development trends. Wind Energy Science 7(1), 387–411 (2022). DOI https://doi.org/10.5194/wes-7-387-2022

2. Tatar, A., Schwingshackl, C.W., Friswell, M.I.: Dynamic behaviour of three-dimensional planetary geared rotor systems. Mechanism and Machine Theory 134, 39–56 (2019). DOI https://doi.org/10.1016/j.mechmachtheory.2018.12.023

3. Mason, J., Homer, M., Eddie Wilson, R.: Mathematical models of gear rattle in roots blower vacuum pumps. Journal of Sound and Vibration 308(3), 431–440 (2007). DOI https://doi.org/10.1016/j.jsv.2007.03.071
4. Garambois, P., Donnard, G., Rigaud, E., Perret-Liaudet, J.: 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 (2017). DOI 10.1016/j.jsv.2017.05.043

5. Welbourn, D.: Fundamental knowledge of gear noise: A survey. In: Proceedings of conf. on Noise and Vibrations of Engines and Transmissions. C177/79, pp. 9-29. (1979)

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

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

8. Karagiannis, K., Pfeiffer, F.: Theoretical and experimental investigations of gear-rattling. Nonlinear Dynamics **2**(5), 367–387 (1991). DOI 10.1007/BF00045670

9. Pfeiffer, F., Prestl, W.: Hammering in diesel-engine driveline systems. Nonlinear Dynamics **5**(4), 477–492 (1994). DOI 10.1007/BF00052455

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

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

12. Margielewicz, J., Gąska, D., Litak, G.: Modelling of the gear backlash. Nonlinear Dynamics **97**(1), 355–368 (2019). DOI https://doi.org/10.1007/s11071-019-04973-z

13. Cao, Z., Chen, Z., Jiang, H.: Nonlinear dynamics of a spur gear pair with force-dependent mesh stiffness. Nonlinear Dynamics **99**(2), 1227–1241 (2020). DOI https://doi.org/10.1007/s11071-019-05348-0

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

15. Mélot, A., Benaïcha, Y., Rigaud, E., Perret-Liaudet, J., Thouverez, F.: Effect of gear topology discontinuities on the nonlinear dynamic response of a multi-degree-of-freedom gear train. Journal of Sound and Vibration **516**, 116495 (2022). DOI https://doi.org/10.1016/j.jsv.2021.116495

16. Shin, D., Palazzolo, A.: Nonlinear analysis of a geared rotor system supported by fluid film journal bearings. Journal of Sound and Vibration **475**, 115269 (2020). DOI https://doi.org/10.1016/j.jsv.2020.115269

17. Azimi, M.: Pitchfork and Hopf bifurcations of geared systems with nonlinear suspension in permanent contact regime. Nonlinear Dynamics **107**(4), 3339–3363 (2022). DOI https://doi.org/10.1007/s11071-021-07110-x

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

19. Yavuz, S.D., Saribay, Z.B., Cigeroglu, E.: Nonlinear dynamic analysis of a drivetrain composed of spur, helical and spiral bevel gears. Nonlinear Dynamics **100**(4), 3145–3170 (2020). DOI https://doi.org/10.1007/s11071-020-05666-8

20. Rigaud, E., Sabot, J.: Effect of Elasticity of Shafts, Bearings, Casing and Couplings on the Critical Rotational Speeds of a Gearbox. In: International Conference on Gears, VDI Berichte, vol. 1230, pp. 833–845. Dresde, Germany (1996)

21. Byrtus, M., Zeman, V.: On modeling and vibration of gear drives influenced by nonlinear couplings. Mechanism and Machine Theory **46**(3), 375–397 (2011). DOI https://doi.org/10.1016/j.mechmachtheory.2010.10.007

22. Pan, W., Li, X., Wang, L., Yang, Z.: Nonlinear response analysis of gear-shaft-bearing system considering tooth contact temperature and random excitations. Applied Mathematical Modelling **68**, 113–136 (2019). DOI https://doi.org/10.1016/j.apm.2018.10.022

23. Al-shyyab, A., Kahraman, A.: 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**(1-2), 151–172 (2005). DOI https://doi.org/10.1016/j.jsv.2004.06.010

24. Al-shyyab, A., Kahraman, A.: 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**(1-2), 417–451 (2005). DOI https://doi.org/10.1016/j.jsv.2003.11.029

25. Yoon, J.Y., Kim, B.: Effect and feasibility analysis of the smoothening functions for clearance-type nonlinearity in a practical driveline system. Nonlinear Dynamics **85**(3), 1651–1664 (2016). DOI https://doi.org/10.1007/s11071-016-2784-3

26. Mélot, A., Rigaud, E., Perret-Liaudet, J.: Bifurcation tracking of geared systems with parameter-dependent internal excitation. Nonlinear Dynamics **107**(1), 413–431 (2022). DOI https://doi.org/10.1007/s11071-021-07018-6

27. Pierre, C., Ferri, A.A., Dowell, E.H.: Multi-Harmonic Analysis of Dry Friction Damped Systems Using an Incremental Harmonic Balance Method. Journal of Applied Mechanics **52**(4), 958–964 (1985). DOI https://doi.org/10.1115/1.3169175

28. Grolet, A., Thouverez, F.: On a new harmonic selection technique for harmonic balance method. Mechanical Systems and Signal Processing **30** (2012). DOI https://doi.org/10.1016/j.ymssp.2012.01.024

29. Süß, D., Jerschl, M., Willner, K.: Adaptive harmonic balance analysis of dry friction damped systems, pp. 405–414 (2016). DOI https://doi.org/10.1007/978-3-319-29739-2_36

30. Gastaldi, C., Berruti, T.M.: A method to solve the efficiency-accuracy trade-off of multi-harmonic balance calculation of structures with friction contacts. International Journal of Non-Linear Mechanics **92**, 25–40 (2017). DOI https://doi.org/10.1016/j.ijnolinmec.2017.03.010

31. Ottewill, J.R., Neild, S.A., Wilson, R.E.: Intermittent gear rattle due to interactions between forcing and manufacturing errors. Journal of Sound and Vibration **321**(3), 913–935 (2009). DOI https://doi.org/10.1016/j.jsv.2008.09.050

32. Battiato, G., Firrone, C., Berruti, T., Epureanu, B.: Reduction and coupling of substructures via gram-schmidt interface modes. Computer Methods in Applied Mechanics and Engineering **336**, 187–212 (2018). DOI https://doi.org/10.1016/j.cma.2018.03.001

33. Yuan, J., El-Haddad, F., Salles, L., Wong, C.: Numerical Assessment of Reduced Order Modeling Techniques for Dynamic Analysis of Jointed Structures With Contact Nonlinearities. Journal of Engineering for Gas Turbines
and Power **141**(3) (2018). DOI https://doi.org/10.1115/1.4041147. 031027

34. Yuan, J., Schwingshackl, C., Wong, C., Salles, L.: On an improved adaptive reduced-order model for the computation of steady-state vibrations in large-scale non-conservative systems with friction joints. *Nonlinear Dynamics* **103**(4), 3283–3300 (2021). DOI https://doi.org/10.1007/s11071-020-05890-2

35. Craig, R.R., Bampton, M.C.C.: Coupling of substruc- tures for dynamic analyses. AIAA Journal **6**(7), 1313– 1319 (1968). DOI https://doi.org/10.2514/3.4741

36. Garambois, P., Perret-Liaudet, J., Rigaud, E.: NVH robust optimization of gear macro and microgeometries us- ing an efficient tooth contact model. *Mechanism and Machine Theory* **117**, 78 – 95 (2017). DOI https://doi.org/10.1016/j.mechmachtheory.2017.07.008

37. Farshidianfar, A., Saghafi, A.: Global bifurcation and chaos analysis in nonlinear vibration of spur gear sys- tems. *Nonlinear Dynamics* **75**(4), 783–806 (2014). DOI https://doi.org/10.1007/s11071-013-1104-4

38. Alcorta, R., Baguet, S., Prabel, B., Piteau, P., Jacquet- Richardet, G.: Period doubling bifurcation analysis and isolated sub-harmonic resonances in an oscilla- tor with asymmetric clearances. *Nonlinear Dynamics* **98**(4), 2939–2960 (2019). DOI https://doi.org/10.1007/s11071-019-05245-6

39. Cameron, T.M., Griffin, J.: An Alternating Fre- quency/Time Domain Method for Calculating the Steady-State Response of Nonlinear Dynamic Systems. *Journal of Applied Mechanics*. (1989). DOI https://doi.org/10.1115/1.3176036

40. Petrov, E.P.: A high-accuracy model reduction for anal- ysis of nonlinear vibrations in structures with contact interfaces. *Journal of Engineering for Gas Turbines and Power* **133**(10) (2011). DOI 10.1115/1.4002810

41. Colaïtis, Y.: Stratégie numérique pour l'analyse qualita- tive des interactions aube/carter. Ph.D. thesis, Polytech- nique Montréal (2021)

42. Ibrahim, R.A.: Vibro-Impact Dynamics. Springer-Verlag Berlin Heidelberg (2009). DOI https://doi.org/10.1007/978-3-642-00275-5

Appendix

<table><caption>Table 2: List of harmonics used in the convergence study</caption><thead><tr><td>Nb of harmonics</td><td>List of harmonics</td></tr></thead><tbody><tr><td>5</td><td>(0, 1, 2, 4, 76)</td></tr><tr><td>6</td><td>(0, 1, 2, 4, 76, 152)</td></tr><tr><td>7</td><td>(0, 1, 2, 4, 76, 152, 228)</td></tr><tr><td>8</td><td>(0, 1, 2, 4, 76, 152, 228, 304)</td></tr><tr><td>10</td><td>(0, 1, 2, 4, 74, 76, 78, 152, 228, 304)</td></tr><tr><td>12</td><td>(0, 1, 2, 4, 74, 76, 78, 150, 152, 154, 228, 304)</td></tr><tr><td>14</td><td>(0, 1, 2, 4, 74, 75, 76, 77, 78, 150, 152, 228, 304)</td></tr><tr><td>18</td><td>(0, 1, 2, 4, 72, 73, 74, 75, 76, 77, 78, 79, 80, 150, 152, 154, 228, 304)</td></tr><tr><td>20</td><td>(0, 1, 2, 4, 72, 73, 74, 75, 76, 77, 78, 79, 80, 148, 150, 152, 154, 156, 228, 304)</td></tr><tr><td>24</td><td>(0, 1, 2, 4, 72, 73, 74, 75, 76, 77, 78, 79, 80, 148, 150, 152, 154, 156, 224, 226, 228, 230, 232, 304)</td></tr></tbody></table>