Adrien Mélot
email: adrien.melot@ec-lyon.fr

Joël Perret-Liaudet

Emmanuel Rigaud

Vibro-impact dynamics of large-scale geared systems

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

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.

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 [START_REF] Nejad | Wind turbine drivetrains: state-of-the-art technologies and future development trends[END_REF], geared turbofan engines [START_REF] Tatar | Dynamic behaviour of three-dimensional planetary geared rotor systems[END_REF] and pumps [START_REF] Mason | Mathematical models of gear rattle in roots blower vacuum pumps[END_REF][START_REF] Garambois | Multiphysics coupling between periodic gear mesh excitation and input/output fluctuating torques: Application to a roots vacuum pump[END_REF], 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 [START_REF] Welbourn | Fundamental knowledge of gear noise: A survey[END_REF]. 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 [START_REF] Rigaud | Investigation of gear rattle noise including visualization of vibro-impact regimes[END_REF]. 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, multiharmonic noise, known as whining noise [START_REF] Carbonelli | Vibro-Acoustic Analysis of Geared Systems-Predicting and Controlling the Whining Noise[END_REF] in case of permanent contact, or rattle [START_REF] Karagiannis | Theoretical and experimental investigations of gear-rattling[END_REF]/hammering noise [START_REF] Pfeiffer | Hammering in diesel-engine driveline systems[END_REF] 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 [START_REF] Nevzat Özgüven | Mathematical models used in gear dynamics-a review[END_REF][START_REF] Kahraman | Non-linear dynamics of a spur gear pair[END_REF][START_REF] Margielewicz | Modelling of the gear backlash[END_REF][START_REF] Cao | Nonlinear dynamics of a spur gear pair with force-dependent mesh stiffness[END_REF], although some consider the shafts with a lumped torsional stiffness [START_REF] Liu | Investigation of nonlinear characteristics of the motor-gear transmission system by trajectory-based stability preserving dimension reduction methodology[END_REF][START_REF] Mélot | Effect of gear topology discontinuities on the nonlinear dynamic response of a multi-degree-offreedom gear train[END_REF] and the bearings [START_REF] Shin | Nonlinear analysis of a geared rotor system supported by fluid film journal bearings[END_REF][START_REF] Azimi | Pitchfork and Hopf bifurcations of geared systems with nonlinear suspension in permanent contact regime[END_REF]. 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 [START_REF] Yavuz | Nonlinear timevarying dynamic analysis of a spiral bevel geared system[END_REF][START_REF] Yavuz | Nonlinear dynamic analysis of a drivetrain composed of spur, helical and spiral bevel gears[END_REF]. Both types of models suffer from a limited predictability as the effects of the housing on the driveline dynamic behaviour is neglected [START_REF] Rigaud | Effect of Elasticity of Shafts, Bearings, Casing and Couplings on the Critical Rotational Speeds of a Gearbox[END_REF]. Furthermore, these models cannot be used to retrieve the displacement field of the housing which is essential to compute the vibro-impact-induced noise.

The most common numerical methods used in nonlinear gear dynamics studies are time integration techniques such as Runge-Kutta (RK) schemes [START_REF] Kahraman | Non-linear dynamics of a spur gear pair[END_REF][START_REF] Byrtus | On modeling and vibration of gear drives influenced by nonlinear couplings[END_REF][START_REF] Pan | Nonlinear response analysis of gear-shaft-bearing system considering tooth contact temperature and random excitations[END_REF] and the Harmonic Balance [START_REF] Al-Shyyab | Non-linear dynamic analysis of a multi-mesh gear train using multi-term harmonic balance method: period-one motions[END_REF][START_REF] Al-Shyyab | Non-linear dynamic analysis of a multi-mesh gear train using multi-term harmonic balance method: sub-harmonic motions[END_REF][START_REF] Yoon | Effect and feasibility analysis of the smoothening functions for clearance-type nonlinearity in a practical driveline system[END_REF][START_REF] Yavuz | Nonlinear timevarying dynamic analysis of a spiral bevel geared system[END_REF][START_REF] Mélot | Effect of gear topology discontinuities on the nonlinear dynamic response of a multi-degree-offreedom gear train[END_REF][START_REF] Mélot | Bifurcation tracking of geared systems with parameterdependent internal excitation[END_REF]. RKbased integrators are found in numerous studies as they provide the means to investigate both periodic and nonperiodic 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 [START_REF] Pierre | Multi-Harmonic Analysis of Dry Friction Damped Systems Using an Incremental Harmonic Balance Method[END_REF], or by dynamically adapting the truncation order [START_REF] Grolet | On a new harmonic selection technique for harmonic balance method[END_REF][START_REF] Süß | Adaptive harmonic balance analysis of dry friction damped systems[END_REF][START_REF] Gastaldi | A method to solve the efficiency-accuracy trade-off of multi-harmonic balance calculation of structures with friction contacts[END_REF].

All excitations acting on the system are also seldom taken into account [START_REF] Ottewill | Intermittent gear rattle due to interactions between forcing and manufacturing errors[END_REF][START_REF] Garambois | Multiphysics coupling between periodic gear mesh excitation and input/output fluctuating torques: Application to a roots vacuum pump[END_REF]. Garambois et al. [START_REF] Garambois | Multiphysics coupling between periodic gear mesh excitation and input/output fluctuating torques: Application to a roots vacuum pump[END_REF] 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 scenarii, 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.

Model description

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 Ω. 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. 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.

f aero ex = ∞ k=1 f aero ex,2k e i2kΩt (1)
In the following, two harmonic components are considered, corresponding to up to harmonic H 4 of the shafts rotation Ω. 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 π/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 Ω 2 sin Ωt, (2)
f unb,y = -m unb Ω 2 cos Ωt, (3)
is uniformly distributed along the shafts with a phase shift equal to π/2 between the input and output shafts. The values considered is the following are m unb = 30 • 10 -6 kg•m for the input shaft and m unb = 20•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: where M and K are respectively the mass and stiffness matrices and u is the vector of displacements of each DoF. f ex is the vector of external periodic forcing and 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.

Mü + Ku + f nl (u) = f ex , (4)

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 [START_REF] Battiato | Reduction and coupling of substructures via gram-schmidt interface modes[END_REF][START_REF] Yuan | Numerical Assessment of Reduced Order Modeling Techniques for Dynamic Analysis of Jointed Structures With Contact Nonlinearities[END_REF][START_REF] Yuan | On an improved adaptive reduced-order model for the computation of steady-state vibrations in large-scale nonconservative systems with friction joints[END_REF]. 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) [START_REF] Craig | Coupling of substructures for dynamic analyses[END_REF] 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: Fixed-interface modes are computed by solving the following eigenvalue problem:

K = K ii K ic K ci K cc and M = M ii M ic M ci M cc . (5)
K ii -ω 2 M ii ϕ ϕ ϕ = 0. (6)
Static constraint modes are computed by applying a unit displacement on each contact DoF, sequentially:

Ψ = -K -1 ii K ic (7)
The reduction matrix T is then formed by concatenating the fixed-interface and static constraint modes:

T = Φ Φ Φ Ψ Ψ Ψ 0 I (8)
A Galerkin projection then yields n r = n c + n CM equations of motion in the reduced subspace written in matrix form:

T T MT Mr q + T T KT Kr q + T T f nl f r nl (Tq) = T T f ex f r ex ( 9
)
where M r , K r are the mass and stiffness matrices of the reduced-order model. f r nl and f r ex are the vectors of nonlinear forces and external forces, respectively.

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 q of the reduced-order model by using a projection vector G of size n r defined in the local frame of the gear pair (see. Fig. 1):

∆ ∆ ∆ = G T q. ( 10
)
Only the restriction of G to the gear DoF

q c = (u 1 , u 2 , u 3 , θ 1 , θ 2 , θ 3 , u 4 , u 5 , u 6 , θ 4 , θ 5 , θ 6
) has nonzero elements, i.e.:

G| qc = (0, 1, tan β, r b,1 tan α tan β, -r b,1 tan β, r b,1 , (11)
0, -1, -tan β, r b,2 tan α tan β, -r b,2 tan β, r b,2 ) T ,
where α and β 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 (θ) 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 timevarying mesh stiffness k m (t). Contact loss can therefore occur at a threshold g(t) (see Fig. 4):

g(t) = b + q s (t) - F s k m (t) , ( 12
)
where b is the constant half backlash and F s is the transmitted load. The nonlinear mesh force acting along the line of action is expressed as:

f LoA nl (q, t) = k m (t) G T q -g H G T q -g + k m (t) G T q + g H -G T q -g , ( 13
)
where H is the Heaviside step function. The generalized nonlinear force vector is thus written as:

f r nl (q, t) = Gf LoA nl (q, t) (14)
The STE and mesh stiffness are computed with an inhouse code based on the Reissner-Mindlin thick plate theory [START_REF] Garambois | NVH robust optimization of gear macro and microgeometries using an efficient tooth contact model[END_REF]. 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) = Hq s k=0 q c,k s cos(kΩt) + q s,k s sin(kΩt) (15)
k m (t) = H km k=0 k c,k m cos(kΩt) + k s,k m sin(kΩt) (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.

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 [START_REF] Kahraman | Non-linear dynamics of a spur gear pair[END_REF][START_REF] Farshidianfar | Global bifurcation and chaos analysis in nonlinear vibration of spur gear systems[END_REF]. 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 k m , i.e.:

K r + k m GG T -ω 2 r,j M r ϕ ϕ ϕ contact r = 0. (17)
It is thus possible to write the modal damping matrix C m r in the reduced modal subspace:

C m r = diag (2ξ j ω r,j ) , (18)
where ω r,j is the eigenfrequency of the j-th mode of the reduced-order model and ξ j = 5% in the following. The damping matrix in the reduced subspace finally reads:

C r = M T r Φ r C mod r Φ T r M r , (19)
where matrix Φ 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

Harmonic balance method

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

q = Re ∞ k=0 qk e ikΩt ≈ Re H k=0 qk e ikΩt , (20)
where 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:

R(q, Ω) = Z(Ω)q + fnl (q) -fex = 0, (21)
where Z(Ω) is the dynamic stiffness matrix duplicated on each considered harmonic:

Z(Ω) = Ω 2 ∇ 2 ⊗ M + Ω∇ ⊗ C + I H+1 ⊗ K. ( 22
)
where ⊗ is the Kronecker product and ∇ the frequency domain differential operator:

∇ = diag(0, i, 2i, . . . , Hi). ( 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 [START_REF] Al-Shyyab | Non-linear dynamic analysis of a multi-mesh gear train using multi-term harmonic balance method: sub-harmonic motions[END_REF][START_REF] Alcorta | Period doubling bifurcation analysis and isolated sub-harmonic resonances in an oscillator with asymmetric clearances[END_REF]. The Fourier coefficients of the nonlinear forces are evaluated using Eq. ( 14) and the alternating frequency/time (AFT) scheme [START_REF] Cameron | An Alternating Frequency/Time Domain Method for Calculating the Steady-State Response of Nonlinear Dynamic Systems[END_REF] whose principle is summarized in Fig. 6 Eq. ( 21) can be solved using a Newton-Raphsonbased iterative solver to yield the solution at a given frequency. In order to compute the response in a given frequency range, the frequency Ω 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

Q(p) = (q Ω) T (p)
based on the previously computed solution point 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 [START_REF] Mélot | Effect of gear topology discontinuities on the nonlinear dynamic response of a multi-degree-offreedom gear train[END_REF][START_REF] Mélot | Bifurcation tracking of geared systems with parameterdependent internal excitation[END_REF]. Fig. 7: Illustration of the arc-length continuation.

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 [START_REF] Petrov | A high-accuracy model reduction for analysis of nonlinear vibrations in structures with contact interfaces[END_REF]. The equations of motion are first partitioned into nonlinear (subscript 'nl') and linear (subscript 'ln') DoF:

Z ln,ln Z ln,nl Z nl,ln Z nl,nl qln qnl + 0 Fnl = Fex,ln Fex,nl , (24)
which yields:

Z r qnl + Fnl = Fr , (25)
where Z r is the Schur complement of the linear part of the dynamic stiffness matrix and Fr is the condensed external forcing vector:

Z r = Z nl,nl -Z nl,ln Z -1 ln,ln Z ln,nl , Fr = Fex,nl -Z nl,ln Z -1 ln,ln Fex,ln . (26)
The Jacobian matrix can be expressed by direct differentiation of Eq. ( 26) with respect to the unknowns:

∂ Ω Z r =∂ Ω Z nl,nl -∂ Ω Z nl,ln Z -1 ln,ln Z ln,nl + Z nl,ln ∂ Ω Z -1 ln,ln Z ln,nl +Z nl,ln Z -1 ln,ln ∂ Ω Z ln,nl , (27)
∂ Ω Fr = ∂ Ω Fex,nl -∂ Ω Z nl,ln Z -1 ln,ln Fex,ln + Z nl,ln ∂ Ω Z -1 ln,ln Fex,ln + Z nl,ln Z -1 ln,ln ∂ Ω Fex,ln . (28)
The above expressions involve the inverse of the linear part of the dynamic stiffness matrix Z -1 ln,ln . 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 Z -1 ln,ln is not computed explicitly. Instead, the terms ∂ Ω Z nl,ln Z -1 ln,ln , Z nl,ln ∂ Ω Z -1 ln,ln and Z nl,ln Z -1 ln,ln are evaluated by solving three linear systems [START_REF] Colaïtis | Stratégie numérique pour l'analyse qualitative des interactions aube/carter[END_REF]:

Z T ln,ln Z nl,ln Z -T ln,ln = Z T nl,ln , Z T ln,ln Z nl,ln ∂ Ω Z -T ln,ln = -∂ Ω Z T ln,ln Z nl,ln Z -T ln,ln , Z T ln,ln ∂ Ω Z nl,ln Z -T ln,ln = ∂ Ω Z T nl,ln .

Detection of grazing bifurcations

The non-smoothness of the nonlinear contact force can give rise to discontinuity-induced bifurcations [START_REF] Ibrahim | Vibro-Impact Dynamics[END_REF]. Of particular importance in this paper are grazing bifurcations, characterized by impacts with zero normal velocity, which mark the onset and disappearance of vibroimpacting 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 ϕ GR defined as:

ϕ GR = 1 -2 H 1 -H G T q -g(t) , (29)
where ⟨•⟩ denotes the mean value over one fundamental period. Before a grazing bifurcation, the gear teeth are in permanent contact and, assuming H(0) = 0, ϕ GR equals 1. As soon as the teeth lose contact, the value of ϕ 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 ϕ 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.

Results and discussion

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 fixedinterface 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 ρ [START_REF] Rigaud | Effect of Elasticity of Shafts, Bearings, Casing and Couplings on the Critical Rotational Speeds of a Gearbox[END_REF] defined for each mode of the underlying linear system by Eq. ( 30):

ρ j = ϕ ϕ ϕ T j k m GG T ϕ ϕ ϕ j ϕ ϕ ϕ T j K + k m GG T ϕ ϕ ϕ j , (30) with:
ϕ ϕ ϕ 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 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 (ρ = 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 Ω = 9 kHz. 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 multiharmonic 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 ∈ 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. 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 ∈ 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.

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. cases. The regions where the response exhibits vibroimpacts 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 Ω = 1085 rpm and Ω = 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 scenarii. One can see that the external excitation leads to an increase of the frequency range exhibiting vibro-impacts, as grazing bifurcations occurring at Ω = 1000 rpm and Ω = 1300 rpm appear more spread. Overall the span of vibro-impact regions increases by 14% when all excitations are considered.

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 Ω = 1070 rpm, Ω = 1570 rpm and Ω = 3000 rpm, corresponding to the excitation of the mode with high- est 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 Ω = 994 rpm marks the onset of a vibro-impact region spanning ∆Ω = 260 rpm until the second grazing bifurcation at Ω = 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 Ω = 1295 rpm. The responses exhibit vibro-impacts until a speed Ω = 2202 rpm where the fourth grazing bifurcation occurs, leading to another small interval without vibro-impacts. At Ω = 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. 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. [START_REF] Nevzat Özgüven | Mathematical models used in gear dynamics-a review[END_REF] at Ω = 2970 rpm, where the most energetic mesh mode is excited by the first harmonic of the internal excitation H Z . Fig. 15a and15b 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 µm which is much smaller than the gear backlash b = 40 µm.

Figure 15c shows the harmonic content of the dynamic mesh load at Ω = 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. [START_REF] Nevzat Özgüven | Mathematical models used in gear dynamics-a review[END_REF] at Ω = 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 Ω = 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 Ω = 6875 rpm. Compared to the results at Ω = 2970 rpm, harmonics H 2 and H 4 appears 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 (≈ 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 and17b shows two close-ups of the time series of the dynamic transmission error at Ω = 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 Ω = 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 Ω = 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: Fig. 17: Time series of the dynamic transmission error at Ω = 1700 rpm where the response exhibits one impact per period (a) and two impacts per period (b). The black dashed line correspond to the gap limit. -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 rotation), -2-IPP (5.25% of the rotation), -1-IPP (4% of the rotation).

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 largescale 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 vibroimpacts 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.

Fig. 1 :

 1 Fig. 1: Local frame of the gear pair.

Fig. 2 :

 2 Fig. 2: Example of aerodynamic forcing. The blue solid line and red solid line correspond to the forcing on the input and ouput shafts, respectively.

Fig. 3 :

 3 Fig. 3: Schematic representation of the system (a) and view of the FE model of the housing (b)

Fig. 4 :

 4 Fig. 4: Nonlinear force model.

Fig. 5 :

 5 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.

Fig. 6 :

 6 Fig. 6: Illustration of the AFT procedure.

Fig. 8 :

 8 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)

Fig. 9 :

 9 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.

Fig. 10 :Fig. 11 :

 1011 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)

Fig. 12 :Fig. 13 :

 1213 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.

Fig. 14 :

 14 Fig. 14: Forced response curve with both internal and external excitations. Blue diamond markers denote grazing bifurcations and blue highlights denote vibroimpacting responses.

Fig. 15 :Fig. 16 :

 1516 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 Ω = 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.

Fig. 18 :

 18 Fig.18: Time series of the dynamic transmission error at Ω = 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 correspond to the gap limit.

Table 1 :

 1 Characteristics of the gear pair

	Name		Gear 1 Gear 2 Unit
	Module	m	1		mm
	Number of teeth	Z	76	76	-
	Helix angle	β	0		deg
	Normal pressure angle	α	20		deg
	Base radius	r b	46.984	46.984	mm
	Profile shift coefficient	x	0	0	-
	Addendum coefficient	h a	1	1	-
	Dedendum coefficient	h d	1.25	1.25	-
	Face width	b f	20		mm
	Each shaft is supported by two ball bearings mod-
	elled by linear axial and radial stiffness elements. 3D
	quadratic tetrahedral elements are used to generate the

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).

Data availability

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

Conflict of interest

The authors declare that they have no conflict of interest.

Appendix