# Bifurcation tracking of geared systems with parameter-dependent internal excitation

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

► To cite this version:

Adrien Mélot, Emmanuel Rigaud, Joël Perret-Liaudet. Bifurcation tracking of geared systems with parameter-dependent internal excitation. Nonlinear Dynamics, 2021, 10.1007/s11071-021-07018-6. hal-03295433v2

HAL Id: hal-03295433

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

Submitted on 15 Jan 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Bifurcation tracking of geared systems with parameter-dependent internal excitation

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

Received: date / Accepted: date

**Abstract** 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.

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

## 1 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 structures [64,65], composite materials used in lightweight structures [22,68], friction in bolted assemblies [21,45], hertzian contact [47] and contact loss in geared systems [6]. 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 [59]. Time domain techniques are readily available in continuation toolboxes such as COCO [12], AUTO [16] (with specialized drivers SLIDECONT [13], HOMCONT [9] and TC-HAT [62]) or MATCONT [28]. 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 [1,2,51,70], rotordynamics [56], brake squeal [11], bladed disks dynamics [52], nonlinear electrical circuits [26,25], fluid-structure interactions [60,19], micro-electro-mechanical systems (MEMS) [36], nonlinear tuned vibrations absorbers [30,14]. Karkar et al. [40] 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 [67] provides a thorough comparison between two algorithms, HBM-PreCo and HBM-ANM. The former consists in performing a suc-

<footnote>A. Mélot · E. Rigaud · J. Perret-Liaudet
Laboratoire de Tribologie et Dynamique des Systèmes, UMR CNRS 5513, Ecole Centrale de Lyon, 36 avenue Guy de Collongue, 69134 Ecully Cedex, France
E-mail: adrien.melot@ec-lyon.fr</footnote>
# Bifurcation tracking of geared systems with parameter-dependent internal excitation

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

Received: date / Accepted: date

**Abstract** 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.

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

## 1 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 structures [64,65], composite materials used in lightweight structures [22,68], friction in bolted assemblies [21,45], hertzian contact [47] and contact loss in geared systems [6]. 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 [59]. Time domain techniques are readily available in continuation toolboxes such as COCO [12], AUTO [16] (with specialized drivers SLIDECONT [13], HOMCONT [9] and TC-HAT [62]) or MATCONT [28]. 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 [1, 2, 51, 70], rotordynamics [56], brake squeal [11], bladed disks dynamics [52], nonlinear electrical circuits [26,25], fluid-structure interactions [60,19], micro-electro-mechanical systems (MEMS) [36], nonlinear tuned vibrations absorbers [30,14]. Karkar et al. [40] 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 [67] provides a thorough comparison between two algorithms, HBM-PreCo and HBM-ANM. The former consists in performing a suc-

<footnote>A. Mélot · E. Rigaud · J. Perret-Liaudet
Laboratoire de Tribologie et Dynamique des Systèmes, UMR CNRS 5513, Ecole Centrale de Lyon, 36 avenue Guy de Collongue, 69134 Ecully Cedex, France
E-mail: adrien.melot@ec-lyon.fr</footnote>
cession of predictions, usually in the direction of the
tangent of the solution branch at a given iteration, fol-
lowed by corrections with a Newton-like solver. The
nonlinear forces are computed by the alternating fre-
quency/time procedure (AFT) [5] which takes advan-
tage 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 [10,33,39]. It is shown that,
while the HBM-ANM algorithm offers interesting per-
formances with smooth nonlinearities, the AFT-based
algorithm is more suitable for problems involving non-
smooth 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 dy-
namic behaviour of the system. Unfortunately, such
analyses often lead to intractable computations. Be-
sides, mechanical systems operating in the nonlinear
regime can possess isolated solutions, i.e, solutions de-
tached from the main solution branch. An elegant so-
lution to these problems consists in following the evo-
lution of bifurcation points, where qualitative changes
of the system's behaviour occur, when a parameter is
varied. Kuether et al. [44] tracked saddle-node bifur-
cations to unveil the relationship between modal inter-
actions in multi-degree-of-freedom systems and isolated
solutions. Furthermore, Alcorta [3] studied the creation
of subharmonic isolas in a single degree-of-freedom os-
cillator with asymmetric clearances by tracking period-
doubling bifurcations. To achieve this, equations char-
acterizing the bifurcations of interest are usually ap-
pended to the equilibrium equations in order to form
either a standard extended system [3, 69, 30] or a min-
imally extended system [4,27]. 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 [58]
proposed a standard extended system to study branch
point bifurcations. Moore and Spence [50] later used
such extended systems to compute saddle-node bifur-
cations. More recently, Xie et al. [69] coupled stan-
dard extended systems with the HBM to follow saddle-
node, branch points and Neimark-Sacker bifurcations
and eventually period-doubling bifurcations [3]. Mini-
mally extended systems, on the other hand, have seen
a much more limited use. They were first introduced by
Griewank and Reddien [31]. Such systems rely on defin-
ing auxiliary quantities that are determined by solving
a linear system of equations. In particular, when the
bifurcation is characterized by the singularity of a jaco-
bian matrix, the expansive computation of a determi-

nant can be bypassed by using a bordering technique
[17], resulting in a single scalar equation added to the
extended system.

The parameters considered in the literature are usu-
ally the damping ratio, forcing amplitude and nonlin-
earity level. One peculiarity of geared systems is the in-
ternal excitation generated by the meshing process and
associated to the transmission error defined as the dif-
ference between the instantaneous position of the driven
gear and the position it would occupy if gears were geo-
metrically perfect and infinitely rigid [66]. It mainly re-
sults from teeth deflections and profile deviations, both
manufacturing errors and intentional profile modifica-
tions [66,53,8]. The transmission error is also respon-
sible 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 inter-
dependent. For instance, tooth deflections and conse-
quently the transmission error and the mesh stiffness
fluctuations vary with the input torque. Furthermore,
a change in the mean value of the mesh stiffness mod-
ifies the modal basis of the underlying linearized gear
system under load. Consequently, modifying of the in-
put torque not only changes the spectral content of the
excitation but also the frequency and amplitude thresh-
old at which contact loss is likely to occur, leading to
vibro-impacts between the active and reverse flanks of
the gear teeth [55].

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 fol-
lows: we first introduce the complex harmonic balance
method (CHBM) and its coupling with minimally ex-
tended 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 pro-
posed methodology are presented and discussed in sec-
tion 4 and section 5 concludes the paper.

## 2 Numerical procedure

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

2.1 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-
cession of predictions, usually in the direction of the
tangent of the solution branch at a given iteration, fol-
lowed by corrections with a Newton-like solver. The
nonlinear forces are computed by the alternating fre-
quency/time procedure (AFT) [5] which takes advan-
tage 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 [10,33,39]. It is shown that,
while the HBM-ANM algorithm offers interesting per-
formances with smooth nonlinearities, the AFT-based
algorithm is more suitable for problems involving non-
smooth 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 dy-
namic behaviour of the system. Unfortunately, such
analyses often lead to intractable computations. Be-
sides, mechanical systems operating in the nonlinear
regime can possess isolated solutions, i.e, solutions de-
tached from the main solution branch. An elegant so-
lution to these problems consists in following the evo-
lution of bifurcation points, where qualitative changes
of the system's behaviour occur, when a parameter is
varied. Kuether et al. [44] tracked saddle-node bifur-
cations to unveil the relationship between modal inter-
actions in multi-degree-of-freedom systems and isolated
solutions. Furthermore, Alcorta [3] studied the creation
of subharmonic isolas in a single degree-of-freedom os-
cillator with asymmetric clearances by tracking period-
doubling bifurcations. To achieve this, equations char-
acterizing the bifurcations of interest are usually ap-
pended to the equilibrium equations in order to form
either a standard extended system [3, 69, 30] or a min-
ally extended system [4, 27]. 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 [58]
proposed a standard extended system to study branch
point bifurcations. Moore and Spence [50] later used
such extended systems to compute saddle-node bifur-
cations. More recently, Xie et al. [69] coupled stan-
dard extended systems with the HBM to follow saddle-
node, branch points and Neimark-Sacker bifurcations
and eventually period-doubling bifurcations [3]. Mini-
mally extended systems, on the other hand, have seen
a much more limited use. They were first introduced by
Griewank and Reddien [31]. Such systems rely on defin-
ing auxiliary quantities that are determined by solving
a linear system of equations. In particular, when the
bifurcation is characterized by the singularity of a jaco-
bian matrix, the expansive computation of a determi-

nant can be bypassed by using a bordering technique
[17], resulting in a single scalar equation added to the
extended system.

The parameters considered in the literature are usu-
ally the damping ratio, forcing amplitude and nonlin-
earity level. One peculiarity of geared systems is the in-
ternal excitation generated by the meshing process and
associated to the transmission error defined as the dif-
ference between the instantaneous position of the driven
gear and the position it would occupy if gears were geo-
metrically perfect and infinitely rigid [66]. It mainly re-
sults from teeth deflections and profile deviations, both
manufacturing errors and intentional profile modifica-
tions [66,53,8]. The transmission error is also respon-
sible 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 inter-
dependent. For instance, tooth deflections and conse-
quently the transmission error and the mesh stiffness
fluctuations vary with the input torque. Furthermore,
a change in the mean value of the mesh stiffness mod-
ifies the modal basis of the underlying linearized gear
system under load. Consequently, modifying of the in-
put torque not only changes the spectral content of the
excitation but also the frequency and amplitude thresh-
old at which contact loss is likely to occur, leading to
vibro-impacts between the active and reverse flanks of
the gear teeth [55].

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 fol-
lows: we first introduce the complex harmonic balance
method (CHBM) and its coupling with minimally ex-
tended 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 pro-
posed methodology are presented and discussed in sec-
tion 4 and section 5 concludes the paper.

## 2 Numerical procedure

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

### 2.1 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:

$$M\ddot{q} + C\dot{q} + Kq + f_{nl}(q) = f_{ex} \quad (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 omitted 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 $\Omega$ in equation (1) as truncated complex Fourier series of order H:

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

where $\tilde{q}$ contains the coefficients of the one-sided Fourier transform. The frequency domain differential operator $\nabla$ is used to obtain the expressions of the velocities and accelerations:

$$\nabla = \mathrm{diag}(0, i, 2i, \dots, Hi) \quad (3)$$

It follows:

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

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

where $\otimes$ 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:

$$\langle f, g \rangle = \frac{2}{T} \int_{0}^{T} f(t) \overline{g(t)} dt \quad (6)$$

This yields a residual which consists of a set of $n(H+1)$ nonlinear algebraic equations

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

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

$$Z(\Omega) = \Omega^2 \nabla^2 \otimes M + \Omega \nabla \otimes C + I_{H+1} \otimes K \quad (8)$$

## 2.2 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 [37,54]. The corresponding proposed contact modelling is described in detail in section 3.2.

The alternating frequency-time (AFT) procedure [5] 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.

<figure><img src="image_1.png" /><figcaption>Fig. 1: Illustration of the implemented AFT procedure.</figcaption></figure>

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
trix form:

$$M\ddot{q} + C\dot{q} + Kq + f_{nl}(q) = f_{ex} \quad (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 omitted 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 $\Omega$ in equation (1) as truncated complex Fourier series of order H:

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

where $\tilde{q}$ contains the coefficients of the one-sided Fourier transform. The frequency domain differential operator $\nabla$ is used to obtain the expressions of the velocities and accelerations:

$$\nabla = \operatorname{diag}(0, i, 2i, \dots, Hi) \quad (3)$$

It follows:

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

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

where $\otimes$ 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:

$$\langle f, g \rangle = \frac{2}{T} \int_{0}^{T} f(t) \overline{g(t)} dt \quad (6)$$

This yields a residual which consists of a set of $n(H+1)$ nonlinear algebraic equations

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

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

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

## 2.2 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 [37,54]. The corresponding proposed contact modelling is described in detail in section 3.2.

The alternating frequency-time (AFT) procedure [5] 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.

<figure><img src="image_2.png" /><figcaption>Fig. 1: Illustration of the implemented AFT procedure.</figcaption></figure>

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.

## 2.3 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 [59] 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 $\tilde{q}(s)$ and frequency $\Omega(s)$ also depend on s. In the following we suppose that iteration $k$ is a converged solution of equation (1). 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 $\Delta\tilde{Q} = (\Delta\tilde{q}, \Delta\Omega)^T$:

$$ \begin{bmatrix} \partial_{\tilde{q}} \mathbf{R} & \partial_{\Omega} \mathbf{R} \\ \Delta \tilde{q}^T & \Delta \Omega \end{bmatrix}_{(k)} \begin{pmatrix} \Delta \tilde{q} \\ \Delta \Omega \end{pmatrix} = \begin{pmatrix} 0 \\ 1 \end{pmatrix} \qquad (9) $$

where $\partial_{\tilde{q}}$ and $\partial_{\Omega}$ 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:

$$ \begin{pmatrix} \tilde{q} \\ \Omega \end{pmatrix}_{(p)} = \begin{pmatrix} \tilde{q} \\ \Omega \end{pmatrix}_{(k)} + \Delta s \begin{pmatrix} \Delta \tilde{q} \\ \Delta \Omega \end{pmatrix} \qquad (10) $$

where $\Delta s$ is the continuation step. Since the frequency $\Omega$ is treated as an unknown, the system composed of $n(H+1)$ equations and $n(H+1)+1$ unknowns is under-constrained. The following constraint equation is therefore appended to equation (7):

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

Starting from the predicted solution, Newton-like corrections bounded by a radius $\Delta 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:

$$ \begin{pmatrix} \tilde{q} \\ \Omega \end{pmatrix}_{(k+1)} = \begin{pmatrix} \tilde{q} \\ \Omega \end{pmatrix}_{(k)} - \begin{bmatrix} \partial_{\tilde{q}} \mathbf{R} & \partial_{\Omega} \mathbf{R} \\ \partial_{\tilde{q}} P & \partial_{\Omega} P \end{bmatrix}_{(k)}^{-1} \mathbf{R}_{ex}^{(k)} \qquad (12) $$

<figure><img src="image_2.png" /><figcaption>Fig. 2: Illustration of the arc-length continuation.</figcaption></figure>

where $\mathbf{R}_{ex}^{(k)} = (\mathbf{R}(\tilde{q}, \Omega), P(\tilde{q}, \Omega, s))^T$ is the extended residual.

Note that system (12) involves complex-valued Fourier coefficients. Solving (12) with a Newton-like method thus requires the extended residual function to be holomorphic [43], 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.

## 2.4 Detection of bifurcation points

Local bifurcations<fnref n="1" /> 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 [32]. It requires a proper post-processing to give accurate results as previously observed in [46]. A periodic solution $\mathbf{q}$ of eq. (1) is first perturbed by an exponential term:

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

This solution is expressed as truncated Fourier series of order $H$ and plugged into equation (1). A Galerkin procedure is then applied and a first order Taylor expansion is performed on the nonlinear terms. This yields a quadratic eigenvalue problem:

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

<footnote n="1">¹ Since global bifurcations are out of the scope of this paper, local bifurcations are henceforth simply referred to as bifurcations.</footnote>
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.

## 2.3 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 [59] 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 $\tilde{q}(s)$ and frequency $\Omega(s)$ also depend on s. In the following we suppose that iteration k is a converged solution of equation (1). 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 $\Delta\tilde{Q} = (\Delta\tilde{q}, \Delta\Omega)^T$:

$$ \begin{bmatrix} \partial_{\tilde{q}} \mathbf{R} & \partial_{\Omega} \mathbf{R} \\ \Delta \tilde{q}^T & \Delta \Omega \end{bmatrix}_{(k)} \begin{pmatrix} \Delta \tilde{q} \\ \Delta \Omega \end{pmatrix} = \begin{pmatrix} 0 \\ 1 \end{pmatrix} \qquad (9) $$

where $\partial_{\tilde{q}}$ and $\partial_{\Omega}$ 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:

$$ \begin{pmatrix} \tilde{q} \\ \Omega \end{pmatrix}_{(p)} = \begin{pmatrix} \tilde{q} \\ \Omega \end{pmatrix}_{(k)} + \Delta s \begin{pmatrix} \Delta \tilde{q} \\ \Delta \Omega \end{pmatrix} \qquad (10) $$

where $\Delta s$ is the continuation step. Since the frequency $\Omega$ is treated as an unknown, the system composed of $n(H+1)$ equations and $n(H+1)+1$ unknowns is under-constrained. The following constraint equation is therefore appended to equation (7):

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

Starting from the predicted solution, Newton-like corrections bounded by a radius $\Delta 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:

$$ \begin{pmatrix} \tilde{q} \\ \Omega \end{pmatrix}_{(k+1)} = \begin{pmatrix} \tilde{q} \\ \Omega \end{pmatrix}_{(k)} - \begin{bmatrix} \partial_{\tilde{q}} \mathbf{R} & \partial_{\Omega} \mathbf{R} \\ \partial_{\tilde{q}} P & \partial_{\Omega} P \end{bmatrix}_{(k)}^{-1} \mathbf{R}_{ex}^{(k)} \qquad (12) $$

<figure><img src="image_3.png" /><figcaption>Fig. 2: Illustration of the arc-length continuation.</figcaption></figure>

where $\mathbf{R}_{ex}^{(k)} = (\mathbf{R}(\tilde{q}, \Omega), P(\tilde{q}, \Omega, s))^T$ is the extended residual.

Note that system (12) involves complex-valued Fourier coefficients. Solving (12) with a Newton-like method thus requires the extended residual function to be holomorphic [43], 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.

## 2.4 Detection of bifurcation points

Local bifurcations<fnref n="1" /> 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 [32]. It requires a proper post-processing to give accurate results as previously observed in [46]. A periodic solution $\mathbf{q}$ of eq. (1) is first perturbed by an exponential term:

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

This solution is expressed as truncated Fourier series of order H and plugged into equation (1). A Galerkin procedure is then applied and a first order Taylor expansion is performed on the nonlinear terms. This yields a quadratic eigenvalue problem:

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

<footnote n="1">¹ Since global bifurcations are out of the scope of this paper, local bifurcations are henceforth simply referred to as bifurcations.</footnote>
with

$$ \begin{aligned} \tilde{\mathbf{M}} &= \mathbf{I}_{2H+1} \otimes \mathbf{M} \\ \tilde{\mathbf{C}} &= \nabla \otimes 2\mathbf{M} + \mathbf{I}_{2H+1} \otimes \mathbf{C} \end{aligned} \qquad (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 [63]:

$$ \left[ \begin{array}{cc} \tilde{\mathbf{C}} & \partial_{\mathbf{q}} \mathbf{R} \\ -\mathbf{I}_{n(2H+1)} & \mathbf{0} \end{array} \right] + \lambda \left[ \begin{array}{cc} \tilde{\mathbf{M}} & \mathbf{0} \\ \mathbf{0} & \mathbf{I}_{n(2H+1)} \end{array} \right] = 0 \qquad (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 [49] 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 $\phi$ whose zeros correspond to a particular bifurcation point are defined. Branch points

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

$$ \phi_{BP,SN} = \det(\partial_{\mathbf{q}} \mathbf{R}) = 0 \qquad (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:

$$ \phi_{SN} = \Delta\Omega \qquad (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 period-doubling (PD) bifurcations. Both bifurcations are characterized by a pair of complex conjugate eigenvalues crossing the imaginary axis. Their detection rely on the bialternate matrix product [29]. For $\mathbf{A}$, $\mathbf{B}$ two square matrices of size $n$, the bialternate product $\mathbf{A} \odot \mathbf{B}$ is defined as

$$ \mathbf{A} \odot \mathbf{B}_{(p,q),(r,s)} = \frac{1}{2} \left( \begin{array}{c|c} |a_{pr} \ a_{ps}| & |b_{pr} \ b_{ps}| \\ \hline |b_{qr} \ b_{qs}| & |a_{qr} \ a_{qs}| \end{array} \right) \qquad (19) $$

with $(1 \le q < p \le n)$ and $(1 \le s < r \le n)$. The product $2\mathbf{A} \odot \mathbf{I}_n$ is a diagonal matrix composed of the sums $\lambda_i + \lambda_j$ of the eigenvalues of matrix $\mathbf{A}$ [27] and is singular when a pair of complex conjugate eigenvalues crosses the imaginary axis. A suitable test function is formulated with the eigenvalue matrix $\mathbf{V}$ obtained with Hill's method:

$$ \phi_{NS,PD} = |2\mathbf{V} \odot \mathbf{I}| \qquad (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 $\Omega/2$. A NS bifurcation is detected otherwise. In the special case of eigenvalues that are real conjugates, a neutral saddle point is detected and can therefore be ignored.

<figure><img src="image_3.png" /><figcaption>Fig. 3: Bifurcation scenarii depending on the Imaginary part of the crossing eigenvalue.</figcaption></figure>
with

$$ \begin{aligned} \tilde{\mathbf{M}} &= \mathbf{I}_{2H+1} \otimes \mathbf{M} \\ \tilde{\mathbf{C}} &= \nabla \otimes 2\mathbf{M} + \mathbf{I}_{2H+1} \otimes \mathbf{C} \end{aligned} \qquad (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 [63]:

$$ \left[ \begin{array}{cc} \tilde{\mathbf{C}} & \partial_{\mathbf{q}} \mathbf{R} \\ -\mathbf{I}_{n(2H+1)} & \mathbf{0} \end{array} \right] + \lambda \left[ \begin{array}{cc} \tilde{\mathbf{M}} & \mathbf{0} \\ \mathbf{0} & \mathbf{I}_{n(2H+1)} \end{array} \right] = 0 \qquad (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 [49] 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 $\phi$ whose zeros correspond to a particular bifurcation point are defined. Branch points

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

$$ \phi_{BP,SN} = \det(\partial_{\mathbf{q}} \mathbf{R}) = 0 \qquad (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:

$$ \phi_{SN} = \Delta\Omega \qquad (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 period-doubling (PD) bifurcations. Both bifurcations are characterized by a pair of complex conjugate eigenvalues crossing the imaginary axis. Their detection rely on the bialternate matrix product [29]. For $\mathbf{A}$, $\mathbf{B}$ two square matrices of size $n$, the bialternate product $\mathbf{A} \odot \mathbf{B}$ is defined as

$$ \mathbf{A} \odot \mathbf{B}_{(p,q),(r,s)} = \frac{1}{2} \left( \begin{array}{c|c} |a_{pr} \ a_{ps}| & |b_{pr} \ b_{ps}| \\ \hline |b_{qr} \ b_{qs}| & |a_{qr} \ a_{qs}| \end{array} \right) \qquad (19) $$

with $(1 \le q < p \le n)$ and $(1 \le s < r \le n)$. The product $2\mathbf{A} \odot \mathbf{I}_n$ is a diagonal matrix composed of the sums $\lambda_i + \lambda_j$ of the eigenvalues of matrix $\mathbf{A}$ [27] and is singular when a pair of complex conjugate eigenvalues crosses the imaginary axis. A suitable test function is formulated with the eigenvalue matrix $\mathbf{V}$ obtained with Hill's method:

$$ \phi_{NS,PD} = |2\mathbf{V} \odot \mathbf{I}| \qquad (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 $\Omega/2$. A NS bifurcation is detected otherwise. In the special case of eigenvalues that are real conjugates, a neutral saddle point is detected and can therefore be ignored.

<figure><img src="image_4.png" /><figcaption>Fig. 3: Bifurcation scenarios depending on the Imaginary part of the crossing eigenvalue.</figcaption></figure>
## 2.5 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 $\mu$. The new vector of unknowns reads:

$$ \tilde{\mathbf{Q}} = (\tilde{\mathbf{q}}^T \quad \mu \quad \Omega)^T \quad (21) $$

An equation characterizing the bifurcation of interest is added to the residual to provide closure:

$$ \mathbf{R}_{\text{track}} = \begin{pmatrix} \mathbf{R}(\tilde{\mathbf{q}}, \mu, \Omega) \\ g(\tilde{\mathbf{q}}, \mu, \Omega) \\ P(\tilde{\mathbf{q}}, \mu, \Omega, s) \end{pmatrix} = \mathbf{0} \quad (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(\tilde{\mathbf{q}}, \Omega, s) = (\Delta\tilde{\mathbf{q}})^T (\Delta\tilde{\mathbf{q}}) + \Delta\mu^2 + \Delta\Omega^2 - \Delta s^2 = 0 \quad (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 [28] and later adapted to the HBM framework by Detroux et al. [15]. It consists in defining a scalar function $g$ that is evaluated by solving a linear system of equation defined as:

$$ \begin{bmatrix} \mathbf{A} & \mathbf{b} \\ \mathbf{d}^\dagger & 0 \end{bmatrix} \begin{pmatrix} \mathbf{w} \\ g \end{pmatrix} = \begin{pmatrix} \mathbf{0} \\ 1 \end{pmatrix} \quad (24) $$

where $\mathbf{A}$ denotes the jacobian matrix $\partial_{\tilde{\mathbf{q}}} \mathbf{R}$ when SN and BP bifurcations are tracked, or $2\mathbf{V} \odot \mathbf{I}_{2n}$ in the case of NS and PD bifurcations. $\mathbf{b}$ and $\mathbf{d}$ are bordering vectors that ought to ensure that system (24) is nonsingular and the superscript $\dagger$ 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 $\mathbf{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 $\mathcal{N}(\mathbf{A})$ of $\mathbf{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 $\mathbf{e}_j$ of $\mathbb{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):

$$ \begin{pmatrix} \tilde{\mathbf{q}} \\ \mu \\ \Omega \end{pmatrix}_{(k+1)} = \begin{pmatrix} \tilde{\mathbf{q}} \\ \mu \\ \Omega \end{pmatrix}_{(k)} - \begin{bmatrix} \partial_{\tilde{\mathbf{q}}} \mathbf{R} & \partial_{\mu} \mathbf{R} & \partial_{\Omega} \mathbf{R} \\ \partial_{\tilde{\mathbf{q}}} g & \partial_{\mu} g & \partial_{\Omega} g \\ \partial_{\tilde{\mathbf{q}}} P & \partial_{\mu} P & \partial_{\Omega} P \end{bmatrix}_{(k)}^{-1} \mathbf{R}_{\text{track}}^{(k)} \quad (25) $$

As in the previous section, the extended residual is split into real and imaginary parts to be able to use jacobian-based solvers. Future work may look into implementing jacobian-free Krylov-Newton solvers such as GMRES [42]. The derivatives of the residual and arc-length equation with respect to the unknown are calculated analytically. In particular, $\partial_{\tilde{\mathbf{q}}} \mathbf{R}$ and $\partial_{\Omega} \mathbf{R}$ are essentially the same as with a classical forced response analysis. The analytical derivations for the additional term $\partial_{\mu} \mathbf{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 [15]:

$$ \partial_{\tilde{\mathbf{q}}} g = -\mathbf{v}^{\dagger} \partial_{\tilde{\mathbf{q}}} \mathbf{A} \mathbf{w} \quad (26) $$

$$ \partial_{\mu} g = -\mathbf{v}^{\dagger} \partial_{\mu} \mathbf{A} \mathbf{w} \quad (27) $$

$$ \partial_{\Omega} g = -\mathbf{v}^{\dagger} \partial_{\Omega} \mathbf{A} \mathbf{w} \quad (28) $$

where $\mathbf{v}$ is computed by solving:

$$ \begin{bmatrix} \mathbf{A} & \mathbf{b} \\ \mathbf{d}^\dagger & 0 \end{bmatrix}^{\dagger} \begin{pmatrix} \mathbf{v} \\ h \end{pmatrix} = \begin{pmatrix} \mathbf{0} \\ 1 \end{pmatrix} \quad (29) $$

The terms $\partial_{\tilde{\mathbf{q}}} \mathbf{A}$, $\partial_{\mu} \mathbf{A}$, $\partial_{\Omega} \mathbf{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

## 3.1 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 $\theta_1$ and $\theta_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:

$$ \begin{cases} I_1 \ddot{\theta}_1 - T_1 + r_{b,1} f_{nl} = 0 \\ I_2 \ddot{\theta}_2 + T_2 - r_{b,2} f_{nl} = 0 \end{cases} \quad (30) $$
## 2.5 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 $\mu$. The new vector of unknowns reads:

$$ \tilde{\mathbf{Q}} = (\tilde{\mathbf{q}}^T \quad \mu \quad \Omega)^T \qquad (21) $$

An equation characterizing the bifurcation of interest is added to the residual to provide closure:

$$ \mathbf{R}_{\text{track}} = \begin{pmatrix} \mathbf{R}(\tilde{\mathbf{q}}, \mu, \Omega) \\ g(\tilde{\mathbf{q}}, \mu, \Omega) \\ P(\tilde{\mathbf{q}}, \mu, \Omega, s) \end{pmatrix} = \mathbf{0} \qquad (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(\tilde{\mathbf{q}}, \Omega, s) = (\Delta\tilde{\mathbf{q}})^T (\Delta\tilde{\mathbf{q}}) + \Delta\mu^2 + \Delta\Omega^2 - \Delta s^2 = 0 \qquad (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 [28] and later adapted to the HBM framework by Detroux et al. [15]. It consists in defining a scalar function $g$ that is evaluated by solving a linear system of equation defined as:

$$ \begin{bmatrix} \mathbf{A} & \mathbf{b} \\ \mathbf{d}^\dagger & 0 \end{bmatrix} \begin{pmatrix} \mathbf{w} \\ g \end{pmatrix} = \begin{pmatrix} \mathbf{0} \\ 1 \end{pmatrix} \qquad (24) $$

where $\mathbf{A}$ denotes the jacobian matrix $\partial_{\tilde{\mathbf{q}}} \mathbf{R}$ when SN and BP bifurcations are tracked, or $2\mathbf{V} \odot \mathbf{I}_{2n}$ in the case of NS and PD bifurcations. $\mathbf{b}$ and $\mathbf{d}$ are bordering vectors that ought to ensure that system (24) is nonsingular and the superscript $\dagger$ 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 $\mathbf{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 $\mathcal{N}(\mathbf{A})$ of $\mathbf{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 $\mathbf{e}_j$ of $\mathbb{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):

$$ \begin{pmatrix} \tilde{\mathbf{q}} \\ \mu \\ \Omega \end{pmatrix}_{(k+1)} = \begin{pmatrix} \tilde{\mathbf{q}} \\ \mu \\ \Omega \end{pmatrix}_{(k)} - \begin{bmatrix} \partial_{\tilde{\mathbf{q}}} \mathbf{R} & \partial_{\mu} \mathbf{R} & \partial_{\Omega} \mathbf{R} \\ \partial_{\tilde{\mathbf{q}}} g & \partial_{\mu} g & \partial_{\Omega} g \\ \partial_{\tilde{\mathbf{q}}} P & \partial_{\mu} P & \partial_{\Omega} P \end{bmatrix}_{(k)}^{-1} \mathbf{R}_{\text{track}}^{(k)} \qquad (25) $$

As in the previous section, the extended residual is split into real and imaginary parts to be able to use jacobian-based solvers. Future work may look into implementing jacobian-free Krylov-Newton solvers such as GMRES [42]. The derivatives of the residual and arc-length equation with respect to the unknown are calculated analytically. In particular, $\partial_{\tilde{\mathbf{q}}} \mathbf{R}$ and $\partial_{\Omega} \mathbf{R}$ are essentially the same as with a classical forced response analysis. The analytical derivations for the additional term $\partial_{\mu} \mathbf{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 [15]:

$$ \partial_{\tilde{\mathbf{q}}} g = -\mathbf{v}^{\dagger} \partial_{\tilde{\mathbf{q}}} \mathbf{A} \mathbf{w} \qquad (26) $$

$$ \partial_{\mu} g = -\mathbf{v}^{\dagger} \partial_{\mu} \mathbf{A} \mathbf{w} \qquad (27) $$

$$ \partial_{\Omega} g = -\mathbf{v}^{\dagger} \partial_{\Omega} \mathbf{A} \mathbf{w} \qquad (28) $$

where $\mathbf{v}$ is computed by solving:

$$ \begin{bmatrix} \mathbf{A} & \mathbf{b} \\ \mathbf{d}^\dagger & 0 \end{bmatrix}^{\dagger} \begin{pmatrix} \mathbf{v} \\ h \end{pmatrix} = \begin{pmatrix} \mathbf{0} \\ 1 \end{pmatrix} \qquad (29) $$

The terms $\partial_{\tilde{\mathbf{q}}} \mathbf{A}$, $\partial_{\mu} \mathbf{A}$, $\partial_{\Omega} \mathbf{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

## 3.1 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 $\theta_1$ and $\theta_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:

$$ \begin{cases} I_1 \ddot{\theta}_1 - T_1 + r_{b,1} f_{nl} = 0 \\ I_2 \ddot{\theta}_2 + T_2 - r_{b,2} f_{nl} = 0 \end{cases} \qquad (30) $$
<figure><img src="image_4.png" /><figcaption>Fig. 4: Flow chart of the bifurcation tracking algorithm</figcaption></figure>

where $T_1$ and $T_2$ are the input and output torques, respectively, $r_{b,k}$ denotes the base radius of gear $k$ and $f_{nl}$ is the nonlinear restoring force acting between the gear teeth. These equations model a semi-definite sys-

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:

$$ \frac{\dot{\theta}_2}{\dot{\theta}_1} = \frac{Z_1}{Z_2} = \frac{r_{b,1}}{r_{b,2}} \qquad (31) $$

In the following, we introduce the transmission error $q$, defined as the relative displacement at the gear mesh:

$$ q = r_{b,1}\theta_1 - r_{b,2}\theta_2 \qquad (32) $$

The equation of motion can therefore be written as:

$$ M_{eq}\ddot{q} + f_{nl}(q, t) = F_s \qquad (33) $$

with

$$ \left\{ \begin{array}{l} M_{eq} = \frac{I_1 I_2}{r_{b,1}^2 I_2 + r_{b,2}^2 I_1} \\ F_s = \frac{T_1}{r_{b,1}} = \frac{T_2}{r_{b,2}} \end{array} \right. \qquad (34) $$

<table><caption>Table 1: Characteristics of the gear pair</caption><thead><tr><th>Name</th><th></th><th>Gear 1</th><th>Gear 2</th><th>Unit</th></tr></thead><tbody><tr><td>Module</td><td>m</td><td>2</td><td></td><td>mm</td></tr><tr><td>Number of teeth</td><td>Z</td><td>50</td><td>50</td><td>-</td></tr><tr><td>Pressure angle</td><td>α</td><td>20</td><td></td><td>deg</td></tr><tr><td>Base radius</td><td>r<sub>b</sub></td><td>46.984</td><td>46.984</td><td>mm</td></tr><tr><td>Profile shift coefficient</td><td>x</td><td>0</td><td>0</td><td>-</td></tr><tr><td>Addendum coefficient</td><td>h<sub>a</sub></td><td>1</td><td>1</td><td>-</td></tr><tr><td>Dedendum coefficient</td><td>h<sub>d</sub></td><td>1.25</td><td>1.25</td><td>-</td></tr><tr><td>Face width</td><td>b<sub>f</sub></td><td>20</td><td></td><td>mm</td></tr><tr><td colspan="5">Tip relief modification</td></tr><tr><td>Length</td><td>l</td><td>1.75</td><td>1.75</td><td>mm</td></tr><tr><td>Amount</td><td>a</td><td>5</td><td>5</td><td>μm</td></tr></tbody></table>
<figure><img src="image_5.png" /><figcaption>Fig. 4: Flow chart of the bifurcation tracking algorithm</figcaption></figure>

where $T_1$ and $T_2$ are the input and output torques, respectively, $r_{b,k}$ denotes the base radius of gear $k$ and $f_{nl}$ is the nonlinear restoring force acting between the gear teeth. These equations model a semi-definite sys-

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:

$$ \frac{\dot{\theta}_2}{\dot{\theta}_1} = \frac{Z_1}{Z_2} = \frac{r_{b,1}}{r_{b,2}} \qquad (31) $$

In the following, we introduce the transmission error $q$, defined as the relative displacement at the gear mesh:

$$ q = r_{b,1}\theta_1 - r_{b,2}\theta_2 \qquad (32) $$

The equation of motion can therefore be written as:

$$ M_{eq}\ddot{q} + f_{nl}(q, t) = F_s \qquad (33) $$

with

$$ \left\{ \begin{array}{l} M_{eq} = \frac{I_1 I_2}{r_{b,1}^2 I_2 + r_{b,2}^2 I_1} \\ F_s = \frac{T_1}{r_{b,1}} = \frac{T_2}{r_{b,2}} \end{array} \right. \qquad (34) $$

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

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 [20,38] 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\xi\sqrt{k_m M_{eq}}$, with $\xi = 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}\ddot{q} + C\dot{q} + f_{nl}(q,t) = F_s \quad (35)$$

## 3.2 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 $\theta_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) - \frac{F_s}{k_m(t)} \quad (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))\mathcal{H}(q-j(t)) \\ + k_m(t)(q+j(t))\mathcal{H}(-q-j(t)) \quad (37)$$

where $\mathcal{H}$ is the Heaviside step function.

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

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 [20,61], square root functions [34,67], arctangent/hyperbolic tangent functions [71] or quintic splines [41] to cite a few. More recently, Margielewicz et al. [48] and Saunders et al. [57] proposed two studies that aimed at assessing the capacity of regularized force laws to capture various nonlinear phenomena including chaotic responses.
<figure><img src="image_6.png" /><figcaption>Fig. 5: Nonlinear gear model.</figcaption></figure>

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 [20,38] 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\xi\sqrt{k_m M_{eq}}$, with $\xi = 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}\ddot{q} + C\dot{q} + f_{nl}(q,t) = F_s \quad (35)$$

## 3.2 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 $\theta_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) - \frac{F_s}{k_m(t)} \quad (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))\mathcal{H}(q-j(t)) \\ + k_m(t)(q+j(t))\mathcal{H}(-q-j(t)) \quad (37)$$

where $\mathcal{H}$ is the Heaviside step function.

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

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 [20,61], square root functions [34,67], arctangent/hyperbolic tangent functions [71] or quintic splines [41] to cite a few. More recently, Margielewicz et al. [48] and Saunders et al. [57] 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 con-vergence, such functions generate other difficulties. First, they induce new sources of errors since the amount of smoothening depends on a user-defined pa-rameter that needs to be fine tuned and can result in either a nonzero force when contact is lost or in in-correct 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. Fi-nally, with appropriate preconditioning, e.g. by using a Jacobi preconditioner as suggested in [43], the HBM is found to offer quite satisfying convergence, thus elimi-nating the need for regularization.

The STE $q_s(t)$ and mesh stiffness $k_m(t)$ being peri-odic signals whose fundamental frequency corresponds to the mesh frequency $\omega_m = Z_1\Omega$, they can be written as truncated Fourier series:

$$ \left\{ \begin{aligned} q_s(t) &= \sum_{k=0}^{H_{q_s}} q_s^{c,k} \cos(k\Omega t) + q_s^{s,k} \sin(k\Omega t) \\ k_m(t) &= \sum_{k=0}^{H_{k_m}} k_m^{c,k} \cos(k\Omega t) + k_m^{s,k} \sin(k\Omega t) \end{aligned} \right. \qquad (38) $$

where $H_{q_s}$ and $H_{k_m}$ are the harmonic truncation order of the STE and the mesh stiffness, respectively. Without loss of generality, we assume $H_{q_s} = H_{k_m}$ in this study.

The STE and mesh stiffness are computed with an in-house code. The complete numerical procedure is detailed in [53]. The STE is evaluated at several successive angular positions $\theta_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 $\theta_1$ for a transmitted load $F_s$ reads:

$$ \left\{ \begin{aligned} \mathbf{H}(\theta_1) \cdot \mathbf{P}(\theta_1) &= q_s(\theta_1) \cdot \mathbf{1} - \mathbf{e}(\theta_1) \\ \mathbf{1}^T \cdot \mathbf{P}(\theta_1) &= F_s \end{aligned} \right. \qquad (39) $$

under the constraints

$$ \left\{ \begin{array}{l} -\mathbf{H}(\theta_1) \cdot \mathbf{P}(\theta_1) + q_s(\theta_1) \cdot \mathbf{1} \geq \mathbf{e}(\theta_1) \\ \mathbf{P}(\theta_1) \geq \mathbf{0} \end{array} \right. \qquad (40) $$

The contact lines are discretized and a compliance matrix $\mathbf{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 [24]. 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 $\mathbf{e}$ and $\mathbf{P}(\theta_1)$ is the distributed load on the contact lines.

As with the STE, the mesh stiffness is defined for every angular position $\theta_1$ of the driving gear. It is computed by differentiating the transmitted load $F_s$ with respect to the STE $q_s$:

$$ k_m(\theta_1) = \frac{\partial F_s}{\partial q_s(\theta_1)} \qquad (41) $$

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

<figure><img src="image_8.png" /><figcaption>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 (—)</figcaption></figure>

Figure 7 shows the time evolutions of the STE and mesh stiffness computed with the above-described procedure. Several computations were carried out with
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 [43], 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 $\omega_m = Z_1\Omega$, they can be written as truncated Fourier series:

$$ \left\{ \begin{aligned} q_s(t) &= \sum_{k=0}^{H_{q_s}} q_s^{c,k} \cos(k\Omega t) + q_s^{s,k} \sin(k\Omega t) \\ k_m(t) &= \sum_{k=0}^{H_{k_m}} k_m^{c,k} \cos(k\Omega t) + k_m^{s,k} \sin(k\Omega t) \end{aligned} \right. \qquad (38) $$

where $H_{q_s}$ and $H_{k_m}$ are the harmonic truncation order of the STE and the mesh stiffness, respectively. Without loss of generality, we assume $H_{q_s} = H_{k_m}$ in this study.

The STE and mesh stiffness are computed with an in-house code. The complete numerical procedure is detailed in [53]. The STE is evaluated at several successive angular positions $\theta_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 $\theta_1$ for a transmitted load $F_s$ reads:

$$ \left\{ \begin{aligned} \mathbf{H}(\theta_1) \cdot \mathbf{P}(\theta_1) &= q_s(\theta_1) \cdot \mathbf{1} - \mathbf{e}(\theta_1) \\ \mathbf{1}^T \cdot \mathbf{P}(\theta_1) &= F_s \end{aligned} \right. \qquad (39) $$

under the constraints

$$ \left\{ \begin{array}{l} -\mathbf{H}(\theta_1) \cdot \mathbf{P}(\theta_1) + q_s(\theta_1) \cdot \mathbf{1} \geq \mathbf{e}(\theta_1) \\ \mathbf{P}(\theta_1) \geq \mathbf{0} \end{array} \right. \qquad (40) $$

The contact lines are discretized and a compliance matrix $\mathbf{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 [24]. 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 $\mathbf{e}$ and $\mathbf{P}(\theta_1)$ is the distributed load on the contact lines.

As with the STE, the mesh stiffness is defined for every angular position $\theta_1$ of the driving gear. It is computed by differentiating the transmitted load $F_s$ with respect to the STE $q_s$:

$$ k_m(\theta_1) = \frac{\partial F_s}{\partial q_s(\theta_1)} \qquad (41) $$

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

<figure><img src="image_9.png" /><figcaption>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 (—)</figcaption></figure>

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 [18]. 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.

### 3.3 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 $\mu$.

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

$$ \left\{ \begin{aligned} q_s(\mu, t) &= \sum_{k=0}^{H_{q_s}} q_s^{c,k}(\mu) \cos(k\Omega t) + q_s^{s,k}(\mu) \sin(k\Omega t) \\ k_m(\mu, t) &= \sum_{k=0}^{H_{k_m}} k_m^{c,k}(\mu) \cos(k\Omega t) + k_m^{s,k}(\mu) \sin(k\Omega t) \end{aligned} \right. \quad (42) $$

During the bifurcation tracking procedure, coefficients $q_s^{c,k}$, $q_s^{s,k}$, $k_m^{c,k}$ and $k_m^{s,k}$ and their derivatives have to be evaluated at a priori unknown values of parameter $\mu$. A third order polynomial fit was found to provide satisfying approximations of the Fourier coefficients which can exhibit a non-monotonic evolution with parameter $\mu$. However, the mean value of the mesh stiffness $\tilde{k}_m^0(\mu)$ is modelled with a hyperbolic tangent function $\tilde{K}_m^0(\mu)$ of the form:

$$ \tilde{K}_m^0(\mu) = \kappa_1 \tanh(\kappa_2\mu) + \kappa_3 \quad (43) $$

where $\kappa_1$, $\kappa_2$, $\kappa_3$ and $\kappa_4$ are four coefficients determined by solving a nonlinear curve-fitting problem in the least-squares sense, that is:

Find $\kappa = (\kappa_1, \kappa_2, \kappa_3, \kappa_4) \in \mathbb{R}^4$ that solves

$$ \min_{\kappa} \left\| \kappa_1 \cdot \tanh(\kappa_2\mu + \kappa_3) + \kappa_4 - \tilde{k}_m^0(\mu) \right\|_2^2 \quad (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 non-physical model due to oscillations known as Runge's phenomenon. As can be seen in Fig. 8 (with parameter $\mu$ 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.

<figure><img src="image_9.png" /><figcaption>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.</figcaption></figure>
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 [18]. 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.

### 3.3 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 $\mu$.

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

$$ egin{cases} q_s(\mu, t) = \sum_{k=0}^{H_{q_s}} q_s^{c,k}(\mu) \cos(k\Omega t) + q_s^{s,k}(\mu) \sin(k\Omega t) \\ k_m(\mu, t) = \sum_{k=0}^{H_{k_m}} k_m^{c,k}(\mu) \cos(k\Omega t) + k_m^{s,k}(\mu) \sin(k\Omega t) \end{cases} \quad (42) $$

During the bifurcation tracking procedure, coefficients $q_s^{c,k}$, $q_s^{s,k}$, $k_m^{c,k}$ and $k_m^{s,k}$ and their derivatives have to be evaluated at a priori unknown values of parameter $\mu$. A third order polynomial fit was found to provide satisfying approximations of the Fourier coefficients which can exhibit a non-monotonic evolution with parameter $\mu$. However, the mean value of the mesh stiffness $\tilde{k}_m^0(\mu)$ is modelled with a hyperbolic tangent function $\tilde{K}_m^0(\mu)$ of the form:

$$ \tilde{K}_m^0(\mu) = \kappa_1 \tanh(\kappa_2\mu) + \kappa_3 \quad (43) $$

where $\kappa_1$, $\kappa_2$, $\kappa_3$ and $\kappa_4$ are four coefficients determined by solving a nonlinear curve-fitting problem in the least-squares sense, that is:

Find $\kappa = (\kappa_1, \kappa_2, \kappa_3, \kappa_4) \in \mathbb{R}^4$ that solves

$$ \min_{\kappa} \left\| \kappa_1 \cdot \tanh(\kappa_2\mu + \kappa_3) + \kappa_4 - \tilde{k}_m^0(\mu) \right\|_2^2 \quad (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 non-physical model due to oscillations known as Runge's phenomenon. As can be seen in Fig. 8 (with parameter $\mu$ 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.

<figure><img src="image_10.png" /><figcaption>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.</figcaption></figure>
## 4 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 ($\omega_m = Z_1\Omega$) of the dynamic response are retained. The vector of harmonics reads $\mathbf{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_{RMS}$ corresponding to the root mean square of the fluctuation of the dynamic transmission error is plotted, that is:

$$DTE_{RMS} = \sqrt{\sum_{k>0} Q_k} \qquad (45)$$

with $Q_k$ the amplitude of the harmonic corresponding to the k-th element of vector $\mathbf{H}$ excluding the 0-th harmonic.

### 4.1 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 = \omega_0/2\pi = 5730$ Hz. This mode is excited by the internal excitation for an operating rotational speed $\Omega = 720$ rad/s and results are displayed in Figure 9.

<figure><img src="image_10.png" /><figcaption>Fig. 9: Forced response of the gear model for an input torque $T = 20$ N·m.</figcaption></figure>

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 $\Omega = 505$ rad/s and $\Omega = 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 $\Omega = \omega_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
## 4 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 ($\omega_m = Z_1\Omega$) of the dynamic response are retained. The vector of harmonics reads $\mathbf{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_{RMS}$ corresponding to the root mean square of the fluctuation of the dynamic transmission error is plotted, that is:

$$DTE_{RMS} = \sqrt{\sum_{k>0} Q_k} \qquad (45)$$

with $Q_k$ the amplitude of the harmonic corresponding to the k-th element of vector $\mathbf{H}$ excluding the 0-th harmonic.

### 4.1 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 = \omega_0/2\pi = 5730$ Hz. This mode is excited by the internal excitation for an operating rotational speed $\Omega = 720$ rad/s and results are displayed in Figure 9.

<figure><img src="image_11.png" /><figcaption>Fig. 9: Forced response of the gear model for an input torque $T = 20$ N·m.</figcaption></figure>

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 $\Omega = 505$ rad/s and $\Omega = 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 $\Omega = \omega_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
<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /><figcaption>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.</figcaption></figure>

to $7 \cdot 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.

## 4.2 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-

<figure><img src="image_13.png" /><figcaption>Fig. 11: Spectral content of the response at $\Omega = 360$ rad/s.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 12: Spectral content of the response on the upper branch of primary resonance at $\Omega = 400$ rad/s showing negligible harmonics.</figcaption></figure>

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 $\Omega = 720$ rad/s with an input torque $T = 20$ N·m and $\Omega = 790$ rad/s with a torque $T = 120$ N·m. This behaviour is induced by
<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /><figcaption>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.</figcaption></figure>

to $7 \cdot 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.

## 4.2 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-

<figure><img src="image_14.png" /><figcaption>Fig. 11: Spectral content of the response at $\Omega = 360$ rad/s.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 12: Spectral content of the response on the upper branch of primary resonance at $\Omega = 400$ rad/s showing negligible harmonics.</figcaption></figure>

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 $\Omega = 720$ rad/s with an input torque $T = 20$ N·m and $\Omega = 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\xi\sqrt{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 [35]. 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.

### 4.3 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 [7,18].

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 $\Omega = 750$ rad/s and is associated to an unstable region spanning about $\Delta\Omega = 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 $\mathbb{R}^{n(2H+1)}$ were found to yield much more accurate results.
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\xi\sqrt{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 [35]. 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.

### 4.3 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 [7,18].

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 $\Omega = 750$ rad/s and is associated to an unstable region spanning about $\Delta\Omega = 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 $\mathbb{R}^{n(2H+1)}$ were found to yield much more accurate results.
<figure><img src="image_15.png" /><figcaption>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 (—).</figcaption></figure>

Another computation can be initialized at the other end of the parameter range ($a = 15 \ \mu 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 \ \mu 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 \ \mu 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_{RMS} \approx 3.5 \ \mu m$ and $a \approx 5.2 \ \mu m$) and than the one associated to the first unstable region ($DTE_{RMS} \approx 3 \ \mu m$ and $a \approx 2.4 \ \mu 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 \cdot m$, using gears modified with a tip relief $a \in [3 \ 5] \ \mu m$ yield acceptable results. This result can also be used during
<figure><img src="image_16.png" /><figcaption>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 (—).</figcaption></figure>

Another computation can be initialized at the other end of the parameter range ($a = 15 \ \mu 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 \ \mu 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 \ \mu 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_{RMS} \approx 3.5 \ \mu m$ and $a \approx 5.2 \ \mu m$) and than the one associated to the first unstable region ($DTE_{RMS} \approx 3 \ \mu m$ and $a \approx 2.4 \ \mu 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 \cdot m$, using gears modified with a tip relief $a \in [3 \ 5] \ \mu m$ yield acceptable results. This result can also be used during
<figure><img src="image_16.png" /><figcaption>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 (—).</figcaption></figure>

the design phase to estimate the required manufactur-
ing tolerances.

Finally, although the method is illustrated on a
SDOF gear model, the proposed methodology can be
employed on MDOF systems including the shafts, bear-
ings 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 [23] or
empirical modifications [72].

5 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 bal-
ance method coupled to an arc-length continuation al-
gorithm and bordering technique to form a minimally
extended system. The system is solved using a Newton-
like method requiring derivatives that can be formu-
lated (semi-)analytically. Future work may explore the
possibility of implementing a jacobian-free Newton-
Krylov solver to alleviate the cumbersome implementa-
tion of parameter-dependent partial derivatives or the
potentially expensive computation via finite differences.
<figure><img src="image_17.png" /><figcaption>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 (—).</figcaption></figure>

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 [23] or empirical modifications [72].

## 5 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 Newton-like 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 multi-degree-of-freedom gear transmissions, meaning that the optimization of tooth profile modifications of large-scale geared systems is now within reach.

**Acknowledgements** This work was performed within the framework of the LabCom LADAGE (LAboratoire de Dynamique des engrenAGEs), created by the LTDS and the Vibratec Company and operated by the French National Research Agency (ANR-14-LAB6-0003). It was also performed within the framework of the LABEX CeLyA (ANR-10-LABX-0060) of Université de Lyon, within the program "Investissements d'Avenir" (ANR-16-IDEX-0005) operated by the French National Research Agency (ANR).

## Conflict of interest

The authors declare that they have no conflict of interest.

## Data availability

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

## References

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

2. Al-shyyab, A., Kahraman, A.: Non-linear dynamic analysis of a multi-mesh gear train using multi-term harmonic balance method: sub-harmonic motions. Journal of Sound and Vibration **279**(1-2), 417-451 (2005). DOI https://doi.org/10.1016/j.jsv.2003.11.029

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

4. Allgower, E.L., Schwetlick, H.: A general view of minimally extended systems for simple bifurcation points. ZAMM - Journal of Applied Mathematics and Mechanics / Zeitschrift für Angewandte Mathematik und Mechanik **77**(2), 83-97 (1997). DOI https://doi.org/10.1002/zamm.19970770203

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

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

7. Carbonelli, A., Perret-Liaudet, J., Rigaud, E., Le Bot, A.: Particle Swarm Optimization as an Efficient Computational Method in order to Minimize Vibrations of Multimesh Gears Transmission. Advances in Acoustics and Vibration **2011**, 195642 (2011). DOI https://doi.org/10.1155/2011/195642. Publisher: Hindawi Publishing Corporation

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

9. Champneys, A.R., Kuznetsov, Y.A., Sandstede, B.: A numerical toolbox for homoclinic bifurcation analysis. Int. J. Bifurcation Chaos Appl.Sci. Eng. **6**(5), 867-887 (1996)

10. Cochelin, B., Vergez, C.: A high order purely frequency-based harmonic balance formulation for continuation of periodic solutions. Journal of Sound and Vibration **324**(1), 243-262 (2009). DOI https://doi.org/10.1016/j.jsv.2009.01.054

11. Coudeyras, N., Nacivet, S., Sinou, J.J.: Periodic and quasi-periodic solutions for multi-instabilities involved in brake squeal. Journal of Sound and Vibration **328**(4), 520-540 (2009). DOI https://doi.org/10.1016/j.jsv.2009.08.017

12. Dankowicz, H., Schilder, F.: An Extended Continuation Problem for Bifurcation Analysis in the Presence of Constraints. Journal of Computational and Nonlinear Dynamics **6**(031003) (2010). DOI https://doi.org/10.1115/1.4002684

13. Dercole, F., Kuznetsov, Y.A.: Slidecont: An auto97 driver for bifurcation analysis of filippov systems. ACM Trans. Math. Softw. **31**(1), 95-119 (2005). DOI https://doi.org/10.1145/1055531.1055536

14. Detroux, T., Habib, G., Masset, L., Kerschen, G.: Performance, robustness and sensitivity analysis of the nonlinear tuned vibration absorber. Mechanical Systems and Signal Processing **60-61**, 799-809 (2015). DOI https://doi.org/10.1016/j.ymssp.2015.01.035

15. Detroux, T., Renson, L., Masset, L., Kerschen, G.: The harmonic balance method for bifurcation analysis of large-scale nonlinear mechanical systems. Computer Methods in Applied Mechanics and Engineering **296**, 18-38 (2015). DOI https://doi.org/10.1016/j.cma.2015.07.017

16. Doedel, E.J., Champneys, A.R., Fairgrieve, T.F., Kuznetsov, Y.A., Sandstede, B., Wang, X.: Auto 97: Continuation and bifurcation software for ordinary differential equations (with homcont) (1997)

17. Doedel, E.J., Govaerts, W., Kuznetsov, Y.A.: Computation of periodic solution bifurcations in odes using bordered systems. SIAM Journal on Numerical Analysis
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 multi-degree-of-freedom gear transmissions, meaning that the optimization of tooth profile modifications of large-scale geared systems is now within reach.

**Acknowledgements** This work was performed within the framework of the LabCom LADAGE (LAboratoire de Dynamique des engrenAGEs), created by the LTDS and the Vibratec Company and operated by the French National Research Agency (ANR-14-LAB6-0003). It was also performed within the framework of the LABEX CeLyA (ANR-10-LABX-0060) of Université de Lyon, within the program "Investissements d'Avenir" (ANR-16-IDEX-0005) operated by the French National Research Agency (ANR).

## Conflict of interest

The authors declare that they have no conflict of interest.

## Data availability

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

## References

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

2. Al-shyyab, A., Kahraman, A.: Non-linear dynamic analysis of a multi-mesh gear train using multi-term harmonic balance method: sub-harmonic motions. Journal of Sound and Vibration **279**(1-2), 417-451 (2005). DOI https://doi.org/10.1016/j.jsv.2003.11.029

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

4. Allgower, E.L., Schwetlick, H.: A general view of minimally extended systems for simple bifurcation points. ZAMM - Journal of Applied Mathematics and Mechanics / Zeitschrift für Angewandte Mathematik und Mechanik **77**(2), 83-97 (1997). DOI https://doi.org/10.1002/zamm.19970770203

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

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

7. Carbonelli, A., Perret-Liaudet, J., Rigaud, E., Le Bot, A.: Particle Swarm Optimization as an Efficient Computational Method in order to Minimize Vibrations of Multimesh Gears Transmission. Advances in Acoustics and Vibration **2011**, 195642 (2011). DOI https://doi.org/10.1155/2011/195642. Publisher: Hindawi Publishing Corporation

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

9. Champneys, A.R., Kuznetsov, Y.A., Sandstede, B.: A numerical toolbox for homoclinic bifurcation analysis. Int. J. Bifurcation Chaos Appl.Sci. Eng. **6**(5), 867-887 (1996)

10. Cochelin, B., Vergez, C.: A high order purely frequency-based harmonic balance formulation for continuation of periodic solutions. Journal of Sound and Vibration **324**(1), 243-262 (2009). DOI https://doi.org/10.1016/j.jsv.2009.01.054

11. Coudeyras, N., Nacivet, S., Sinou, J.J.: Periodic and quasi-periodic solutions for multi-instabilities involved in brake squeal. Journal of Sound and Vibration **328**(4), 520-540 (2009). DOI https://doi.org/10.1016/j.jsv.2009.08.017

12. Dankowicz, H., Schilder, F.: An Extended Continuation Problem for Bifurcation Analysis in the Presence of Constraints. Journal of Computational and Nonlinear Dynamics **6**(031003) (2010). DOI https://doi.org/10.1115/1.4002684

13. Dercole, F., Kuznetsov, Y.A.: Slidecont: An auto97 driver for bifurcation analysis of filippov systems. ACM Trans. Math. Softw. **31**(1), 95-119 (2005). DOI https://doi.org/10.1145/1055531.1055536

14. Detroux, T., Habib, G., Masset, L., Kerschen, G.: Performance, robustness and sensitivity analysis of the nonlinear tuned vibration absorber. Mechanical Systems and Signal Processing **60-61**, 799-809 (2015). DOI https://doi.org/10.1016/j.ymssp.2015.01.035

15. Detroux, T., Renson, L., Masset, L., Kerschen, G.: The harmonic balance method for bifurcation analysis of large-scale nonlinear mechanical systems. Computer Methods in Applied Mechanics and Engineering **296**, 18-38 (2015). DOI https://doi.org/10.1016/j.cma.2015.07.017

16. Doedel, E.J., Champneys, A.R., Fairgrieve, T.F., Kuznetsov, Y.A., Sandstede, B., Wang, X.: Auto 97: Continuation and bifurcation software for ordinary differential equations (with homcont) (1997)

17. Doedel, E.J., Govaerts, W., Kuznetsov, Y.A.: Computation of periodic solution bifurcations in odes using bordered systems. SIAM Journal on Numerical Analysis
41(2), 401-435 (2003). DOI https://doi.org/10.1137/S0036142902400779

18. Driot, N., Rigaud, E., Sabot, J., Perret-Liaudet, J.: Allocation of Gear Tolerances to Minimize Gearbox Noise Variability. Acta Acustica united with Acustica 87, 67–76 (2001)

19. Ekici, K., Kielb, R.E., Hall, K.C.: The effect of aerodynamic asymmetries on turbomachinery flutter. Journal of Fluids and Structures 36, 1–17 (2013). DOI https://doi.org/10.1016/j.jfluidstructs.2012.08.009

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

21. Festjens, H., Chevallier, G., Dion, J.: Nonlinear model order reduction of jointed structures for dynamic analysis. Journal of Sound and Vibration 333(7), 2100–2113 (2014). DOI https://doi.org/10.1016/j.jsv.2013.11.039

22. Ganapathi, M., Patel, B., Makhecha, D.: Nonlinear dynamic analysis of thick composite/sandwich laminates using an accurate higher-order theory. Composites Part B: Engineering 35(4), 345–355 (2004). DOI https://doi.org/10.1016/S1359-8368(02)00075-6

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

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

25. Gilmore, R.J., Steer, M.B.: Nonlinear circuit analysis using the method of harmonic balance—a review of the art. ii. advanced concepts. International Journal of Microwave and Millimeter-Wave Computer-Aided Engineering 1(2), 159–180 (1991). DOI https://doi.org/10.1002/mmce.4570010205

26. Gilmore, R.J., Steer, M.B.: Nonlinear circuit analysis using the method of harmonic balance—a review of the art. part i. introductory concepts. International Journal of Microwave and Millimeter-Wave Computer-Aided Engineering 1(1), 22–37 (1991). DOI https://doi.org/10.1002/mmce.4570010104

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

28. Govaerts, W., Kuznetsov, Y., Dhooge, A.: Numerical continuation of bifurcations of limit cycles in matlab. SIAM Journal ON Scientific Computing 27(1), 231–252 (2005)

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

30. Grenat, C., Baguet, S., Lamarque, C.H., Dufour, R.: A multi-parametric recursive continuation method for nonlinear dynamical systems. Mechanical Systems and Signal Processing 127, 276–289 (2019). DOI https://doi.org/10.1016/j.ymssp.2019.03.011

31. Griewank, A., Reddien, G.W.: Characterization and computation of generalized turning points. SIAM Journal on Numerical Analysis 21(1), 176–185 (1984). DOI https://doi.org/10.1137/0721012

32. Groll, G.V., Ewins, D.: The harmonic balance method with arc-length continuation in rotor/stator contact problems. Journal of Sound and Vibration 241(2), 223 – 233 (2001). DOI https://doi.org/10.1006/jsvi.2000.3298

33. Guillot, L., Cochelin, B., Vergez, C.: A Taylor series-based continuation method for solutions of dynamical systems. Nonlinear Dynamics 98(4), 2827–2845 (2019). DOI https://doi.org/10.1007/s11071-019-04989-5

34. Hilali, Youssef, Braikat, Bouazza, Lahmam, Hassane, Damil, Noureddine: An implicit algorithm for the dynamic study of nonlinear vibration of spur gear system with backlash. Mechanics & Industry 19(3), 310 (2018). DOI https://doi.org/10.1051/meca/2017006

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

36. Kacem, N., Baguet, S., Hentz, S., Dufour, R.: Computational and quasi-analytical models for non-linear vibrations of resonant mems and nems sensors. International Journal of Non-Linear Mechanics 46(3), 532–542 (2011). DOI https://doi.org/10.1016/j.ijnonlinmec.2010.12.012

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

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

39. Karkar, S., Cochelin, B., Vergez, C.: A high-order, purely frequency based harmonic balance formulation for continuation of periodic solutions: The case of non-polynomial nonlinearities. Journal of Sound and Vibration 332(4), 968–977 (2013). DOI https://doi.org/10.1016/j.jsv.2012.09.033

40. Karkar, S., Cochelin, B., Vergez, C.: A comparative study of the harmonic balance method and the orthogonal collocation method on stiff nonlinear systems. Journal of Sound and Vibration 333(12), 2554–2567 (2014). DOI https://doi.org/10.1016/j.jsv.2014.01.019

41. Kim, T.C., Rook, T., Singh, R.: Effect of smoothening functions on the frequency response of an oscillator with clearance non-linearity. Journal of Sound and Vibration 263, 665–678 (2003). DOI https://doi.org/10.1016/S0022-460X(02)01469-4

42. Knoll, D., Keyes, D.: Jacobian-free newton–krylov methods: a survey of approaches and applications. Journal of Computational Physics 193(2), 357–397 (2004). DOI https://doi.org/10.1016/j.jcp.2003.08.010

43. Krack, M., Gross, J.: Harmonic Balance for Nonlinear Vibration Problems (2019). DOI https://doi.org/10.1007/978-3-030-14023-6

44. Kuether, R., Renson, L., Detroux, T., Grappasonni, C., Kerschen, G., Allen, M.: Nonlinear normal modes, modal interactions and isolated resonance curves. Journal of Sound and Vibration 351, 299–310 (2015). DOI https://doi.org/10.1016/j.jsv.2015.04.035

45. Lacayo, R., Pesaresi, L., Groß, J., Fochler, D., Armand, J., Salles, L., Schwingshackl, C., Allen, M., Brake, M.: Nonlinear modeling of structures with bolted joints: A comparison of two approaches based on a time-domain and frequency-domain solver. Mechanical Systems and Signal Processing 114, 413–438 (2019). DOI https://doi.org/10.1016/j.ymssp.2018.05.033

46. Lazarus, A., Thomas, O.: A harmonic-based method for computing the stability of periodic solutions of dynamical
41(2), 401-435 (2003). DOI https://doi.org/10.1137/S0036142902400779

18. Driot, N., Rigaud, E., Sabot, J., Perret-Liaudet, J.: Allocation of Gear Tolerances to Minimize Gearbox Noise Variability. Acta Acustica united with Acustica **87**, 67–76 (2001)

19. Ekici, K., Kielb, R.E., Hall, K.C.: The effect of aerodynamic asymmetries on turbomachinery flutter. Journal of Fluids and Structures **36**, 1–17 (2013). DOI https://doi.org/10.1016/j.jfluidstructs.2012.08.009

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

21. Festjens, H., Chevallier, G., Dion, J.: Nonlinear model order reduction of jointed structures for dynamic analysis. Journal of Sound and Vibration **333**(7), 2100–2113 (2014). DOI https://doi.org/10.1016/j.jsv.2013.11.039

22. Ganapathi, M., Patel, B., Makhecha, D.: Nonlinear dynamic analysis of thick composite/sandwich laminates using an accurate higher-order theory. Composites Part B: Engineering **35**(4), 345–355 (2004). DOI https://doi.org/10.1016/S1359-8368(02)00075-6

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

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

25. Gilmore, R.J., Steer, M.B.: Nonlinear circuit analysis using the method of harmonic balance—a review of the art. ii. advanced concepts. International Journal of Microwave and Millimeter-Wave Computer-Aided Engineering **1**(2), 159–180 (1991). DOI https://doi.org/10.1002/mmce.4570010205

26. Gilmore, R.J., Steer, M.B.: Nonlinear circuit analysis using the method of harmonic balance—a review of the art. part i. introductory concepts. International Journal of Microwave and Millimeter-Wave Computer-Aided Engineering **1**(1), 22–37 (1991). DOI https://doi.org/10.1002/mmce.4570010104

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

28. Govaerts, W., Kuznetsov, Y., Dhooge, A.: Numerical continuation of bifurcations of limit cycles in matlab. SIAM Journal ON Scientific Computing **27**(1), 231–252 (2005)

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

30. Grenat, C., Baguet, S., Lamarque, C.H., Dufour, R.: A multi-parametric recursive continuation method for nonlinear dynamical systems. Mechanical Systems and Signal Processing **127**, 276–289 (2019). DOI https://doi.org/10.1016/j.ymssp.2019.03.011

31. Griewank, A., Reddien, G.W.: Characterization and computation of generalized turning points. SIAM Journal on Numerical Analysis **21**(1), 176–185 (1984). DOI https://doi.org/10.1137/0721012

32. Groll, G.V., Ewins, D.: The harmonic balance method with arc-length continuation in rotor/stator contact problems. Journal of Sound and Vibration **241**(2), 223 – 233 (2001). DOI https://doi.org/10.1006/jsvi.2000.3298

33. Guillot, L., Cochelin, B., Vergez, C.: A Taylor series-based continuation method for solutions of dynamical systems. Nonlinear Dynamics **98**(4), 2827–2845 (2019). DOI https://doi.org/10.1007/s11071-019-04989-5

34. Hilali, Youssef, Braikat, Bouazza, Lahmam, Hassane, Damil, Noureddine: An implicit algorithm for the dynamic study of nonlinear vibration of spur gear system with backlash. Mechanics & Industry **19**(3), 310 (2018). DOI https://doi.org/10.1051/meca/2017006

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

36. Kacem, N., Baguet, S., Hentz, S., Dufour, R.: Computational and quasi-analytical models for non-linear vibrations of resonant mems and nems sensors. International Journal of Non-Linear Mechanics **46**(3), 532–542 (2011). DOI https://doi.org/10.1016/j.ijnonlinmec.2010.12.012

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

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

39. Karkar, S., Cochelin, B., Vergez, C.: A high-order, purely frequency based harmonic balance formulation for continuation of periodic solutions: The case of non-polynomial nonlinearities. Journal of Sound and Vibration **332**(4), 968–977 (2013). DOI https://doi.org/10.1016/j.jsv.2012.09.033

40. Karkar, S., Cochelin, B., Vergez, C.: A comparative study of the harmonic balance method and the orthogonal collocation method on stiff nonlinear systems. Journal of Sound and Vibration **333**(12), 2554–2567 (2014). DOI https://doi.org/10.1016/j.jsv.2014.01.019

41. Kim, T.C., Rook, T., Singh, R.: Effect of smoothening functions on the frequency response of an oscillator with clearance non-linearity. Journal of Sound and Vibration **263**, 665–678 (2003). DOI https://doi.org/10.1016/S0022-460X(02)01469-4

42. Knoll, D., Keyes, D.: Jacobian-free newton–krylov methods: a survey of approaches and applications. Journal of Computational Physics **193**(2), 357–397 (2004). DOI https://doi.org/10.1016/j.jcp.2003.08.010

43. Krack, M., Gross, J.: Harmonic Balance for Nonlinear Vibration Problems (2019). DOI https://doi.org/10.1007/978-3-030-14023-6

44. Kuether, R., Renson, L., Detroux, T., Grappasonni, C., Kerschen, G., Allen, M.: Nonlinear normal modes, modal interactions and isolated resonance curves. Journal of Sound and Vibration **351**, 299–310 (2015). DOI https://doi.org/10.1016/j.jsv.2015.04.035

45. Lacayo, R., Pesaresi, L., Groß, J., Fochler, D., Armand, J., Salles, L., Schwingshackl, C., Allen, M., Brake, M.: Nonlinear modeling of structures with bolted joints: A comparison of two approaches based on a time-domain and frequency-domain solver. Mechanical Systems and Signal Processing **114**, 413–438 (2019). DOI https://doi.org/10.1016/j.ymssp.2018.05.033

46. Lazarus, A., Thomas, O.: A harmonic-based method for computing the stability of periodic solutions of dynamical
systems. Comptes Rendus Mécanique **338**(9), 510–517 (2010). DOI https://doi.org/10.1016/j.crmecrme.2010.07.020

47. Ma, Q.L., Kahraman, A., Perret-Liaudet, J., Rigaud, E.: An Investigation of Steady-State Dynamic Response of a Sphere-Plane Contact Interface With Contact Loss. Journal of Applied Mechanics **74**(2), 249–255 (2006). DOI https://doi.org/10.1115/1.2190230

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

49. Moore, G.: Floquet theory as a computational tool. SIAM Journal on Numerical Analysis **42**(6), 2522–2568 (2005)

50. Moore, G., Spence, A.: The calculation of turning points of nonlinear equations. SIAM Journal on Numerical Analysis **17**(4), 567–576 (1980)

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

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

53. Rigaud, E., Barday, D.: Modelling and Analysis of Static Transmission Error. Effect of Wheel Body Deformation and Interactions between Adjacent Loaded Teeth. In: 4th World Congress on Gearing and Power Transmission, pp. 1961–1972. Paris, France (1999)

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

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

56. Salles, L., Staples, B., Hoffmann, N., Schwingshackl, C.: Continuation techniques for analysis of whole aeroengine dynamics with imperfect bifurcations and isolated solutions. Nonlinear Dynamics **86**(3), 1897–1911 (2016). DOI https://doi.org/10.1007/s11071-016-3003-y

57. Saunders, B.E., Vasconcellos, R., Kuether, R.J., Abdelkefi, A.: Insights on the continuous representations of piecewise-smooth nonlinear systems: limits of applicability and effectiveness. Nonlinear Dynamics (2021). DOI https://doi.org/10.1007/s11071-021-06436-w

58. Seydel, R.: Numerical computation of branch points in nonlinear equations. Numerische Mathematik **33**(3), 339–352 (1979). DOI https://doi.org/10.1007/BF01398649

59. Seydel, R.: Practical Bifurcation and Stability Analysis. Springer-Verlag New York (2010). DOI https://doi.org/10.1007/978-1-4419-1740-9

60. Shahrzad, P., Mahzoon, M.: Limit cycle flutter of airfoils in steady and unsteady flows. Journal of Sound and Vibration **256**(2), 213–225 (2002). DOI https://doi.org/10.1006/jsvi.2001.4113

61. Si-yu, C., Jin-yuan, T.: Study on a new nonlinear parametric excitation equation: Stability and bifurcation. Journal of Sound and Vibration **318**(4), 1109–1118 (2008). DOI https://doi.org/10.1016/j.jsv.2008.04.055

62. Thota, P., Dankowicz, H.: Tc-hat: A novel toolbox for the continuation of periodic trajectories in hybrid dynamical systems. SIAM J. Appl. Dyn. Syst. **7**, 1283–1322 (2008)

63. Tisseur, F., Meerbergen, K.: The quadratic eigenvalue problem. SIAM Review **43**(2), 235–286 (2001)

64. Vizzaccaro, A., Givois, A., Longobardi, P., Shen, Y., Deü, J.F., Salles, L., Touzé, C., Thomas, O.: Non-intrusive reduced order modelling for the dynamics of geometrically nonlinear flat structures using three-dimensional finite elements. Computational Mechanics **66**(6), 1293–1319 (2020). DOI https://doi.org/10.1007/s00466-020-01902-5

65. Vizzaccaro, A., Shen, Y., Salles, L., Blahoš, J., Touzé, C.: Direct computation of nonlinear mapping via normal form for reduced-order models of finite element nonlinear structures. Computer Methods in Applied Mechanics and Engineering **384**, 113957 (2021). DOI https://doi.org/10.1016/j.cma.2021.113957

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

67. Woiwode, L., Balaji, N.N., Kappauf, J., Tubita, F., GuilLOT, L., Vergez, C., Cochelin, B., Grolet, A., Krack, M.: Comparison of two algorithms for harmonic balance and path continuation. Mechanical Systems and Signal Processing **136**, 106503 (2020). DOI https://doi.org/10.1016/j.ymssp.2019.106503

68. Wu, Z., Li, H., Friswell, M.I.: Advanced nonlinear dynamic modelling of bi-stable composite plates. Composite Structures **201**, 582–596 (2018). DOI https://doi.org/10.1016/j.compstruct.2018.06.072

69. Xie, L., Baguet, S., Prabel, B., Dufour, R.: Bifurcation tracking by harmonic balance method for performance tuning of nonlinear dynamical systems. Mechanical Systems and Signal Processing **88**, 445–461 (2017). DOI https://doi.org/10.1016/j.ymssp.2016.09.037

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

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

72. Öztürk, V.Y., Cigeroglu, E., Özgüven, H.N.: Ideal tooth profile modifications for improving nonlinear dynamic response of planetary gear trains. Journal of Sound and Vibration **500**, 116007 (2021). DOI https://doi.org/10.1016/j.jsv.2021.116007
systems. Comptes Rendus Mécanique **338**(9), 510–517 (2010). DOI https://doi.org/10.1016/j.crmecrme.2010.07.020

47. Ma, Q.L., Kahraman, A., Perret-Liaudet, J., Rigaud, E.: An Investigation of Steady-State Dynamic Response of a Sphere-Plane Contact Interface With Contact Loss. Journal of Applied Mechanics **74**(2), 249–255 (2006). DOI https://doi.org/10.1115/1.2190230

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

49. Moore, G.: Floquet theory as a computational tool. SIAM Journal on Numerical Analysis **42**(6), 2522–2568 (2005)

50. Moore, G., Spence, A.: The calculation of turning points of nonlinear equations. SIAM Journal on Numerical Analysis **17**(4), 567–576 (1980)

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

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

53. Rigaud, E., Barday, D.: Modelling and Analysis of Static Transmission Error. Effect of Wheel Body Deformation and Interactions between Adjacent Loaded Teeth. In: 4th World Congress on Gearing and Power Transmission, pp. 1961–1972. Paris, France (1999)

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

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

56. Salles, L., Staples, B., Hoffmann, N., Schwingshackl, C.: Continuation techniques for analysis of whole aeroengine dynamics with imperfect bifurcations and isolated solutions. Nonlinear Dynamics **86**(3), 1897–1911 (2016). DOI https://doi.org/10.1007/s11071-016-3003-y

57. Saunders, B.E., Vasconcellos, R., Kuether, R.J., Abdelkefi, A.: Insights on the continuous representations of piecewise-smooth nonlinear systems: limits of applicability and effectiveness. Nonlinear Dynamics (2021). DOI https://doi.org/10.1007/s11071-021-06436-w

58. Seydel, R.: Numerical computation of branch points in nonlinear equations. Numerische Mathematik **33**(3), 339–352 (1979). DOI https://doi.org/10.1007/BF01398649

59. Seydel, R.: Practical Bifurcation and Stability Analysis. Springer-Verlag New York (2010). DOI https://doi.org/10.1007/978-1-4419-1740-9

60. Shahrzad, P., Mahzoon, M.: Limit cycle flutter of airfoils in steady and unsteady flows. Journal of Sound and Vibration **256**(2), 213–225 (2002). DOI https://doi.org/10.1006/jsvi.2001.4113

61. Si-yu, C., Jin-yuan, T.: Study on a new nonlinear parametric excitation equation: Stability and bifurcation. Journal of Sound and Vibration **318**(4), 1109–1118 (2008). DOI https://doi.org/10.1016/j.jsv.2008.04.055

62. Thota, P., Dankowicz, H.: Tc-hat: A novel toolbox for the continuation of periodic trajectories in hybrid dynamical systems. SIAM J. Appl. Dyn. Syst. **7**, 1283–1322 (2008)

63. Tisseur, F., Meerbergen, K.: The quadratic eigenvalue problem. SIAM Review **43**(2), 235–286 (2001)

64. Vizzaccaro, A., Givois, A., Longobardi, P., Shen, Y., Deü, J.F., Salles, L., Touzé, C., Thomas, O.: Non-intrusive reduced order modelling for the dynamics of geometrically nonlinear flat structures using three-dimensional finite elements. Computational Mechanics **66**(6), 1293–1319 (2020). DOI https://doi.org/10.1007/s00466-020-01902-5

65. Vizzaccaro, A., Shen, Y., Salles, L., Blahoš, J., Touzé, C.: Direct computation of nonlinear mapping via normal form for reduced-order models of finite element nonlinear structures. Computer Methods in Applied Mechanics and Engineering **384**, 113957 (2021). DOI https://doi.org/10.1016/j.cma.2021.113957

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

67. Woiwode, L., Balaji, N.N., Kappauf, J., Tubita, F., GuilLOT, L., Vergez, C., Cochelin, B., Grolet, A., Krack, M.: Comparison of two algorithms for harmonic balance and path continuation. Mechanical Systems and Signal Processing **136**, 106503 (2020). DOI https://doi.org/10.1016/j.ymssp.2019.106503

68. Wu, Z., Li, H., Friswell, M.I.: Advanced nonlinear dynamic modelling of bi-stable composite plates. Composite Structures **201**, 582–596 (2018). DOI https://doi.org/10.1016/j.compstruct.2018.06.072

69. Xie, L., Baguet, S., Prabel, B., Dufour, R.: Bifurcation tracking by harmonic balance method for performance tuning of nonlinear dynamical systems. Mechanical Systems and Signal Processing **88**, 445–461 (2017). DOI https://doi.org/10.1016/j.ymssp.2016.09.037

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

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

72. Öztürk, V.Y., Cigeroglu, E., Özgüven, H.N.: Ideal tooth profile modifications for improving nonlinear dynamic response of planetary gear trains. Journal of Sound and Vibration **500**, 116007 (2021). DOI https://doi.org/10.1016/j.jsv.2021.116007