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

Emmanuel Rigaud

Joël Perret-Liaudet

Bifurcation tracking of geared systems with parameterdependent internal excitation

Keywords: Bifurcation analysis, Harmonic balance method, Nonlinear gear dynamics, Numerical continuation, Vibro-impacts, Gear backlash

We herein propose an algorithm for tracking smooth bifurcations of nonlinear systems with interdependent parameters. The approach is based on a complex formulation of the well-known Harmonic Balance Method (HBM). Hill's method is used to assess the stability of the computed forced response curves and a minimally extended system is built to allow for the parametric continuation of the detected bifurcation points. The feasibility of coupling HBM-based minimally extended systems and arclength continuation algorithms is established and demonstrated. The method offers an efficient way of determining the stability regions of the system. The methodology is applied on a spur gear pair model including the backlash nonlinearity and subjected to transmission error and mesh stiffness fluctuation whose harmonic contents depend on several parameters that do not appear explicitly in the equations of motion.

Introduction

The design of modern mechanical systems warrants the consideration of nonlinearities and their impact on the dynamic behaviour of the system. This constitutes a challenge, as sources of nonlinearity are numerous. One can, for instance, cite distributed geometric nonlinearities induced by large deflections in slender struc-tures [START_REF] Vizzaccaro | Nonintrusive reduced order modelling for the dynamics of geometrically nonlinear flat structures using threedimensional finite elements[END_REF][START_REF] Vizzaccaro | Direct computation of nonlinear mapping via normal form for reduced-order models of finite element nonlinear structures[END_REF], composite materials used in lightweight structures [START_REF] Ganapathi | Nonlinear dynamic analysis of thick composite/sandwich laminates using an accurate higher-order theory[END_REF][START_REF] Wu | Advanced nonlinear dynamic modelling of bi-stable composite plates[END_REF], friction in bolted assemblies [START_REF] Festjens | Nonlinear model order reduction of jointed structures for dynamic analysis[END_REF][START_REF] Lacayo | Nonlinear modeling of structures with bolted joints: A comparison of two approaches based on a time-domain and frequency-domain solver[END_REF], hertzian contact [START_REF] Ma | An Investigation of Steady-State Dynamic Response of a Sphere-Plane Contact Interface With Contact Loss[END_REF] and contact loss in geared systems [START_REF] Cao | Nonlinear dynamics of a spur gear pair with force-dependent mesh stiffness[END_REF]. This often results in complicated dynamics as nonlinear behaviours heavily deviate from their linear counterparts. Among the various phenomena that nonlinear systems can exhibit, one can mention periodic and quasi-periodic responses, subharmonic oscillations, chaos, coexisting solutions and so forth.

A nonlinear dynamic analysis usually consists in following the evolution of periodic solutions with the excitation frequency. The most commonly encountered techniques couple a time-domain method, such as the shooting method, or a frequency domain method, such as the harmonic balance (HBM), with a continuation algorithm [START_REF] Seydel | Practical Bifurcation and Stability Analysis[END_REF]. Time domain techniques are readily available in continuation toolboxes such as COCO [START_REF] Dankowicz | An Extended Continuation Problem for Bifurcation Analysis in the Presence of Constraints[END_REF], AUTO [START_REF] Doedel | Auto 97: Continuation and bifurcation software for ordinary differential equations[END_REF] (with specialized drivers SLIDECONT [START_REF] Dercole | Slidecont: An auto97 driver for bifurcation analysis of filippov systems[END_REF], HOMCONT [START_REF] Champneys | A numerical toolbox for homoclinic bifurcation analysis[END_REF] and TC-HAT [START_REF] Thota | Tc-hat: A novel toolbox for the continuation of periodic trajectories in hybrid dynamical systems[END_REF]) or MATCONT [START_REF] Govaerts | Numerical continuation of bifurcations of limit cycles in matlab[END_REF]. These softwares are based on a Petrov-Galerkin method, namely the orthogonal collocation. The HBM has seen a widespread use in fields as diverse as gear dynamics [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] Mélot | Effect of gear topology discontinuities on the nonlinear dynamic response of a multi-degree-offreedom gear train[END_REF][START_REF] Yavuz | Nonlinear timevarying dynamic analysis of a spiral bevel geared system[END_REF], rotordynamics [START_REF] Salles | Continuation techniques for analysis of whole aeroengine dynamics with imperfect bifurcations and isolated solutions[END_REF], brake squeal [START_REF] Coudeyras | Periodic and quasi-periodic solutions for multi-instabilities involved in brake squeal[END_REF], bladed disks dynamics [START_REF] Quaegebeur | Energy transfer between nodal diameters of cyclic symmetric structures exhibiting polynomial nonlinearities: Cyclic condition and analysis[END_REF], nonlinear electrical circuits [START_REF] Gilmore | Nonlinear circuit analysis using the method of harmonic balance-a review of the art. part i. introductory concepts[END_REF][START_REF] Gilmore | Nonlinear circuit analysis using the method of harmonic balance-a review of the art. ii. advanced concepts[END_REF], fluid-structure interactions [START_REF] Shahrzad | Limit cycle flutter of airfoils in steady and unsteady flows[END_REF][START_REF] Ekici | The effect of aerodynamic asymmetries on turbomachinery flutter[END_REF], micro-electro-mechanical systems (MEMS) [START_REF] Kacem | Computational and quasi-analytical models for non-linear vibrations of resonant mems and nems sensors[END_REF], nonlinear tuned vibrations absorbers [START_REF] Grenat | A multi-parametric recursive continuation method for nonlinear dynamical systems[END_REF][START_REF] Detroux | Performance, robustness and sensitivity analysis of the nonlinear tuned vibration absorber[END_REF]. Karkar et al. [START_REF] Karkar | A comparative study of the harmonic balance method and the orthogonal collocation method on stiff nonlinear systems[END_REF] proposed a comparison between the orthogonal collocation and the HBM. They proved that the latter was less memory-demanding and offered better convergence, especially when applied to nonsmooth nonlinearities.

One of the difficulties of the harmonic balance method lies in expressing the nonlinear forces in the frequency domain. Woiwode [START_REF] Woiwode | Comparison of two algorithms for harmonic balance and path continuation[END_REF] provides a thorough comparison between two algorithms, HBM-PreCo and HBM-ANM. The former consists in performing a suc-cession of predictions, usually in the direction of the tangent of the solution branch at a given iteration, followed by corrections with a Newton-like solver. The nonlinear forces are computed by the alternating frequency/time procedure (AFT) [START_REF] Cameron | An Alternating Frequency/Time Domain Method for Calculating the Steady-State Response of Nonlinear Dynamic Systems[END_REF] which takes advantage of the existence of a closed-form expression of the nonlinear forces in the time domain. The latter, on the other hand, is a pure frequency domain formulation of the HBM which relies on expanding the solution curve as high-order Taylor series [START_REF] Cochelin | A high order purely frequencybased harmonic balance formulation for continuation of periodic solutions[END_REF][START_REF] Guillot | A Taylor seriesbased continuation method for solutions of dynamical systems[END_REF][START_REF] Karkar | A high-order, purely frequency based harmonic balance formulation for continuation of periodic solutions: The case of non-polynomial nonlinearities[END_REF]. It is shown that, while the HBM-ANM algorithm offers interesting performances with smooth nonlinearities, the AFT-based algorithm is more suitable for problems involving nonsmooth nonlinearities such as contact.

Due to the multiplicity of possible phenomena and parameters involved, parametric analyses are required to gain a thorough understanding of the nonlinear dynamic behaviour of the system. Unfortunately, such analyses often lead to intractable computations. Besides, mechanical systems operating in the nonlinear regime can possess isolated solutions, i.e, solutions detached from the main solution branch. An elegant solution to these problems consists in following the evolution of bifurcation points, where qualitative changes of the system's behaviour occur, when a parameter is varied. Kuether et al. [START_REF] Kuether | Nonlinear normal modes, modal interactions and isolated resonance curves[END_REF] tracked saddle-node bifurcations to unveil the relationship between modal interactions in multi-degree-of-freedom systems and isolated solutions. Furthermore, Alcorta [START_REF] Alcorta | Period doubling bifurcation analysis and isolated sub-harmonic resonances in an oscillator with asymmetric clearances[END_REF] studied the creation of subharmonic isolas in a single degree-of-freedom oscillator with asymmetric clearances by tracking perioddoubling bifurcations. To achieve this, equations characterizing the bifurcations of interest are usually appended to the equilibrium equations in order to form either a standard extended system [START_REF] Alcorta | Period doubling bifurcation analysis and isolated sub-harmonic resonances in an oscillator with asymmetric clearances[END_REF][START_REF] Xie | Bifurcation tracking by harmonic balance method for performance tuning of nonlinear dynamical systems[END_REF][START_REF] Grenat | A multi-parametric recursive continuation method for nonlinear dynamical systems[END_REF] or a minimally extended system [START_REF] Allgower | A general view of minimally extended systems for simple bifurcation points[END_REF][START_REF] Govaerts | Numerical bifurcation analysis for odes[END_REF]. The former technique is quite common in the bifurcation tracking literature. The first attempts at numerically following singular bifurcations date back to the late 1970s. Seydel [START_REF] Seydel | Numerical computation of branch points in nonlinear equations[END_REF] proposed a standard extended system to study branch point bifurcations. Moore and Spence [START_REF] Moore | The calculation of turning points of nonlinear equations[END_REF] later used such extended systems to compute saddle-node bifurcations. More recently, Xie et al. [START_REF] Xie | Bifurcation tracking by harmonic balance method for performance tuning of nonlinear dynamical systems[END_REF] coupled standard extended systems with the HBM to follow saddlenode, branch points and Neimark-Sacker bifurcations and eventually period-doubling bifurcations [START_REF] Alcorta | Period doubling bifurcation analysis and isolated sub-harmonic resonances in an oscillator with asymmetric clearances[END_REF]. Minimally extended systems, on the other hand, have seen a much more limited use. They were first introduced by Griewank and Reddien [START_REF] Griewank | Characterization and computation of generalized turning points[END_REF]. Such systems rely on defining auxiliary quantities that are determined by solving a linear system of equations. In particular, when the bifurcation is characterized by the singularity of a jacobian matrix, the expansive computation of a determi-nant can be bypassed by using a bordering technique [START_REF] Doedel | Computation of periodic solution bifurcations in odes using bordered systems[END_REF], resulting in a single scalar equation added to the extended system.

The parameters considered in the literature are usually the damping ratio, forcing amplitude and nonlinearity level. One peculiarity of geared systems is the internal excitation generated by the meshing process and associated to the transmission error defined as the difference between the instantaneous position of the driven gear and the position it would occupy if gears were geometrically perfect and infinitely rigid [START_REF] Welbourn | Fundamental knowledge of gear noise: A survey[END_REF]. It mainly results from teeth deflections and profile deviations, both manufacturing errors and intentional profile modifications [START_REF] Welbourn | Fundamental knowledge of gear noise: A survey[END_REF][START_REF] Rigaud | Modelling and Analysis of Static Transmission Error. Effect of Wheel Body Deformation and Interactions between Adjacent Loaded Teeth[END_REF][START_REF] Carbonelli | Vibro-Acoustic Analysis of Geared Systems-Predicting and Controlling the Whining Noise[END_REF]. The transmission error is also responsible for the fluctuation of the mesh stiffness between the mating gear teeth. This generates a multi-harmonic parametric internal excitation of the system. Numerous parameters governing the dynamic response are interdependent. For instance, tooth deflections and consequently the transmission error and the mesh stiffness fluctuations vary with the input torque. Furthermore, a change in the mean value of the mesh stiffness modifies the modal basis of the underlying linearized gear system under load. Consequently, modifying of the input torque not only changes the spectral content of the excitation but also the frequency and amplitude threshold at which contact loss is likely to occur, leading to vibro-impacts between the active and reverse flanks of the gear teeth [START_REF] Rigaud | Effect of Elasticity of Shafts, Bearings, Casing and Couplings on the Critical Rotational Speeds of a Gearbox[END_REF].

To the author's best knowledge, the present paper is the very first attempt at employing bifurcation tracking to study nonlinear gear dynamics. It is structured as follows: we first introduce the complex harmonic balance method (CHBM) and its coupling with minimally extended systems and arclength continuation algorithms in section 2. The numerical gear model and specifics of how the interdependant parameters are handled are detailed in section 3. The results obtained with the proposed methodology are presented and discussed in section 4 and section 5 concludes the paper.

Numerical procedure

The objective of this section is to propose an original numerical procedure which allows the tracking of bifurcation points.

Complex Harmonic Balance Method

In the general case, the equation of motion of a n degrees-of-freedom (DoF) model can be written in ma-trix form:

Mq + C q + Kq + f nl (q) = f ex (1)
where M, C, K are respectively the mass, damping and stiffness matrices and q is the vector of displacements of each DoF. f nl and f ex are vectors containing the nonlinear forces and external periodic forcing, respectively. Note that obvious time dependence is ommitted in the notations for clarity.

The following original complex formulation of the harmonic balance method (CHBM) relies on expanding the periodic solution q of fundamental frequency Ω in equation ( 1) as truncated complex Fourier series of order H:

q = Re ∞ k=0 qk e ikΩt ≈ Re H k=0 qk e ikΩt (2)
where q contains the coefficients of the one-sided Fourier transform. The frequency domain differential operator ∇ is used to obtain the expressions of the velocities and accelerations :

∇ = diag(0, i, 2i, . . . , Hi) (3)
It follows:

q = [TΩ∇ ⊗ I n ] q (4) q = TΩ 2 ∇ 2 ⊗ I n q (5)
where ⊗ is the Kronecker product and I n the identity matrix of size n. The complex exponential functions form an orthogonal Fourier basis T on which the solution is decomposed. The resulting expression is then plugged into eq. ( 1). The passage from the time domain to the frequency domain is done by eliminating the time variable via a Galerkin projection on the same Fourier basis T with the following hermitian inner product:

f, g = 2 T T 0 f (t)g(t)dt (6)
This yields a residual which consists of a set of n(H +1) nonlinear algebraic equations R(q, Ω) = Z(Ω)q + fnl (q) -fex = 0

where Z(Ω) is the dynamic stiffness matrix duplicated on each considered harmonic:

Z(Ω) = Ω 2 ∇ 2 ⊗ M + Ω∇ ⊗ C + I H+1 ⊗ K (8)

Treatment of the nonlinear forces

One of the main challenges of using the HBM consists in determining the Fourier coefficients of the nonlinear forces. The target application of this study is a geared system. Motion and force are transmitted by contacts between mating gear teeth. An internal excitation is generated by the meshing process or the cyclic disturbance of the circular motion of input shaft. Furthermore, backlash (clearance) is introduced through manufacturing tolerances in order to allow for assembly and operation. Under particular operating conditions and beyond an excitation threshold, contact loss between gear teeth may occur, the gear response is then characterized by impacts between the active and/or the reverse tooth flanks, leading to a highly nonlinear and discontinuous dynamic mesh force [START_REF] Kadmiri | Experimental and numerical analysis of automotive gearbox rattle noise[END_REF][START_REF] Rigaud | Investigation of gear rattle noise including visualization of vibro-impact regimes[END_REF]. The corresponding proposed contact modelling is described in detail in section 3.2.

The alternating frequency-time (AFT) procedure [START_REF] Cameron | An Alternating Frequency/Time Domain Method for Calculating the Steady-State Response of Nonlinear Dynamic Systems[END_REF] is implemented in order to evaluate the Fourier coefficients of the nonlinear forces. In practice the nonlinear response is sampled at specific time instants and back and forth transformations between the time and frequency domains are performed. The nonlinear force law being discontinuous, this can induce sampling errors as the fixed discretization does not ensure a good approximation of the location of the discontinuity. The approach adopted in this study is the classical AFT where the DFT is replaced by the fast Fourier transform (FFT). This allows for a higher sampling rate, thus reducing the discretization-induced errors, while keeping the computational effort reasonable. The AFT algorithm is illustrated in Fig. 1. The Fourier coefficients of the solution at a given iteration are converted to the time domain by an inverse fast Fourier transform (iFFT). This time domain solution is used to compute the evolution of the nonlinear forces and their derivatives over one fundamental period of vibration. The resulting signal is then transferred back into the frequency domain by a FFT.

Continuation of the solution with respect to the excitation frequency

At a given frequency, the forced response curve of a nonlinear system may exhibit multiple solutions. From a purely geometrical standpoint, this phenomenon is associated with turning points, where the curve folds and reverses direction with respect to the excitation frequency. The frequency parameter is therefore freed during the analysis to overcome this issue. The proposed methodology couples the arc-length continuation procedure [START_REF] Seydel | Practical Bifurcation and Stability Analysis[END_REF] to the HBM.

The fundamental idea of the arc-length continuation is to parameterize the sought solution curve by the curvilinear abscissa s. As a result, the displacement q(s) and frequency Ω(s) also depend on s. In the following we suppose that iteration k is a converged solution of equation [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]. When that is not the case, a Newton-like solver is used with an initial guess until convergence is reached. Given a converged solution at iteration k, the solution at iteration k + 1 is first estimated by computing the tangent vector ∆ Q = (∆q, ∆Ω)

T :

∂ qR ∂ Ω R ∆q T ∆Ω (k) ∆q ∆Ω = 0 1 (9)
where ∂ q and ∂ Ω represent the partial derivatives with respect to the vector of Fourier coefficient and frequency, respectively. The predicted solution, denoted by the subscript (p), is given by:

q Ω (p) = q Ω (k) + ∆s ∆q ∆Ω ( 10
)
where ∆s is the continuation step. Since the frequency Ω is treated as an unknown, the system composed of n(H +1) equations and n(H +1)+1 unkwowns is underconstrained. The following constraint equation is therefore appended to equation ( 7):

P (q, Ω, s) = (∆q) T (∆q) + ∆Ω 2 -∆s 2 = 0 (11)
Starting from the predicted solution, Newton-like corrections bounded by a radius ∆s of a hypersphere centered on the solution at iteration k (see Fig. 2) are performed in order to iteratively compute the converged solution to this augmented system of equations: where R (k) ex = (R(q, Ω), P (q, Ω, s))

q Ω (k+1) = q Ω (k) - ∂ qR ∂ Ω R ∂ qP ∂ Ω P -1 (k) R (k) ex (12)
T is the extended residual.

Note that system (12) involves complex-valued Fourier coefficients. Solving [START_REF] Dankowicz | An Extended Continuation Problem for Bifurcation Analysis in the Presence of Constraints[END_REF] with a Newton-like method thus requires the extended residual function to be holomorphic [START_REF] Krack | Harmonic Balance for Nonlinear Vibration Problems[END_REF], which is not usually the case. The complex extended residual is therefore split into real and imaginary parts during the correction steps to be able to resort to real arithmetics.

Detection of bifurcation points

Local bifurcations1 correspond to sudden changes of the system behaviour. They are associated with a change of local stability and possibly mark the appearance of a different sort of solution. To study such changes, a stability analysis is carried out using Hill's method [START_REF] Groll | The harmonic balance method with arc-length continuation in rotor/stator contact problems[END_REF]. It requires a proper post-processing to give accurate results as previously observed in [START_REF] Lazarus | A harmonic-based method for computing the stability of periodic solutions of dynamical systems[END_REF]. A periodic solution q of eq. ( 1) is first perturbed by an exponential term:

x 0 (t) = q(t) + p(t)e Λt ( 13
)
This solution is expressed as truncated Fourier series of order H and plugged into equation [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]. A Galerkin procedure is then applied and a first order Taylor expansion is performed on the nonlinear terms. This yields a quadratic eigenvalue problem :

Λ 2 M + Λ C + ∂ qR p = 0 (14)
with

M = I 2H+1 ⊗ M C = ∇ ⊗ 2M + I 2H+1 ⊗ C (15)
Note that the real harmonic balance formalism is used in the stability analysis since the complex form is solely used during the computation of the residual and the AFT algorithm. The above equation is then transformed into a linear eigenvalue problem to reduce the computational effort [START_REF] Tisseur | The quadratic eigenvalue problem[END_REF]:

C ∂ qR -I n(2H+1) 0 + λ M 0 0 I n(2H+1) = 0 (16)
Of the complete set of 2n(2H + 1) eigenvalues, only 2n are expected to have a physical signification, the remaining solutions being artifacts stemming from the multiplicity of harmonics. A robust post-processing of the computed eigenvalues is therefore necessary. Moore [START_REF] Moore | Floquet theory as a computational tool[END_REF] demonstrated that sorting the eigenvalues by their imaginary parts and keeping only the 2n eigenvalues with smallest imaginary part yielded accurate results for a sufficiently high harmonic truncation order. A converged solution point is eventually deemed stable if the real part of all retained eigenvalues is negative, and unstable otherwise. The bifurcation involved in the change of stability is deduced from the path of the eigenvalues as their real part switch from a negative value to a positive one (see Fig. 3). However, in the framework of bifurcation tracking, test functions φ whose zeros correspond to a particular bifurcation point are defined. Branch points Fig. 3: Bifurcation scenarii depending on the Imaginary part of the crossing eigenvalue.

(BP) and saddle-node (SN) bifurcations are characterized by the singularity of the jacobian matrix of the residual (7) that is:

φ BP,SN = det (∂ qR) = 0 (17)
Since both bifurcations are detected with this test function, a second criterion is used to discriminate SN from BP bifurcations. It is possible to differentiate SN and BP bifurcations by computing the rank of the jacobian matrix augmented with the vector corresponding to the derivatives with respect to the continuation parameter. However, when a SN bifurcation point is encountered, the solution branch folds and turns back with respect to this parameter. Such bifurcations are therefore identified by monitoring the sign of the frequency-component of the tangent vector at each iteration as it is faster to check for a change of sign of a scalar than to compute the rank of a matrix:

φ SN = ∆Ω (18)
A third test function is defined to detect Neimark-Sacker (NS) bifurcations, corresponding to the creation of a branch of quasi-periodic solutions, and perioddoubling (PD) bifurcations. Both bifurcations are characterized by a pair of complex conjugate eigenvalues crossing the imaginary axis. Their detection rely on the bialternate matrix product [START_REF] Govaerts | Matrix manifolds and the jordan structure of the bialternate matrix product[END_REF]. For A, B two square matrices of size n, the bialternate product A B is defined as

A B (p,q),(r,s) = 1 2 a pr a ps b qr b qs + b pr b ps a qr a qs (19)
with (1 ≤ q < p ≤ n) and (1 ≤ s < r ≤ n). The product 2A I n is a diagonal matrix composed of the sums λ i + λ j of the eigenvalues of matrix A [START_REF] Govaerts | Numerical bifurcation analysis for odes[END_REF] and is singular when a pair of complex conjugate eigenvalues crosses the imaginary axis. A suitable test function is formulated with the eigenvalue matrix V obtained with Hill's method:

φ N S,P D = 2V I (20)
In order to differentiate between NS and PD bifurcations, the imaginary part of the crossing eigenvalues is monitored. In case of a PD bifurcation, the imaginary part equals Ω/2. A NS bifurcation is detected otherwise. In the special case of eigenvalues that are real conjugates, a neutral saddle point is detected and can therefore be ignored.

Tracking of a bifurcation point with respect to an additional parameter

The system of equation is augmented in order to follow the evolution of a bifurcation with respect to parameter µ. The new vector of unknowns reads:

Q = qT µ Ω T ( 21
)
An equation characterizing the bifurcation of interest is added to the residual to provide closure:

R track =   R(q, µ, Ω) g(q, µ, Ω) P (q, µ, Ω, s)   = 0 ( 22
)
where P corresponds to the equation of a hypersphere, similar to that described in Sect. 2.3, in a n(2H + 1) + 2-dimensional space:

P (q, Ω, s) = (∆q) T (∆q) + ∆µ 2 + ∆Ω 2 -∆s 2 = 0 (23)
and g is a scalar function that characterizes the bifurcation of interest. As bifurcation points can be detected by monitoring the singularity of a matrix depending on the type of studied bifurcation, the proposed methodology relies on a bordering technique similar to that implemented in MATCONT [START_REF] Govaerts | Numerical continuation of bifurcations of limit cycles in matlab[END_REF] and later adapted to the HBM framework by Detroux et al. [START_REF] Detroux | The harmonic balance method for bifurcation analysis of large-scale nonlinear mechanical systems[END_REF]. It consists in defining a scalar function g that is evaluated by solving a linear system of equation defined as:

A b d † 0 w g = 0 1 ( 24
)
where A denotes the jacobian matrix ∂ qR when SN and BP bifurcations are tracked, or 2V I 2n in the case of NS and PD bifurcations. b and d are bordering vectors that ought to ensure that system (24) is nonsingular and the superscript † denotes the hermitian transpose. Special attention must be devoted to the choice of bordering vectors used during the tracking procedure as they greatly influence the quality of the results. For saddle-node bifurcations, vectors spanning the null space of matrix are commonly used. The singular value decomposition (SVD) of a matrix A is used to determine the span of its null space as the right singular vectors associated with singular values equal to zero provide an orthonormal basis of the null space N (A) of A. However, during the continuation and Newton corrections, the jacobian matrix may not be singular. In that case, the null vectors are approximated using the singular value decomposition (SVD) of the matrix.

Canonical basis vectors e j of R n(2H+1) , i.e. n(2H+1)-dimensional vectors whose j-th component is equal to one and all the others equal to zero, are also good candidates for bordering vectors depending on the application or the type of bifurcation being studied.

Corrections are then performed with a Newton-like solver (Fig. 4):

  q µ Ω   (k+1) =   q µ Ω   (k) -   ∂ qR ∂ µ R ∂ Ω R ∂ qg ∂ µ g ∂ Ω g ∂ qP ∂ µ P ∂ Ω P   -1 (k) R (k) track (25)
As in the previous section, the extended residual is split into real and imaginary parts to be able to use jacobianbased solvers. Future work may look into implementing jacobian-free Krylov-Newton solvers such as GM-RES [START_REF] Knoll | Jacobian-free newton-krylov methods: a survey of approaches and applications[END_REF]. The derivatives of the residual and arc-length equation with respect to the unknown are calculated analytically. In particular, ∂ qR and ∂ Ω R are essentially the same as with a classical forced response analysis.

The analytical derivations for the additional term ∂ µ R are detailed in the next section and depend on the considered dynamic model. The partial derivatives of the bordering equation with respect to the unknowns are expressed as [START_REF] Detroux | The harmonic balance method for bifurcation analysis of large-scale nonlinear mechanical systems[END_REF]:

∂ qg = -v † ∂ qAw (26)
∂ µ g = -v † ∂ µ Aw (27)
∂ Ω g = -v † ∂ Ω Aw ( 28
)
where v is computed by solving:

A b d † 0 † v h = 0 1 (29)
The terms ∂ qA, ∂ µ A, ∂ Ω A are computed by finite difference to avoid cumbersome analytical developments for every new tracking parameter. However, the computations are parallelized to speed up the code.

3 Model description

Mathematical model

The studied system corresponds to a reverse spur gear pair, with a gear ratio equal to 1:1. The characteristics of the gear pair are summarized in Table 1. The model assumes that the gear blanks are assimilated to 2 rigid disks with lumped inertia I 1 and I 2 . It has two degrees of freedom corresponding to the gear rotations θ 1 and θ 2 . They are connected by the nonlinear restoring force acting at the gear mesh, as depicted in Fig. 5. The equations of motion read as follows: tem with a rigid-body mode corresponding to the ideal input/output law of the gear pair in order to eliminate the rigid body mode:

I 1 θ1 -T 1 + r b,1 f nl = 0 I 2 θ2 + T 2 -r b,2 f nl = 0 (30)
θ2 θ1 = Z 1 Z 2 = r b,1 r b,2 (31)
In the following, we introduce the transmission error q, defined as the relative displacement at the gear mesh:

q = r b,1 θ 1 -r b,2 θ 2 (32)
The equation of motion can therefore be written as:

M eq q + f nl (q, t) = F s (33)
with where M eq is an equivalent mass associated with the lumped inertia I 1 and I 2 and F is the applied force associated with the input and output torques. The sources of damping in gear transmission are numerous and difficult to evaluate. It is therefore usual [START_REF] Farshidianfar | Global bifurcation and chaos analysis in nonlinear vibration of spur gear systems[END_REF][START_REF] Kahraman | Non-linear dynamics of a spur gear pair[END_REF] to define an equivalent viscous damping to model an average energy dissipation in the system, although this does not account changes in the state of contact. An equivalent viscous damping term C = 2ξ k m M eq , with ξ = 5% and k m the mean value of the mesh stiffness (sect. 3.2), is introduced at this stage. We consider a linear damping force that is not influenced by the contact loss nonlinearity. The complete equation of motion thus reads:

   M eq = I1I2 r 2 b,1 I2+r 2 b,2 I1 F s = T1 r b,1 = T2 r b,2 (34)
M eq q + C q + f nl (q, t) = F s (35)

Contact modelling

Depending on the operating conditions, the above described model can exhibit either:

a linear behaviour when the amplitude of the dynamic response is smaller than the static deflection, single-sided impacts when oscillations are sufficient to generate intermittent contact but small enough to never cross the backlash, double-sided impacts when the vibration amplitude is large enough that impacts on the reverse flank of the adjacent tooth occur.

The nonlinear force f nl representing the contact force between gear teeth, or lack thereof, is modelled by a piecewise (PW) linear function. The transmitted torque induces a static deflection q s , known as the static transmission error (STE), which depends on the angular position of the driving wheel θ 1 and is defined at a low rotational speed by eq. ( 32). This generates a periodic displacement excitation for non-zero rotational speeds. Linearizing the contact force around the static equilibrium leads to the definition of a periodic mesh stiffness k m (t) which in turn results in contact loss possibly occurring at a threshold j(t) (see Fig. 6):

j(t) = b + q s (t) - F s k m (t) ( 36
)
where b is the constant half backlash. The nonlinear mesh force can therefore be written as:

f nl (q, t) = k m (t) (q -j(t)) H (q -j(t)) + k m (t) (q + j(t)) H (-q -j(t)) (37)
where H is the Heaviside step function. One might argue that using a PW linear representation of the contact force might lead to numerical difficulties, especially convergence issues. In practice, PW smooth force laws are often regularized with third order polynomials [START_REF] Farshidianfar | Global bifurcation and chaos analysis in nonlinear vibration of spur gear systems[END_REF][START_REF] Si-Yu | Study on a new nonlinear parametric excitation equation: Stability and bifurcation[END_REF], square root functions [START_REF] Hilali | An implicit algorithm for the dynamic study of nonlinear vibration of spur gear system with backlash[END_REF][START_REF] Woiwode | Comparison of two algorithms for harmonic balance and path continuation[END_REF], arctangent/hyperbolic tangent functions [START_REF] Yoon | Effect and feasibility analysis of the smoothening functions for clearance-type nonlinearity in a practical driveline system[END_REF] or quintic splines [START_REF] Kim | Effect of smoothening functions on the frequency response of an oscillator with clearance non-linearity[END_REF] to cite a few. More recently, Margielewicz et al. [START_REF] Margielewicz | Modelling of the gear backlash[END_REF] and Saunders et al. [START_REF] Saunders | Insights on the continuous representations of piecewise-smooth nonlinear systems: limits of applicability and effectiveness[END_REF] proposed two studies that aimed at assessing the capacity of regularized force laws to capture various nonlinear phenomena including chaotic responses.

Although regularization techniques may help convergence, such functions generate other difficulties. First, they induce new sources of errors since the amount of smoothening depends on a user-defined parameter that needs to be fine tuned and can result in either a nonzero force when contact is lost or in incorrect predictions of the contact forces. Besides, when the contact is stiff, regularization techniques often prove suboptimal since small errors in the displacement may lead to severe inaccuracies in the computed force. Finally, with appropriate preconditioning, e.g. by using a Jacobi preconditioner as suggested in [START_REF] Krack | Harmonic Balance for Nonlinear Vibration Problems[END_REF], the HBM is found to offer quite satisfying convergence, thus eliminating the need for regularization.

The STE q s (t) and mesh stiffness k m (t) being periodic signals whose fundamental frequency corresponds to the mesh frequency ω m = Z 1 Ω, they can be written 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) k m (t) = H km k=0 k c,k m cos(kΩt) + k s,k m sin(kΩt) (38)
where H qs and H km are the harmonic truncation order of the STE and the mesh stiffness, respectively. Without loss of generality, we assume H qs = H km in this study.

The STE and mesh stiffness are computed with an in-house code. The complete numerical procedure is detailed in [START_REF] Rigaud | Modelling and Analysis of Static Transmission Error. Effect of Wheel Body Deformation and Interactions between Adjacent Loaded Teeth[END_REF]. The STE is evaluated at several successive angular positions θ 1 with contact established along contact lines that are evaluated prior to the computation. The location of these theoretical contact lines is determined by means of a kinematic analysis carried out for each angular position of the driving wheel. This allows for the estimation of the location of the contact point on the gear teeth. The equation describing the tooth contact at an angular position θ 1 for a transmitted load F s reads:

H (θ 1 ) .P (θ 1 ) = q s (θ 1 ) .1 -e (θ 1 ) 1 T .P (θ 1 ) = F s (39)
under the constraints -H (θ 1 ) .P (θ 1 ) + q s (θ 1 ) .1 ≥ e (θ 1 )

P (θ 1 ) ≥ 0 (40)
The contact lines are discretized and a compliance matrix H is built using an analytical tooth bending model based on the Reissner-Mindlin thick plate theory coupled to a Ritz-Galerkin approximation. The interested reader can find more details on the approach in [START_REF] Garambois | NVH robust optimization of gear macro and microgeometries using an efficient tooth contact model[END_REF]. This matrix links the displacement stemming from the global deformation of the teeth to the applied forces. A local deformation estimated by Hertz theory is introduced a posteriori. Potential manufacturing errors, tooth profile modifications are introduced in vector e and P (θ 1 ) is the distributed load on the contact lines.

As with the STE, the mesh stiffness is defined for every angular position θ 1 of the driving gear. It is computed by differentiating the transmitted load F s with respect to the STE q s : Figure 7 shows the time evolutions of the STE and mesh stiffness computed with the above-described procedure. Several computations were carried out with torques ranging from 20 N•m to 200 N•m and results are displayed over two mesh periods. One can see that the mean value of the STE increases with the input torque. Besides, the peak-to-peak value of the static transmission error decreases for torques up to about 120 N•m and increases for higher input torques. This minimum is linked to the length and amount of tip relief applied to the gear tooth profile (Table . 1) and corresponds to the optimal torque for which the excitation is minimum [START_REF] Driot | Allocation of Gear Tolerances to Minimize Gearbox Noise Variability[END_REF]. Figure 7 also shows that the mean value and the fluctuations of the mesh stiffness no longer evolve for torques higher than 80 N•m.

k m (θ 1 ) = ∂F s ∂q s (θ 1 ) (41)

Internal excitations

As stated in the introduction, the internal excitations q s (t) and k m (t) are quite sensitive to several operating and design parameters such as the input torque or the tooth profile modifications. The main limitation is that the STE q s (t) and the mesh stiffness k m (t) are usually used as input data in the dynamic simulations of a geared system. Their dependence on these parameters therefore does not appear explicitly in the equations of motion. This issue has to be addressed since the bifurcation tracking procedure requires the terms of the equations of motion and their derivatives be evaluated for each value of the tracking parameter µ.

In the algorithm used in this paper, the time evolutions of the internal excitation are reconstructed from their Fourier coefficients in the AFT procedure:

         q s (µ, t) = Hq s k=0 q c,k s (µ) cos(kΩt) + q s,k s (µ) sin(kΩt) k m (µ, t) = H km k=0 k c,k m (µ) cos(kΩt) + k s,k m (µ) sin(kΩt) (42)
During the bifurcation tracking procedure, coefficients q c,k s , q s,k s , k c,k m and k s,k m and their derivatives have to be evaluated at a priori unknown values of parameter µ. A third order polynomial fit was found to provide satisfying approximations of the Fourier coefficients which can exhibit a non-monotonic evolution with parameter µ. However, the mean value of the mesh stiffness k0 m (µ) is modelled with a hyperbolic tangent function K0 m (µ) of the form:

K0 m (µ) = κ 1 tanh(κ 2 µ) + κ 3 ( 43
)
where κ 1 κ 2 , κ 3 and κ 4 are four coefficients determined by solving a nonlinear curve-fitting problem in the least-squares sense, that is:

Find κ = (κ 1 , κ 2 , κ 3 , κ 4 ) ∈ R 4 that solves min κ κ 1 • tanh(κ 2 µ + κ 3 ) + κ 4 -k0 m (µ) 2 2 ( 44
)
Such a function is chosen for physical reasons. Indeed, the mean value of the mesh stiffness is expected to increase with the input torque. Using a third order polynomial for this harmonic can therefore result in a nonphysical model due to oscillations known as Runge's phenomenon. As can be seen in Fig. 8 (with paramenter µ corresponding to the applied torque), a hyperbolic tangent fit prevents such oscillations and leads to a more accurate model. Note that this issue could potentially also be circumvented with the use of spline interpolation.

Results and discussion

The following section aims at presenting the results obtained with the proposed methodology applied to the above-described gear model. The analysis is carried out for an operating rotational speed range of [0 8000] rpm.

In the following, we focus on the primary resonance of the SDOF gear pair. A single harmonic is therefore retained both for the STE and mesh stiffness functions. It should however be stressed that the proposed methodology can be easily extended to additional harmonics if secondary resonances are of interest. For all computations, the mean value as well as 7 harmonics of the mesh frequency (ω m = Z 1 Ω) of the dynamic response are retained. The vector of harmonics reads H = (0 50 100 150 200 250 300 350) T . Since the investigated model consists of a single DoF, the sampling rate for the AFT procedure was chosen relatively high and set to 2 12 . The solution for each computed point being a n(2H +1)+2-dimensional vector, we define here a scalar quantity which is used to plot results. In the following, the evolution of DTE RM S corresponding to the root mean square of the fluctuation of the dynamic transmission error is plotted, that is:

DTE RM S = k>0 Q k (45)
with Q k the amplitude of the harmonic corresponding to the k-th element of vector H excluding the 0-th harmonic.

Description of the dynamic response

We first analyse the effect of the input torque on the resulting dynamics. Increasing the input torque leads to a larger static deflection of the gear teeth. Contact loss is therefore less likely to occur, up to a point where only large torque fluctuations can make the teeth lose contact. As with all continuation procedures, the bifurcation tracking computation requires an initial guess. This initial guess is computed thanks to a preliminary forced response computation at one end of the torque interval. Since the initial guess must be close to a bifurcation point, the preliminary computation ought to be performed with a parameter for which the system's behaviour is expected to be nonlinear. For instance, we are interested in the dynamic behaviour of the gear pair in a torque range of [20 120] N•m. Since a high torque tends to prevent contact loss, the forced response computation is carried out at the lowest torque of the parameter interval T = 20 N•m. The corresponding mean value of the mesh stiffness is equal to 447 N•µm which results in a natural frequency of the underlying linear system equal to f 0 = ω 0 /2π = 5730 Hz. This mode is excited by the internal excitation for an operating rotational speed Ω = 720 rad/s and results are displayed in Figure 9. Figure 9 shows that the system exhibits a pronounced softening behaviour around the primary resonance. This is due to the fact that the response amplitude is such that the gear teeth do not stay in contact over a full vibration cycle. This leads to an intermittent contact state and a smaller apparent mesh stiffness. It should be noted that the contact loss is severe enough that two saddle-node bifurcations emerge on the main solution branch at Ω = 505 rad/s and Ω = 400 rad/s, respectively. This is of particular importance since this gives rise to a hysteretic behaviour, i.e. the well-known amplitude-jump phenomenon, that may induce a sudden appearance of large-amplitude vibrations of the geared system. Another important aspect is that the system only exhibits single-sided impacts, i.e, vibro-impacts only occur on the active tooth flank and the amplitude of oscillation is never large enough to cross the (half) backlash (see Fig. 10). If that were the case, the solution branch around the primary resonance would fold back towards the high frequencies, resulting in a hardening behaviour.

One can also note that, although only one harmonic is considered in the internal excitation, a barely noticeable resonance arises around Ω = ω 0 /100 = 360 rad/s due to the convolution between the harmonic of the STE and the mesh stiffness (Fig. 11). To ensure that the computed solution has converged, another computation is carried out considering harmonics orders up to 7 • Z and ten sidebands around each harmonic of the mesh frequency. Figure 12 shows the harmonic content of a point located on the upper branch of the primary resonance, where the harmonic distortion due to the nonlinearity is expected to be largest. It is clear that only harmonics multiple of the number of gear teeth contribute to the response, therefore validating the choice of retained harmonics. Besides, odd harmonics are negligible due to the asymmetric modelling of the gear backlash.

Bifurcation tracking with respect to the input torque

Either detected saddle-node bifurcation can be chosen as a starting point for the tracking procedure. Depend- ing on the discretization of the forced response curve and ensuing quality of the initial guess, more iterations than that needed for the remaining of the tracking may be required. As explained is section 2.5, the null vectors of the jacobian matrix of the HBM residual computed via a SVD were used as bordering vectors. Figure 13 displays the saddle-node tracking curve and forced response curves for a few different values of input torque. Note that these additional curves are not required for the analysis and only used in this paper for clarity and validation.

Figure 13 shows that the location of the primary resonance shifts towards the high frequencies as the input torque is increased. Indeed, the resonance is excited at a rotational speed of about Ω = 720 rad/s with an input torque T = 20 N•m and Ω = 790 rad/s with a torque T = 120 N•m. This behaviour is induced by the torque-dependence of the average mesh stiffness. Note that the frequency at which the primary resonance occurs is similar for input torques T = 95 N•m and T = 100 N•m. This is coherent with Fig. 8a which shows a linear increase of the average mesh stiffness followed by a plateau around 540 N/µm. From figure 13, it is interesting to note that the maximum amplitude increases with the input torque until it starts to dwindle beyond T = 45 N•m. Despite the increasing maximum amplitude, the extent of the region of instability, bounded by the solid red curve, does not grow both in the torque-frequency and amplitude-frequency planes. since the lower amplitude-limit of the SN curve increases much faster than the upper limit.

One can see that not only does the resonance peak straightens but the maximum amplitude also suddenly decreases for high input torques. This highlights the interdependence between the applied torque and the internal excitation: a larger static load not only results in a larger static deflection but also has a significant influence on the harmonic content of both the STE and mesh stiffness. Besides, modifying the mean value of the mesh stiffness alters the damping term C = 2ξ k m M eq (see Fig. 8a). We can thus distinguish two torque ranges:

-Up to about T = 80 N•m, the mean value of the mesh stiffness varies almost linearly with the input torque. The evolution of the unstable region is therefore governed by two phenomena. Firstly, the larger static deflection of the gear teeth induced by the higher static load increases the amplitude of oscillation required to lose contact. Secondly, since the damping term is expressed as a square root of the average mesh stiffness, the system is expected to exhibit more energy dissipation, leading to smaller vibration amplitudes. -For torques higher than T = 80 N•m, the mean value of the mesh stiffness shows a horizontal asymptote and the damping no longer increases. Physically speaking, this comes from the fact that contact is fully established and the global tooth deflection overcomes the local deflection of the tooth profile deviations. As a result, the global dynamic behaviour is now only governed by the load-induced deflection.

It is important to stress that tracking saddle-node bifurcations in a geared system does not give the boundaries of vibro-impact responses as such responses are bounded by grazing bifurcations, i.e. impacting periodic orbits with a zero normal velocity [START_REF] Ibrahim | Vibro-Impact Dynamics[END_REF]. Tracking SN bifurcations can however provide an approximation: the larger the unstable region, the larger the frequency interval exhibiting vibro-impacts.

Besides, bifurcation tracking offers a very convenient way of computing the stability regions of the sys-tem as the analysis only requires a few preliminary computations of the STE and mesh stiffness fluctuations to build the model, contrary to having to finely discretize the torque range and perform STE and mesh stiffness fluctuations computations for each considered parameter value. Future work could look into developing a fully coupled static and dynamic solver which would eliminate the need for any preliminary computation.

Bifurcation tracking with respect to a tooth profile modification

The following section aims at presenting the analysis of the sensibility of the nonlinear dynamic response to an implicit parameter which is of utmost importance in the design of gears. Tooth profile modifications are commonly employed to modify the internal excitation and diminish the fluctuations of the static transmission error for a target operating torque [START_REF] Carbonelli | Particle Swarm Optimization as an Efficient Computational Method in order to Minimize Vibrations of Multimesh Gears Transmission[END_REF][START_REF] Driot | Allocation of Gear Tolerances to Minimize Gearbox Noise Variability[END_REF].

However, mechanical systems often operate in a given torque range which obviously includes the torque for which the STE fluctuation is minimized. Such designs usually overlook the excitation stemming from the fluctuating mesh stiffness. Besides, the sensibility of the nonlinear dynamic response to variations in the profile modifications, either intentional or coming from manufacturing defects, is seldom studied.

In the following, tip relief modifications ranging from a = 0 µm, corresponding to the unmodified tooth profile, to a = 15 µm are investigated. A torque T = 100 N•m is applied to the system. The previous section showed that such a torque is at the limit of the region of instability for an amount of tip relief of a = 5 µm. This section investigates the robustness of this design with respect to a variability of the parameter.

As in the previous section, a preliminary forced response computation is carried out at one end of the parameter interval. In our case, the unmodified gear pair (a = 0 µm) is used as a starting point. As shown in Fig. 14, the primary resonance lies around Ω = 750 rad/s and is associated to a unstable region spanning about ∆Ω = 100 rad/s between two SN bifurcations. Computing the associated SN curve reveals that increasing the amount of tip relief reduces the distance between the two SN bifurcations and the size of the associated unstable region. Just before reaching an amount of tip relief a = 3 µm, the two SN bifurcation merge together and disappear, leading to a stable response. Note that in this application, using null vectors in the bordered system (24) proved unsatisfactory. Canonical basis vectors of R n(2H+1) were found to yield much more accurate results. Another computation can be initialized at the other end of the parameter range (a = 15 µm) to verify whether increasing the amount of tip relief does not result in the appearance of an additional unstable region. The RMS value of the dynamic transmission error at low frequency is much higher than with the unmodified tooth profile, indicating a larger static deflection of the gear teeth. Interestingly, the forced response curve for an amount of tip relief of a = 15 µm reveals a heavily nonlinear behaviour with a vibration amplitude about twice as high as the response with an unmodified tooth profile.

This indicates that there exist a subset of parameter values that allows for a stable response over the whole frequency range. Computing the SN tracking curve, shows that increasing the amount of tip relief a beyond 5 µm prove suboptimal as another region of instability, and associated softening behaviour is created. One can also note that the point of emergence of SN bifurcations lies at a higher amplitude (DTE RM S ≈ 3.5 µm and a ≈ 5.2 µm) and than the one associated to the first unstable region (DTE RM S ≈ 3 µm and a ≈ 2.4 µm), confirming that the static deflection is larger.

From these results, we can define a set of amount of tip relief that ensures a stable response over the whole range of rotational speed of the system. Considering the studied system and an input torque T = 100 N•m, using gears modified with a tip relief a ∈ [3 5] µm yield acceptable results. This result can also be used during the design phase to estimate the required manufacturing tolerances.

Finally, although the method is illustrated on a SDOF gear model, the proposed methodology can be employed on MDOF systems including the shafts, bearings and housing flexibility and even large scale gear transmissions (with appropriate model order reduction) to optimize the tooth profile modifications of systems potentially operating in a nonlinear regime associated to contact loss and vibro-impact responses. This is a particularly important result, as profile modifications are usually designed using static computations [START_REF] Garambois | Multiphysics coupling between periodic gear mesh excitation and input/output fluctuating torques: Application to a roots vacuum pump[END_REF] or empirical modifications [START_REF] Öztürk | Ideal tooth profile modifications for improving nonlinear dynamic response of planetary gear trains[END_REF].

Concluding remarks

A frequency-domain method for tracking bifurcations and determining the domain of stability of nonlinear systems was introduced. The proposed methodology is based on a complex formulation of the harmonic balance method coupled to an arc-length continuation algorithm and bordering technique to form a minimally extended system. The system is solved using a Newtonlike method requiring derivatives that can be formulated (semi-)analytically. Future work may explore the possibility of implementing a jacobian-free Newton-Krylov solver to alleviate the cumbersome implementation of parameter-dependent partial derivatives or the potentially expensive computation via finite differences.

Bifurcation tracking analysis was applied to a geared system for the very first time. The algorithm was adapted to the peculiarities of such systems, namely the parameter-dependant internal excitation and damping (through the mean value of the mesh stiffness). The methodology can be easily extended to an arbitrary number of harmonics and be coupled with other means of computing the STE and mesh stiffness, e.g. analytical of FE-based methods. The feasibility of developing a fully-coupled solver could also be investigated.

Non-intrusive model order reduction techniques can easily be used to study the domains of stability of multidegree-of-freedom gear transmissions, meaning that the optimization of tooth profile modifications of large-scale geared systems is now within reach.

Fig. 1 :

 1 Fig. 1: Illustration of the implemented AFT procedure.

Fig. 2 :

 2 Fig. 2: Illustration of the arc-length continuation.

Fig. 4 :

 4 Fig. 4: Flow chart of the bifurcation tracking algorithm

Fig. 5 :

 5 Fig. 5: Nonlinear gear model.

Fig. 6 :

 6 Fig. 6: Nonlinear force model.

Fig. 7 :

 7 Fig. 7: Time evolution of the static transmission error and mesh stiffness computed by solving Eq. (39) and reconstructed with 6 harmonics for torque values ranging from 20 N•m ( ) to 200 N•m ( )

Fig. 8 :

 8 Fig.8: Model of the first harmonics: mean value of the mesh stiffness (a), real part of the 1st harmonic of the mesh stiffness (c), imaginary part of the 1st harmonic of the mesh stiffness (e), mean value of the STE (b), real part of the 1st harmonic of the STE (d), imaginary part of the 1st harmonic of the STE (f). The blue solid line ( ) and the red dots ( ) correspond to the fitted model and the computations with the in-house code, respectively.

Fig. 9 :

 9 Fig. 9: Forced response of the gear model for an input torque T = 20 N•m.

Fig. 10 :

 10 Fig. 10: Dynamic transmission error on top of the primary resonance after the saddle-node bifurcation (a) and close-up (b). The DTE and the gap limit are plotted as a blue solid line ( ) and black dashed line ( ), respectively.

Fig. 11 :Fig. 12 :

 1112 Fig. 11: Spectral content of the response at Ω = 360 rad/s.

Fig. 13 :

 13 Fig. 13: Forced response curves and saddle-node bifurcation tracking curve. Coloured solid and dashed lines in the amplitude-frequency planes correspond to stable responses and unstable responses, respectively. Red circle markers ( ) denote saddle-node bifurcations. The saddle-node tracking curve is represented as a red solid line ( ).

Fig. 14 :

 14 Fig. 14: Forced response curves and saddle-node bifurcation tracking curve. Coloured solid and dashed lines in the amplitude-frequency planes correspond to stable responses and unstable responses, respectively. Red circle markers ( ) denote saddle-node bifurcations. The saddle-node tracking curve is represented as a red solid line ( ).

Table 1 :

 1 Characteristics of the gear pair

	Name		Gear 1 Gear 2 Unit
	Module	m	2		mm
	Number of teeth	Z	50	50	-
	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
	Tip relief modification
	Length	l	1.75	1.75	mm
	Amount	a	5	5	µm

Since global bifurcations are out of the scope of this paper, local bifurcations are henceforth simply referred to as bifurcations.

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.


Enciu Petre
email: petre.enciu@g2elab.grenoble-inp.fr

Gerbaud Laurent

Wurtz Frédéric

Automatic Differentiation Applied for Optimization of Dynamical Systems

Simulation is ubiquitous in many scientific areas. Applied for dynamic systems usually by employing differential equations, it gives the time evolution of system states. In order to solve such problems, numerical integration algorithms are often required. Automatic Differentiation (AD) is introduced as a powerful technique to compute derivatives of functions given in the form of computer programs in a high level programming language such as FORTRAN, C or C++. Such technique fits perfectly in combination with gradient based optimization algorithms, provided that the derivatives are valued with no truncation or cancellation error. This paper intends to use Automatic Differentiation employed for numerical integration schemes of dynamical systems simulating electromechanical actuators. Then, the resulting derivatives are used for sizing such devices by means of gradient based constrained optimization.

I. INTRODUCTION

Sizing by optimization is nowadays of major interest since it provides a fast and reliable way to achieve, with low manufacturing costs, desired performances for products lacking of optimality usually by means of minimizing a cost function. We are particularly interested by constrained gradient based optimization using Sequential Quadratic Programming (SQP) algorithms [START_REF] Gloud | A second derivative SQP method: local convergence[END_REF]. Such algorithms require accurately valued derivatives of the objective function. This may be the origin of serious problems provided that often such functions may result from complex numerical algorithms. We are particularly interested in this paper by those objective functions resulting from numerical integration of Initial Value Problems (IVPs) of Ordinary Differential Equations (ODEs) simulating the motion of an active body actuated by the electromagnetic force in the context of electromechanical actuators.

A good compromise in the optimization context is Automatic Differentiation (AD) that is a term applied for a technique able to compute derivatives of functions described by computer programs. That is, this paper only uses AD for sizing dynamical actuators by means of gradient based constrained optimization. In particular, AD will be applied using ADOL-C tool [START_REF] Griewank | Algorithm 755: ADOL-C: A Package for the Automatic Differentiation of Algorithms Written in C/C++[END_REF].

II. OPTIMIZATION PROBLEM

This paper considers the particular optimal problem of an IVP formulated as in (1) to ( 5):

(1)

where denotes the state with its associated initial values . In the paper, two formulations of the state system are intentionally specified. The formulation in (2) represents an autonomous system, meaning that the time variable does not appear in the differential equation, while the formulation in (3) refers to a non-autonomous system.

is the constraint design parameters set and denotes the control. The objective function depends on the reached final states and parameters. Equality (5) represents the simulation end criterion in Fig. 1, meaning that the simulation stops when a state reaches a prescribed final state,

. This implies the existence of the final time or response time, , depending implicitly on parameters. Note that makes part of design parameters.

ݔ ݔ ݐ ݐ integration stop

Fig. 1. Integration end criterion

The gradient based optimization algorithms applied for the optimization problem in (1) require the gradients of the objective function. These are evaluated like in [START_REF] Atienza | A Methodology for the Sizing and the Optimization of an Electromagnetic Release[END_REF] provided that the final states depend also on parameters: [START_REF] Atienza | A Methodology for the Sizing and the Optimization of an Electromagnetic Release[END_REF] Also, one may calculate the partial derivatives of the response time with respect to parameters set. So, the response time is carried out in optimization as a constrained parameter in addition to formulation in (1).

III. AUTOMATIC DIFFERENTIATION

Automatic Differentiation is introduced as a powerful technique that computes error-free derivatives, up to machine precision, of functions described as computer programs in high-level languages such FORTRAN or C/C++. In [START_REF] Barrio | Performance of Taylor Series method for ODEs/DAEs[END_REF] a rich list of tools implementing AD is provided. Therefore, an AD tool could be a library that instruments a user program in order to be differentiated. Such tools require minor modifications on the initial source and they are implemented in packages like ADOL-C that is subject to this paper. In general they are using the operator overloading capabilities of certain programming languages such C++ and FORTRAN95. In order to value the partial derivatives in (1) one may employ ADOL-C over a numerical scheme integrating the ODE system in (2) or (3). The paper is then subject to two numerical integration strategies. First, it applies AD over an adaptive step size Rung-Kutta (RK) scheme as in (7):

(7)

where denotes a slope estimation and is the integration step which depends on design parameters for an adaptive step size scheme.

Recent studies [START_REF] Walter | Automatic Differentiation of explicit Runge-Kutta methods for optimal control[END_REF] were carried out for differentiating such schemes. The difference in [START_REF] Walter | Automatic Differentiation of explicit Runge-Kutta methods for optimal control[END_REF] is that the response time is prescribed in advance at a fixed value. Our approach intends to make use of it as a constrained design parameter carried out further in optimization, so, its corresponding derivatives are to be valued as explained before.

Secondly, truncated Taylor Series (TS), as in [START_REF] Barrio | Performance of Taylor Series method for ODEs/DAEs[END_REF] are applied to advance the solution of the ODE system in ( 2) or (3) over a time interval as in ( 8):

(8) where denotes the order Taylor coefficient. Paper [START_REF] Barrio | Performance of Taylor Series method for ODEs/DAEs[END_REF] provides numerical solutions for adaptive step size schemes for ODE solvers using Taylor expansion. Interesting here is that AD is used for solving the dynamic system, provided that ADOL-C is capable to value high-order Taylor coefficients of the autonomous system in (2) supposing that is sufficiently smooth. In the non-autonomous case like in (3) a special version of the ADOL-C routine responsible for Taylor coefficients valuation is applied. The differentiation of such integration schemes is made by using special drivers implemented in ADOL-C.

The differentiation of a RK integration scheme in (7) tends to be slower since the slope usually is represented by a complex algorithm in the case of schemes up to second degree. Contrary, the differentiation of ( 8) is faster since it represents a sum expansion. However, here a cost of high-order derivatives computation should be paid up to second degree. In the full paper comparisons will be carried out for both schemes in terms of efficiency as also as helpful aspects regarding the AD of such numerical integration algorithms.

IV. OPTIMIZATION GOAL

The benchmark in [START_REF] Atienza | A Methodology for the Sizing and the Optimization of an Electromagnetic Release[END_REF] of the electromechanical actuator modeling a circuit breaker in Fig. 2 is proposed for sizing by gradient constrained optimization. When the switch is turned off, the vacuum force produced by the magnet equilibrate the spring force. The simulation starts when the switch turns on. The electromagnetic force created by the coil cancels partially the magnet force. Consequently, the plunger will move, starting from initial position, , toward the upper bound, .

c o i l ܨ ௦ ሬሬሬሬሬԦ ܨ ሬሬሬሬሬሬሬԦ ܴ ݇ ܨ ሬሬሬሬሬሬሬሬሬԦ ܥ ݐ ൌ Ͳ ܧ ݖ ݖ ௫ ܪ ௫ ܮ ௫
Fig. 2. Electromechanical actuator

The dynamical system of the proposed device combines both the equations of the electrical circuit feeding the coil and the movement equations. The states are:

The response time is found from the end criterion in [START_REF] Barrio | Performance of Taylor Series method for ODEs/DAEs[END_REF], that is satisfied when the mobile plunger is bounded at . A multi objective optimization problem raises form this particular case. These objectives are given in table I. The design parameters in [START_REF] Griewank | Algorithm 755: ADOL-C: A Package for the Automatic Differentiation of Algorithms Written in C/C++[END_REF] are represented by all geometrical parameters of the studied benchmark. The optimization results will be presented in the full paper.

V. CONCLUSION This paper presents a particular optimization problem on a benchmark dealing with state variables in Ordinary Differential Equations. Runge-Kutta and Taylor expansion integration schemes are used to approximate these states. Both schemes are differentiated by employing Automatic Differentiation in order to value the gradients needed by SQP algorithms.