L Xie

S Baguet

B Prabel

R Dufour

Bifurcation tracking by Harmonic Balance Method for performance tuning of nonlinear dynamical systems

Keywords: Harmonic balance method, Periodic solutions, Quasiperiodic solutions, Continuation, Bifurcation tracking, Limit point, Neimark-Sacker bifurcation, Floquet exponents, Nonlinear energy sink, Jeffcott rotor

The aim of this paper is to provide an efficient frequency-domain method for bifurcation analysis of nonlinear dynamical systems. The proposed method consists in directly tracking the bifurcation points when a system parameter such as the excitation or nonlinearity level is varied. To this end, a so-called extended system comprising the equation of motion and an additional equation characterizing the bifurcation of interest is solved by means of the Harmonic Balance Method coupled with an arc-length continuation technique. In particular, an original extended system for the detection and tracking of Neimark-Sacker (secondary Hopf) bifurcations is introduced. By applying the methodology to a nonlinear energy sink and to a rotor-stator rubbing system, it is shown that the bifurcation tracking can be used to efficiently compute the boundaries of stability and/or dynamical regimes, i.e., safe operating zones.

Introduction

Industrial requirements in terms of security, cost reduction and increased performance push designers, manufacturers and operators to create more and more advanced technological equipment in which nonlinearities are now common. In this context, understanding and controlling nonlinear effects due to contact, large deflections, links or components such as bearings or friction dampers is an important issue. Resulting nonlinear systems can exhibit complex dynamical behaviours with specific features such as multi-solutions for a single value of the system parameters, amplitude or frequency jumps, internal resonances, period-doubling, quasi-periodic or chaotic motions [START_REF] Nayfeh | Nonlinear oscillations[END_REF][START_REF] Guckenheimer | Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields[END_REF][START_REF] Thompson | Nonlinear Dynamics and Chaos[END_REF][START_REF] Strogatz | Nonlinear Dynamics and Chaos: With Applications to Physics[END_REF]. However, for a given system, the systematic study of all these phenomena and their possible occurrence is generally out of reach because of the large number of parameters to be considered and the limited available computational resources. An overall understanding of the system's dynamics can nevertheless be obtained through the computation of periodic solutions, forced response curves and associated bifurcations.

The literature comprises various numerical methods for the direct computation of periodic solutions which can be classified into two main categories, namely time domain and frequency domain approaches. The shooting method [START_REF] Kerschen | Computation of nonlinear normal modes through shooting and pseudo-arclength computation[END_REF] and orthogonal collocation [START_REF] De Boor | Collocation at gaussian points[END_REF] which rely on solving a nonlinear boundary value problem are two popular time domain approaches. Orthogonal collocation is implemented for instance in AUTO [START_REF] Doedel | AUTO-07P: Continuation and Bifurcation Software for Ordinary Differential Equations[END_REF] and MATCONT [START_REF] Dhooge | Matcont: A matlab package for numerical bifurcation analysis of odes[END_REF] softwares. In the frequency domain, the most commonly used method is certainly the harmonic balance method (HBM) which consists in approximating the unknown state variables by means of truncated Fourier series. Since nonlinearities cannot be directly computed in the frequency domain, the standard HBM is usually coupled with the alternating frequency-time (AFT) scheme [START_REF] Cameron | An alternating frequency/time domain method for calculating the steady-state response of nonlinear dynamic systems[END_REF] which computes the nonlinear terms in the time domain and subsequently their Fourier coefficients. The AFT scheme is very popular due to its easy implementation, its computational efficiency and its ability to handle almost any type of nonlinearities. Over the past decades, the HBM has been extended to quasi-periodic solutions [START_REF] Schilder | Fourier methods for quasi-periodic oscillations[END_REF][START_REF] Guskov | Multi-dimensional harmonic balance applied to rotor dynamics[END_REF][START_REF] Peletan | Quasi-periodic harmonic balance method for rubbing self-induced vibrations in rotor-stator dynamics[END_REF][START_REF] Zhou | A variable-coefficient harmonic balance method for the prediction of quasi-periodic response in nonlinear systems[END_REF] and many improvements have been proposed, such as adaptive schemes that improve the performance by selecting only the harmonics of interest [START_REF] Jaumouille | An adaptive harmonic balance method for predicting the nonlinear dynamic responses of mechanical systems -application to bolted structures[END_REF][START_REF] Grolet | On a new harmonic selection technique for harmonic balance method[END_REF] as well as methods to handle systems with many distinct states [START_REF] Krack | A high-order harmonic balance method for systems with distinct states[END_REF] and strong or non-smooth nonlinearities [START_REF] Nacivet | A dynamic lagrangian frequencytime method for the vibration of dry-friction-damped systems[END_REF][START_REF] Schreyer | A mixed shooting harmonic balance method for unilaterally constrained mechanical systems[END_REF]. For the computation of forced response curves, i.e. the following of periodic solutions when a control parameter is varied, the HBM is coupled with a continuation technique, e.g. the arclength continuation based on tangent prediction steps and orthogonal corrections [START_REF] Seydel | Practical Bifurcation and Stability Analysis[END_REF][START_REF]Numerical Continuation Methods for Dynamical Systems[END_REF] or the so-called asymptotic numerical method [START_REF] Cochelin | A high order purely frequency-based harmonic balance formulation for continuation of periodic solutions[END_REF].

In an engineering context, the local stability of periodic solutions is often computed when following the response curve since it distinguishes between solutions that may or not be experimentally observed. Several algorithms operating either in the time or frequency domain are available [START_REF] Von Groll | The harmonic balance method with arc-length continuation in rotor/stator contact problems[END_REF][START_REF] Peletan | A comparison of stability computational methods for periodic solution of nonlinear problems with application to rotordynamics[END_REF]. The detection of bifurcations points is more rarely performed. However, their computation is of prime interest. For instance, a limit point (also called fold bifurcation) indicates a change of stability and is responsible for amplitude jumps that can lead to significant and possibly dramatic changes in the system response. A Neimark-Sacker (secondary Hopf) bifurcation corresponds to a change of motion regime and indicates the transition from a periodic to a quasiperiodic motion.

Consequently, the parametric analysis of bifurcations can be used to understand the effects of nonlinear phenomena and to determine the boundaries of stability and/or dynamical regimes, i.e., safe operating zones. The resulting bifurcation map is an efficient tool for designers in order to identify the relevant parameters ruling the system's behaviour and to choose appropriate sets of parameters that lead to optimal runs. A simplified approach for this parametric analysis consists in calculating the whole response curves for several values of a chosen parameter, and collect all the detected bifurcations. However, this approach is very expensive and produces unnecessary results since only bifurcation points are of interest. A more efficient approach consists in detecting a starting bifurcation point for a fixed value of the parameter of interest, then in directly tracking the path of bifurcations while this parameter is varied.

Two approaches exist for the precise computation of bifurcation points. The first one is based on the use of so-called standard extended systems and consists in introducing one or more additional equations characterizing the bifurcation. The second approach relies on minimally extended systems and bordering techniques in which only one scalar function is added. The direct calculation of limit points of nonlinear equations depending on a parameter was first introduced by Seydel [START_REF] Seydel | Numerical computation of branch points in ordinary differential equations[END_REF][START_REF] Seydel | Numerical computation of branch points in nonlinear equations[END_REF], Moore and Spence [START_REF] Moore | The calculation of turning points of nonlinear equations[END_REF] using standard extended systems. Many authors also utilized this approach for the direct calculation of critical points for post-buckling finite element problems [START_REF] Wriggers | A general procedure for the direct computation of turning and bifurcation points[END_REF][START_REF] Cardona | Evaluation of simple bifurcation points and post-critical path in large finite rotation problems[END_REF][START_REF] Lopez | Post-critical analysis of structures with a nonlinear pre-buckling state in the presence of imperfections[END_REF]. It was recently combined with HBM by Petrov [START_REF] Petrov | Analysis of bifurcations in multiharmonic analysis of nonlinear forced vibrations of gas-turbine engine structures with friction and gaps[END_REF] for the detection of branch point bifurcations, where two branches of solutions intersect, and branch-switching along curves of periodic solutions. The direct calculation of limit points by means of minimally extended systems was introduced in [START_REF] Griewank | Characterization and computation of generalized turning points[END_REF] and subsequently used and improved by many authors [START_REF] Govaerts | Numerical Methods for Bifurcations of Dynamical Equilibria[END_REF][START_REF] Battini | Improved minimal augmentation procedure for the direct computation of critical points[END_REF]. The computation of Hopf bifurcations for dynamical systems by means of standard extended systems originates from the work of Jepson [START_REF] Jepson | Numerical Hopf bifurcation[END_REF]. Several variations and improvements have then been developed by Griewank and Reddien [START_REF] Griewank | The calculation of hopf points by a direct method[END_REF] or Roose et al. [START_REF] Roose | An algorithm for the computation of hopf bifurcation points in comparison with other methods[END_REF][START_REF] Roose | A direct method for the computation of hopf bifurcation points[END_REF] among others. This type of algorithm is frequently used in fluid mechanics to detect instabilities when the Reynolds number reaches critical values [START_REF] Brezillon | A numerical algorithm coupling a bifurcating indicator and a direct method for the computation of hopf bifurcation points in fluid mechanics[END_REF]. Such standard extended systems are implemented in AUTO [START_REF] Doedel | AUTO-07P: Continuation and Bifurcation Software for Ordinary Differential Equations[END_REF] and LOCA [START_REF] Salinger | Bifurcation tracking algorithms and software for large scale applications[END_REF] softwares. The computation of Hopf bifurcations by means of minimally extended systems is detailed in [START_REF] Govaerts | Numerical Methods for Bifurcations of Dynamical Equilibria[END_REF][40] [START_REF] Werner | Computation of hopf bifurcation with bordered matrices[END_REF]. These minimally extended systems are implemented in MATCONT [START_REF] Dhooge | Matcont: A matlab package for numerical bifurcation analysis of odes[END_REF] software. A comprehensive review of the methods suitable for detecting bifurcations can be found in [START_REF] Seydel | Practical Bifurcation and Stability Analysis[END_REF] while in [START_REF] Hussein | On efficient methods for detecting hopf bifurcation with applications to power system instability prediction[END_REF] authors focus on Hopf bifurcations.

The numerical continuation of bifurcation points is much less addressed in the literature. The continuation of paths of limit points of nonlinear equations having two parameters was first investigated by Jepson and Spence [START_REF] Jepson | Folds in solutions of two parameter systems and their calculation. part i[END_REF] with standard extended systems. In a mechanical context, it was later used for studying the sensitivity of critical buckling loads to imperfections [START_REF] Eriksson | Numerical analysis of complex instability behaviour using incremental-iterative strategies[END_REF][START_REF] Baguet | Stability of thin-shell structures and imperfection sensitivity analysis with the asymptotic numerical method[END_REF][START_REF] Rezaiee-Pajand | Stability boundaries of two-parameter non-linear elastic structures[END_REF]. In MATCONT, the continuation of codimension-1 bifurcations of dynamical systems is performed by means of minimally extended systems. In [START_REF] Detroux | The harmonic balance method for bifurcation analysis of large-scale nonlinear mechanical systems[END_REF], Detroux et al. combined this approach with the HBM for the tracking of limit point, branch point and Neimark-Sacker bifurcations of large-scale mechanical systems. In this paper, we combine HBM and standard extended systems. We already used this approach in [START_REF] Xie | Numerical tracking of limit points for direct parametric analysis in nonlinear rotordynamics[END_REF] in the case of limit points. Here, this work is extended to all types of codimension-1 bifurcations. In particular, we build on the work of Griewank and Reddien [START_REF] Griewank | The calculation of hopf points by a direct method[END_REF] in order to propose an efficient algorithm for the computation and the tracking of Neimark-Sacker bifurcations.

The paper is organized as follows. The formulation of the harmonic balance method for the continuation of periodic solutions is presented in Section 2. The stability analysis and the characterization of the bifurcations are based on the Floquet exponents obtained from a quadratic eigenvalue problem as described in Section 3. The extended systems used for the computation of the bifurcations are then detailed in Section 4, with emphasis on computational issues such as the efficient calculation of the derivatives involved in the Newton-Raphson iterations. The direct tracking of these bifurcations, i.e., the continuation of bifurcation curves is addressed in Section 5. The performance of the proposed approach is demonstrated in Section 6 on two nonlinear dynamical problems: a nonlinear vibration absorber and a nonlinear Jeffcott rotor. Finally, conclusions are drawn in the last section.

Equilibrium path

Harmonic Balance Method

A forced nonlinear dynamical system with n degrees of freedom (DOFs) governed by the following set of equations of motion is considered

r(x, ω, t) = M ẍ(t) + C ẋ(t) + Kx(t) + f nl (x, ẋ) -p(ω, t) = 0 (1)
with x(t) the displacement vector containing the n DOFs, M, C, K the generalized mass, damping and stiffness matrices, f nl the vector of nonlinear forces, p(ω, t) the periodic excitation forces and ω the excitation frequency. Harmonic balance method is used for its efficiency compared to time integration to predict periodic solutions. Assuming that the harmonic response has the same period τ = 2π/ω as the excitation, the displacements, nonlinear forces and excitation forces can be represented as Fourier series up to order H

x(t) = X 0 + H k=0 X k c cos(kωt) + X k s sin(kωt) f nl (x, ẋ) = F 0 + H k=0 F k c cos(kωt) + F k s sin(kωt) (2)
p(ω, t) = P 0 + H k=0 P k c cos(kωt) + P k s sin(kωt)
where

X = X 0 T , X 1 c T , X 1 s T , . . . , X H c T , X H s T T F nl = F 0 T , F 1 c T , F 1 s T , . . . , F H c T , F H s T T (3)
P = P 0 T , P 1 c T , P 1 s T , . . . , P H c T , P H s T T
are vectors of Fourier coefficients of size L = n × (2H + 1) . Equation ( 2) can be rewritten in the compact form

x(t) = (T(ωt) ⊗ I n )X f nl (x, ẋ) = (T(ωt) ⊗ I n )F nl (4) p(ω, t) = (T(ωt) ⊗ I n )P
where the linear operator T stands for the inverse Fourier transform and is made of trigonometric functions

T(ωt) = [1 cos(ωt) sin(ωt) . . . cos(Hωt) sin(Hωt)] ( 5
)
I n is the n × n identity matrix and ⊗ is the Kronecker tensor product. Similarly, introducing the derivative operator

∇ = diag(0, ∇ 1 , ..., ∇ j , ..., ∇ H ) with ∇ j = j 0 1 -1 0 (6)
the velocity and acceleration vectors can be written as follows

ẋ(t) = ω(T(ωt) ⊗ I n )(∇ ⊗ I n )X = ω[(T(ωt)∇) ⊗ I n ]X (7) ẍ(t) = ω 2 (T(ωt) ⊗ I n )(∇ 2 ⊗ I n )X = ω 2 [(T(ωt)∇ 2 ) ⊗ I n ]X
A Galerkin procedure is then applied to the equations of motion, which consists in projecting Eq.( 1) on the trigonometric basis T(ωt) with respect to the scalar product

f, g = 2 τ τ 0 f.g dt (8)
By introducing Eqs.( 4) and ( 7) into the resulting equations and using the orthogonality property T(ωt), T(ωt) = I 2H+1 , the nonlinear system of n equations (1) in the time domain is transformed into a nonlinear algebraic system of dimension L in the frequency domain

R(X, ω) = Z(ω)X + F nl (X) -P = 0 (9)
where

Z(ω) = ω 2 ∇ 2 ⊗ M + ω∇ ⊗ C + I 2H+1 ⊗ K = diag(K, Z 1 , ..Z j , ..Z H ) (10)
Z j = K -j 2 ω 2 M ωC -ωC K -j 2 ω 2 M (11)
An incremental-iterative Newton-Raphson procedure is used to obtain solutions of Eq.( 9), with corrections δX given by

R k X δX = -R k (12)
X k+1 = X k + δX (13)
where R X stands for the derivative of R with respect to X and the superscript k indicates an evaluation with updated variables at iteration k, i.e.

R k = R(X k ) R k X = ∂R ∂X X=X k (14)

Continuation technique and path-following

Nonlinear systems often have several possible responses for a given excitation frequency ω. The pseudo-arc length continuation method [START_REF] Crisfield | A fast incremental/iterative solution procedure that handles snap-through[END_REF] combined with the above-mentioned algorithm permits following the solution branch beyond limit points in order to obtain both stable and unstable solutions of the response curve. It is based on tangent prediction and orthogonal corrections [START_REF] Seydel | Practical Bifurcation and Stability Analysis[END_REF]. For this purpose, not only X but also ω is considered as an unknown. Once a starting point (X 0 , ω 0 ) on the solution curve has been obtained by solving Eq.( 12), a tangent vector t = (∆X, ∆ω) is obtained from

[R X R ω ] ∆X ∆ω = 0 (15)
where R ω = ∂R/∂ω. A normalization condition t = ∆X T ∆X + ∆ω 2 = 1 is added to make the tangent vector unitary and the prediction step gives a solution at distance ∆s along direction t

X 1 ω 1 = X 0 ω 0 + ∆s t (16)
Then, this prediction is corrected iteratively in the direction orthogonal to t until the convergence criterion R k ≤ ε is satisfied, with ε a user-defined accuracy. The corrections δX k+1 and δω k+1 are given by the system of L + 1 equations

R k X R k ω ∆X T ∆ω δX k+1 δω k+1 = -R k 0 ( 17
)
and the solution at iteration k + 1 is obtained with

X k+1 = X k + δX k+1 ω k+1 = ω k + δω k+1 (18)
Finally, an adaptive step-length algorithm [START_REF] Seydel | Practical Bifurcation and Stability Analysis[END_REF] based on the number of iterations in the corrector step is used to ensure robust continuation.

Computation of nonlinear terms and derivatives

Nonlinear forces are usually much easier to evaluate in the time domain than in the frequency domain. Therefore, the Fourier coefficients of nonlinear efforts F nl and their derivative ∂F nl ∂X involved in the Newton-Raphson iterations are obtained by Alternating Frequency-Time (AFT) method [START_REF] Cameron | An alternating frequency/time domain method for calculating the steady-state response of nonlinear dynamic systems[END_REF]. The AFT algorithm uses discrete Fourier transforms (DFT) to compute the nonlinear forces in the time domain and then switches back to the frequency domain

X DFT -1 -----→ x(t), ẋ(t) ---→ f nl (x, ẋ) DFT ---→ F nl (X) (19)
Defining N uniformly spaced time samples t i = i ∆t, i = 1 . . . N, such that ∆t = τ/N with N ≥ 2H + 1 and the corresponding vectors of N time samples of the displacements and nonlinear forces

x = [x(t 1 ), . . . , x(t N )] T fnl = [ f nl (t 1 ), . . . , f nl (t N )] T (20)
and using Eq.( 4) yields the following inverse-DFT relations between the vectors of time samples and the Fourier coefficients [START_REF] Narayanan | A frequency domain based numerical-analytical method for non-linear dynamical systems[END_REF] x = [T(ωt 1 ) ⊗ I n , . . . , T(ωt

N ) ⊗ I n ] T X = (Γ ⊗ I n )X x = [T(ωt 1 ) ⊗ I n , . . . , T(ωt N ) ⊗ I n ] T (∇ ⊗ I n )X = ω(Γ ⊗ I n )(∇ ⊗ I n )X = ω[(Γ∇) ⊗ I n ]X (21)
Once fnl has been computed in the time domain using x and x, the following DFT relation is used to go back to the frequency domain

F nl = (Γ ⊗ I n ) -1 fnl = (Γ -1 ⊗ I n ) fnl (22)
Remarkably, matrices Γ and Γ -1 do not depend on ω and thus can be conveniently written as follows after introducing

θ i = ωt i = 2πi/N Γ =             1 cos θ 1 sin θ 1 . . . cos Hθ 1 sin Hθ 1 . . . . . . . . . . . . . . . 1 cos θ N sin θ N . . . cos Hθ N sin Hθ N             (23)
Γ -1 = 1 N                             1 . . . 1 2 cos θ 1 . . . 2 cos θ N 2 sin θ 1 . . . 2 sin θ N . . . . . . 2 cos Hθ 1 . . . 2 cos Hθ N 2 sin Hθ 1 . . . 2 sin Hθ N                             (24)
Using the chain rule, the derivative of F nl can then be obtained

∂F nl ∂X = ∂F nl ∂ fnl ∂ fnl ∂ x ∂ x ∂X + ∂F nl ∂ fnl ∂ fnl ∂ x ∂ x ∂X = (Γ -1 ⊗ I n ) ∂ fnl ∂ x (Γ ⊗ I n ) + (Γ -1 ⊗ I n ) ∂ fnl ∂ x ω[(Γ∇) ⊗ I n ] (25)
with the nN × nN block-diagonal matrices of time samples of f nl derivatives

∂ fnl ∂ x = diagblk ∂ f nl ∂x t=t 1 , . . . , ∂ f nl ∂x t=t N ∂ fnl ∂ x = diagblk ∂ f nl ∂ ẋ t=t 1 , . . . , ∂ f nl ∂ ẋ t=t N (26)
In the case of systems with distinct states, the differentiation of the nonlinear force F nl requires a specific treatment as detailed in [START_REF] Krack | A high-order harmonic balance method for systems with distinct states[END_REF].

The Jacobian R X is eventually given by

R X = Z(ω) + ∂F nl ∂X (27)
The expression of R ω is immediate from Eqs.( 9) and ( 10)

R ω = Z ω X (28)
where

Z ω = ∂Z ∂ω = 2ω∇ 2 ⊗ M + ∇ ⊗ C (29)

Stability analysis and bifurcation points

The continuation technique permits to find several solutions on the response curve for one given excitation. However, some of them are not stable and can not be experimentally reproduced. In addition, only periodic solutions are found with the HBM described in the previous section. As a result, stability analysis and detection of bifurcations indicating changes of dynamical regime are essential tools in numerical simulations.

Stability of solutions can be examined by means of Hill's method which is variant of Floquet theory for frequencydomain methods such as HBM [START_REF] Von Groll | The harmonic balance method with arc-length continuation in rotor/stator contact problems[END_REF][START_REF] Peletan | A comparison of stability computational methods for periodic solution of nonlinear problems with application to rotordynamics[END_REF]. Hill's method consists in perturbing the equations of motion (1) around a periodic solution and applying the same Galerkin procedure as in Section 2.1 for HBM. The following quadratic eigenvalue problem is obtained

(R X + Λ∆ 1 + Λ 2 ∆ 2 )φ = 0 (30)
where φ = φ 1 + iφ 2 are complex eigenvectors, R X is the Jacobian defined by Eq.( 27) and

∆ 1 = 2ω∇ ⊗ M + I 2H+1 ⊗ C = diag C, C 2ωM -2ωM C , . . . , C 2HωM -2HωM C (31)
∆ 2 = I 2H+1 ⊗ M (32)
It should be noted that 2L = 2n(2H + 1) complex eigenvalues Λ are obtained from Eq.( 30), of which only the 2n ones with the smallest imaginary part are expected to have a physical meaning and correspond to the Floquet exponents [START_REF] Moore | Floquet theory as a computational tool[END_REF] denoted λ i , i = 1 . . . 2n in the following. By definition, a singular point appears when a Floquet exponent λ is null (see Fig. 1), i.e. when the determinant ϕ = det(R X ) vanishes. The indicator ϕ is thus monitored during the continuation of the response curve and a change of sign of ϕ between two consecutive points indicates the presence of a limit or branch point. The type of solution can be classified as follows [START_REF] Seydel | Practical Bifurcation and Stability Analysis[END_REF][START_REF] Govaerts | Numerical Methods for Bifurcations of Dynamical Equilibria[END_REF] where φ is the eigenvector associated with the null eigenvalue Λ = 0 of R X . A Neimark-Sacker (NS) or secondary Hopf bifurcation indicates a transition from periodic to quasi-periodic regime (see Fig. 2) for solutions of the timedomain equation [START_REF] Nayfeh | Nonlinear oscillations[END_REF]. Its analogue in the framework of the frequency-domain HBM is a Hopf bifurcation which indicates a transition from a fixed point (solution with constant Fourier coefficients) to a limit cycle (solution with periodic Fourier coefficients). Considering this equivalence, both denominations will be used indifferently in the following. A NS bifurcation occurs when a complex conjugate pair of Floquet exponents crosses the imaginary axis of the complex plane with λ = ±iκ, κ π/τ (see Fig. 1). The presence of a NS bifurcation along the equilibrium curve can thus be detected by monitoring the evolution of the Floquet exponents. Indicators which vanish at a NS bifurcation are available [START_REF] Seydel | Practical Bifurcation and Stability Analysis[END_REF][START_REF] Govaerts | Numerical Methods for Bifurcations of Dynamical Equilibria[END_REF][START_REF] Kuznetsov | Elements of applied bifurcation theory[END_REF]

Regular point if ϕ = det(R X ) 0 and R T ω φ 0 Limit point (LP) / Fold bifurcation if ϕ = det(R X ) = 0 and R T ω φ 0 Branch point (BP) if ϕ = det(R X ) = 0 and R T ω φ = 0 (33)
ϕ NS 1 = 1≤ j<i≤n λ i + λ j ϕ NS 2 = det(2 A I 2n ) (34)
where A is the diagonal matrix of the 2n Floquet exponents A = diag(λ 1 , . . . , λ 2n ) and stands for the bialternate product. These two indicators are equivalent because matrix 2A I 2n has eigenvalues λ i + λ j . However, one has to be careful when using them since they vanish not only for NS bifurcations but also for a pair of opposite real Floquet exponents λ = ±κ indicating a neutral saddle point.

Localization of bifurcation points

As stated in the introduction, two approaches exist for the precise computation of bifurcation points. The first one relies on so-called standard extended systems and consists in introducing one or more additional equations characterizing the bifurcation [START_REF] Seydel | Practical Bifurcation and Stability Analysis[END_REF][START_REF] Moore | The calculation of turning points of nonlinear equations[END_REF]. The second approach relies on bordering techniques and minimally extended systems in which only one scalar function is added [START_REF] Govaerts | Numerical Methods for Bifurcations of Dynamical Equilibria[END_REF][START_REF] Kuznetsov | Elements of applied bifurcation theory[END_REF][START_REF] Detroux | The harmonic balance method for bifurcation analysis of large-scale nonlinear mechanical systems[END_REF]. The approach based on standard extended systems is used in the present paper. The main motivation for this choice is the ease of implementation since the additional equations involved in standard extended systems have the same form and thus involve the same derivatives whatever the type of bifurcation, as shown in the next sections. The counterpart is an increase in the size of the problem, as compared to minimally extended systems, but the extra computational cost can be reduced with the use of a block elimination algorithm [START_REF] Keller | Numerical solution of bifurcation and nonlinear eigenvalue problems[END_REF].

The following procedure is used to locate bifurcation points. When a change of sign of indicators ϕ or ϕ NS is detected between two consecutive solution points, the point with the smallest absolute value of the indicator is used as a starting point (X 0 , ω 0 ) for the fully extended system.

Localization of limit points (LP) / Fold bifurcations

To locate limit points precisely, it is computationally more efficient to use the equation R X y = 0, obtained by substituting Λ = 0 into [START_REF] Petrov | Analysis of bifurcations in multiharmonic analysis of nonlinear forced vibrations of gas-turbine engine structures with friction and gaps[END_REF], rather than the test function ϕ = det(R X ) = 0 for the additional equation of the extended system. The resulting 2L + 1 extended system for localizing LPs takes the form

G(Y) =           R(X, ω) R X φ φ T φ -1           = 0 2L+1 (35)
with Y = (X, φ, ω). The first equation defines the equilibrium path, the second one characterizes the limits points and the last one normalizes the eigenvector φ. This extended system was first introduced by Moore and Spence [START_REF] Moore | The calculation of turning points of nonlinear equations[END_REF]. System [START_REF] Griewank | The calculation of hopf points by a direct method[END_REF] can be solved by a Newton-Raphson method, with corrections at iteration k given by

G Y (Y k ) δY = -G(Y k ) ( 36
)
Y k+1 = Y k + δY (37)
with δY = (δX, δφ, δω). Eq.( 36) can be written explicitly as

          R k X 0 L×L R k ω (R X φ) k X R k X (R X φ) k ω 0 T L 2φ k T 0                     δX δφ δω           = -           R k R k X φ k φ k T φ k -1           (38)
In practice, the eigenvector corresponding to the smallest eigenvalue of R X (X 0 , ω 0 ) can be used as a starting value φ 0 for φ. The 2L + 1 extended system (38) is regular at limit points and can be inverted. If L is large, a block elimination algorithm [START_REF] Keller | Numerical solution of bifurcation and nonlinear eigenvalue problems[END_REF] can be used in order to deal with sub-systems of size L only, involving R X , and thus reduce the computational work. In this case, one has to deal with R X becoming singular when approaching the limit point [START_REF] Cardona | Evaluation of simple bifurcation points and post-critical path in large finite rotation problems[END_REF][START_REF] Felippa | Traversing critical points by penalty springs[END_REF].

Localization of branch points (BP)

Based on the classification [START_REF] Battini | Improved minimal augmentation procedure for the direct computation of critical points[END_REF], the following additional equation in introduced in order to discriminate branch points from limit points R T ω φ = 0 (39)

Contrary to limit points, system ( 38) is singular at branch points. In order to avoid the singularity when approaching branch points, a penalty term γe j is introduced where e j is a unit vector with j-th component equal to 1 and γ is an auxiliary variable [START_REF] Felippa | Traversing critical points by penalty springs[END_REF]. In practice, the initial value initial value for γ is assumed to be 0 and the eigenvector corresponding to the smallest eigenvalue of R X (X 0 , ω 0 ) can be used as a starting value φ 0 for φ. Consequently, the 2L + 2 extended system to be solved during the Newton-Raphson procedure is given by

               R k X 0 L×L R ω e j (R X φ) k X R k X (R X φ) k ω 0 L 0 T L 2φ k T 0 0 0 T L R T ω 0 0                              δX δφ δω δγ               = -                 R k + γ k e j R k X φ k φ k T φ k -1 R k ω T φ k                 (40)
Again, a block elimination algorithm can be used in order to solve sub-systems of size L involving R X only, with a particular attention on R X becoming progressively ill-conditioned [START_REF] Petrov | Analysis of bifurcations in multiharmonic analysis of nonlinear forced vibrations of gas-turbine engine structures with friction and gaps[END_REF].

Localization of Neimark-Sacker (NS) bifurcations

The original extended system presented here takes inspiration from the approach of Griewank and Reddien [START_REF] Griewank | The calculation of hopf points by a direct method[END_REF] which has been adapted to work within the framework of the HBM.

By substituting Λ = ±iκ and φ = φ 1 ±iφ 2 into (30) and identifying real and imaginary parts, two equations characterizing NS bifurcations are obtained. Along with equilibrium and normalization equations, they provide the 3L + 2 extended system for locating NS bifurcations

K(X, φ 1 , φ 2 , κ, ω) =                    R(X, ω) R X φ 1 -κ∆ 1 φ 2 -κ 2 ∆ 2 φ 1 R X φ 2 + κ∆ 1 φ 1 -κ 2 ∆ 2 φ 2 q T φ 1 φ T 1 φ 1 -1                    = 0 3L+2 (41)
where q a constant vector with non zero projection on span{φ 1 , φ 2 }. It is worth noting that other normalizations are possible, which can affect the algorithm behaviour. For instance, replacing q T φ 1 = 0 and φ T 1 φ 1 = 1 by q T φ 1 = 0 and q T φ 2 = 1 in (41) leads to other solutions. In this case, for κ = 0, one has φ 1 = 0 and the augmented system is reduced to [START_REF] Griewank | The calculation of hopf points by a direct method[END_REF] and therefore limit points are obtained [START_REF] Griewank | The calculation of hopf points by a direct method[END_REF] [START_REF] Roose | A direct method for the computation of hopf bifurcation points[END_REF]. Obviously, these normalizations should be avoided if only NS bifurcations are tracked, but they can be of interest for switching from a branch of LPs to a branch of NS bifurcations.

Once a starting point (X 0 , ω 0 ) on the equilibrium path in the neighbourhood of the NS bifurcation is available, the eigenvector corresponding to the complex eigenvalue Λ = iκ of R X (X 0 , ω 0 ) can be used as a starting value for φ = φ 1 +iφ 2 . This initial solution is then corrected by means of Newton-Raphson iterations resulting from the linearization of Eq.( 41)

                    R k X 0 L×L 0 L×L 0 L R k ω (R X φ 1 ) k X R k X -(κ k ) 2 ∆ 2 -κ k ∆ k 1 -∆ k 1 φ k 2 -2κ k ∆ 2 φ k 1 (R X φ 1 ) k ω -κ k (∆ 1 φ 2 ) k ω (R X φ 2 ) k X κ k ∆ k 1 R k X -(κ k ) 2 ∆ 2 +∆ k 1 φ k 1 -2κ k ∆ 2 φ k 2 (R X φ 2 ) k ω + κ k (∆ 1 φ 1 ) k ω 0 T L q T 0 T L 0 0 0 T L 2φ k 1 T 0 T L 0 0                                        δX δφ 1 δφ 2 δκ δω                    = -                     R k R k X φ k 1 -κ k ∆ k 1 φ k 2 -(κ k ) 2 ∆ k 2 φ k 1 R k X φ k 2 + κ k ∆ k 1 φ k 1 -(κ k ) 2 ∆ k 2 φ k 2 q T φ k 1 φ k 1 T φ k 1 -1                     (42)
Again, a block elimination algorithm can be used in order to solve sub-systems of size L only. This time, no particular procedure is necessary for R X since it is non-singular at NS bifurcations.

Computation of the additional derivatives involved the extended systems

Expressions for the derivatives R X and R ω are given in Eqs. ( 27) and [START_REF] Cardona | Evaluation of simple bifurcation points and post-critical path in large finite rotation problems[END_REF]. Since the vector of Fourier coefficients φ does not depend on ω, the derivative (R X φ) ω is given by

(R X φ) ω = (R X ) ω φ = Z ω φ (43)
with Z ω calculated with Eq. [START_REF] Lopez | Post-critical analysis of structures with a nonlinear pre-buckling state in the presence of imperfections[END_REF]. Similarly, the derivatives involving ∆ 1 are obtained from Eq.(31) as follows

(∆ 1 φ) ω = (∆ 1 ) ω φ = 2(∇ ⊗ M)φ (44)
The derivative (R X φ) X can be computed by means of the AFT method much in the same way as for

∂F nl ∂X in Section 2.3 (R X φ) X = (Γ -1 ⊗ I n ) ∂r x ϕ ∂ x (Γ ⊗ I n ) ( 45
)
with the nN × nN block-diagonal matrix of time samples

∂r x ϕ ∂ x = diagblk ∂r x ϕ ∂x t=t 1 , . . . , ∂r x ϕ ∂x t=t N ( 46
)
where r x is the derivative of the equilibrium residual (1) in the time domain and the eigenmode ϕ in the time domain can be obtained from φ = [ϕ(t 1 ), . . . , ϕ(t

N )] T = (Γ ⊗ I n )φ (47)
(r x ϕ) x can be obtained analytically in very simple cases only. For complicated systems or systems with many DOFs, it can be approximated by finite differences [START_REF] Wriggers | A general procedure for the direct computation of turning and bifurcation points[END_REF] (r x ϕ) x = ∂r x ϕ ∂x

1 ε x [r x (x + ε x ϕ) -r x (x)] (48)
Similarly, (R X φ) X can also be computed directly in the frequency domain as

(R X φ) X 1 ε X [R X (X + ε X φ) -R X (X)] (49)
In Eqs. ( 48) and ( 49), the proper scaling of the perturbations ε x ϕ and ε X φ with respect to x and X is essential to ensure good accuracy and robustness of the numerical differentiation. This is accomplished by means of the following formulas in which a suitable value for η has to be chosen. Using η = 10 -6 has been found to provide good results for a large range of applications.

ε x = η x ϕ + η ε X = η X φ + η (50)

Parametric continuation of bifurcation points

The following computational procedure is used to follow a branch of bifurcation points :

1. For a fixed value α 0 of a system parameter, nonlinearity or excitation level for instance, follow the equilibrium path (see Section 2.2) until a bifurcation point (X 0 , ω 0 , α 0 ) is encountered and precisely detected (Section 4), which serves as a starting point.

2. Free α and follow the path of bifurcations by a continuation method.

The last step is carried out by adding a pseudo-arc length constraint equation to one of the extended systems (38)(40) [START_REF] Hussein | On efficient methods for detecting hopf bifurcation with applications to power system instability prediction[END_REF] depending on the type on bifurcation to be followed. As in Section 2.2, a tangent vector t = (∆X, ∆ω, ∆α) is first computed, then the next bifurcation point if obtained after Newton-Raphson corrections orthogonal to t.

In the case of LP tracking, these corrections are solutions of the 2L + 2 extended system

               R k X 0 L×L R k ω R k α (R X φ) k X R k X (R X φ) k ω (R X φ) k α 0 T L 2φ k T 0 0 ∆X T 0 T L ∆ω ∆α                              δX δφ δω δα               = -                R k R k X φ k φ k T φ k -1 0                (51)
and the two additional derivatives R α and (R X φ) α are given by

R α = ∂Z ∂α X + ∂F nl ∂α - ∂P ∂α (R X φ) α = ∂Z ∂α + ∂ 2 F nl ∂X∂α φ (52)
Analytical expressions for these two derivatives depend on the system parameter α under consideration and vary with each case. In order to compute them conveniently for any system parameter, finite differences can be used with

ε α = η(|α| + η) and η = 10 -6 R α (X, ω, α) 1 ε α [R(X, ω, α + ε α ) -R(X, ω, α)] ( 53
) (R X φ) α (X, ω, α) 1 ε α [R X (X, ω, α + ε α ) -R X (X, ω, α)]φ (54)
Similarly, in the case of BP tracking, Newton-Raphson corrections are solutions of the 2L + 3 extended system

                     R k X 0 L×L R ω e j R k α (R X φ) k X R k X (R X φ) k ω 0 L (R X φ) k α 0 T L 2φ k T 0 0 0 0 T L R T ω 0 0 0 ∆X T 0 T L ∆ω 0 ∆α                                         δX δφ δω δγ δα                    = -                      R k + γ k e j R k X φ k φ k T φ k -1 R k ω T φ k 0                      (55)
Finally, in the case of NS tracking, Newton-Raphson corrections are solutions of the 3L + 3 extended system

                          R k X 0 L×L 0 L×L 0 L R k ω R k α (R X φ 1 ) k X R k X -(κ k ) 2 ∆ 2 -κ k ∆ k 1 -∆ k 1 φ k 2 -2κ k ∆ 2 φ k 1 (R X φ 1 ) k ω -κ k (∆ 1 φ 2 ) k ω (R X φ 1 ) k α -κ k (∆ 1 ) k α φ k 2 -(κ k ) 2 (∆ 2 ) k α φ k 1 (R X φ 2 ) k X κ k ∆ k 1 R k X -(κ k ) 2 ∆ 2 +∆ k 1 φ k 1 -2κ k ∆ 2 φ k 2 (R X φ 2 ) k ω + κ k (∆ 1 φ 1 ) k ω (R X φ 2 ) k α + κ k (∆ 1 ) k α φ k 1 -(κ k ) 2 (∆ 2 ) k α φ k 2 0 T L q T 0 T L 0 0 0 0 T L 2φ k 1 T 0 T L 0 0 0 ∆X T 0 T L 0 T L ∆ω 0 ∆α                           ×                          δX δφ 1 δφ 2 δκ δω δα                          = -                           R k R k X φ k 1 -κ k ∆ k 1 φ k 2 -(κ k ) 2 ∆ k 2 φ k 1 R k X φ k 2 + κ k ∆ k 1 φ k 1 -(κ k ) 2 ∆ k 2 φ k 2 q T φ k 1 φ k 1 T φ k 1 -1 0                           (56)
where the derivatives (∆ 1 ) α and (∆ 2 ) α are not null only if the system parameter α is involved in matrix M or C.

Applications

The proposed bifurcation tracking technique is applied to a nonlinear energy sink (NES) and a nonlinear Jeffcott rotor in order to demonstrate its capabilities and performance. These applications also aim at showing that bifurcation tracking can be used as an efficient tool for the design and performance tuning of nonlinear dynamical systems.

Nonlinear energy sink (NES)

The first application concerns a linear oscillator with an attached nonlinear energy sink (NES) under harmonic external forcing. Unlike common linear vibration absorbers, NES are efficient in a wide range of frequencies and over a range of amplitudes of the external forcing [START_REF] Vakakis | Nonlinear targeted energy transfer in mechanical and structural systems[END_REF]. This is due mainly to a qualitatively different dynamical behaviour in close vicinity of the main resonance, where the NES can exhibit quasi-periodic rather than periodic response [START_REF] Gendelman | Quasiperiodic energy pumping in coupled oscillators under periodic forcing[END_REF][START_REF] Gendelman | Quasi-periodic response regimes of linear oscillator coupled to nonlinear energy sink under periodic forcing[END_REF]. Indeed, the energy of vibration is transferred to the NES and damped out in quasi-periodic regime, leading to energy level reduction and attenuation of vibrations of the primary oscillator.

In the case of strongly nonlinear vibration absorber, it was shown in [START_REF] Starosvetsky | Attractors of harmonically forced linear oscillator with attached nonlinear energy sink. ii: Optimization of a nonlinear vibration absorber[END_REF] that, depending on system parameters, weak or strong quasi-periodic regimes occur, both of them providing efficient vibration absorption. Thus, in order to perform a careful tuning, it is essential to determine the range of parameters in which quasi-periodic beating responses are possible. Since a weak quasi-periodic response relates to a periodic response loosing its stability via NS bifurcation, the boundaries for the quasi-periodic response can be obtained by the parametric continuation of NS bifurcations introduced in this paper.

The system under consideration consists of a linear oscillator and a small essentially nonlinear attachment, as described in [START_REF] Starosvetsky | Attractors of harmonically forced linear oscillator with attached nonlinear energy sink. ii: Optimization of a nonlinear vibration absorber[END_REF]. The linear oscillator is the primary structure with mass m 1 and the nonlinear attachment with small mass m2 is the NES. They are coupled through a pure cubic elastic force and a viscous damper. For simplicity, the stiffness of the linear oscillator is chosen such that its resonance frequency is equal to 1. The system is governed by the following equations ẍ1 + λ( ẋ1 -ẋ2 )

+ x 1 + k nl (x 1 -x 2 ) 3 = Acos(ωt) ẍ2 + λ( ẋ2 -ẋ1 ) + k nl (x 2 -x 1 ) 3 = 0 ( 57
)
with x 1 and x 2 the displacements of the linear oscillator and the absorber (NES), respectively, = m 2 /m 1 << 1 the mass ratio, λ the damping coefficient, k nl the nonlinear stiffness coefficient and A the amplitude of the external force. In the following, = 0.1 and λ = 0.4 are fixed while k nl and A are varied for tuning purpose. First, A = 0.3 is kept fixed and the dynamics of the coupled nonlinear system is investigated for several values of k nl . The forced response is computed by HBM and continuation technique as described in Section 2. Convergence of the results is obtained with H = 3 harmonics and verified by comparison with time integration. Stability analysis and detection of bifurcations are performed according to Sections 3 and 4. The amplitude x 1 of the linear oscillator is plotted in Figure 3, where solid lines stand for stable periodic solutions, dashed lines for unstable ones and triangle markers for NS bifurcations. For k nl =1, a periodic response is obtained for all frequencies. For greater values of k nl , the NS bifurcations clearly define frequency zones in which the quasi-periodic beating response is the only stable response. The system response in these quasi-periodic zones can be easily computed by time integration. Since its amplitude is modulated, the total average system energy

E tot = ẋ2 1 2 + ẋ2 2 2 + x 2 1 2 + k nl (x 1 -x 2 ) 4 4 t ( 58
)
is used instead, with the average performed over a period of time 50 times longer than the modulation period. The total system energy E tot vs. excitation frequency ω is plotted in Fig. 4. In this figure, the unstable periodic responses have been replaced by stable quasi-periodic ones in the zones of quasi-periodic beating. It can be observed that the level of total system energy is reduced when k nl is increased, therefore demonstrating the absorber efficiency.

Obtaining the boundaries for the quasi-periodic beating response by this way is a tedious task, since the computation of frequency response and the detection of bifurcations have to be repeated for each value of k nl in the range of interest. An efficient alternative consists in using the direct bifurcation tracking presented in Section 5. To do so, a response curve is first computed for an arbitrary fixed value of k nl , e.g. k nl =3, until the first NS bifurcation is detected. Then, k nl is considered as a new unknown and this NS bifurcation is used as a starting point for the direct continuation of the curve of NS bifurcations. Doing so, the boundaries for the quasi-periodic beating response are obtained with only one computation. This curve is represented as a 3D-plot in Fig. 5. The frequency responses plotted in Fig. 5 for some values of k nl were not used for the computations. They are presented here to make the interpretation of the 3D-plot easier and to show that the direct NS tracking works as expected and is accurate.

Projections of the NS tracking curve on the E totk nl and k nl -ω planes are plotted in Fig. 6. They provide useful information for designing the NES. The areas inside and outside the NS tracking curve correspond to quasi-periodic beating responses and to periodic responses without beating, respectively. It can be seen from Fig. 6(a) that the total system energy for beating responses decreases monotonically with k nl increasing. Moreover, for a given value of k nl , the system energy level is comprised between the two values given by the NS tracking curve. From Fig. 6(b), the range of frequencies in which a beating response occurs can be directly read for each value of k nl . It can also be observed that in zone (A) for k nl < 1.92 there is no possibility of a quasi-periodic beating response. In zone (B) for 1.92 < k nl < 7.78, there is a wide range of frequencies centered on the resonance. This operating zone is optimal and the highest value of k nl in this zone (k nl 7.78) appears to be the optimal one. In zone (C) for 7.78 < k nl < 11.28, the quasi-periodic area is divided into two smaller areas separated by a periodic area. The left-hand side quasi-periodic area is really narrow and the right-hand side one is not centered on the resonance. As a consequence, this is not an optimal operating zone. In zone (D) for k nl > 11.28, the quasi-periodic area is not centered on the resonance. In addition, the frequency responses for these values of k nl are very complex and there may be other stable solutions than the expected quasi-periodic response. Thus, this is not a safe operating zone. The amplitude A of the external force is now used as the varying parameter while other system parameters =0.1, λ=0.4 and k nl =5 are kept fixed. The total system energy vs. excitation frequency ω is plotted in Fig. 7 for A=0.1, 0.2, 0.3, 0.4. As expected, it can be observed that, away from the resonance where the response is periodic, the level of total system energy increases with A. However, in the vicinity of the resonance where quasi-periodic beating response occurs, the average energy level is almost the same whatever the value of A, therefore demonstrating the constant absorber efficiency over a range of amplitudes of the external forcing. In order to determine the boundaries of the quasi-periodic beating area, NS tracking is performed, starting from the first detected NS bifurcation on the frequency curve for A=0.3. The resulting tracking curve in presented in the 3D-plot of Fig. 8 and its projections on the E tot -A and A -ω planes are plotted in Fig. 9. Again, four operating zones can be identified.

In zone (A) for A < 0.186 there is no possibility of a quasi-periodic beating response. Zone (B), for 0.186 < A < 0.374, defines an optimal operating zone since it contains a wide range of frequencies centred on the resonance. In zone (C) for 0.374 < A < 0.45, the main quasi-periodic area is not centred any more on the resonance. As a consequence, this is not an optimal operating zone. Zone (D), for A > 0.45, is not a safe operating zone because the quasi-periodic area is not centred on the resonance and there may be other stable solutions than the expected quasi-periodic response.

From the computations presented in this example, it appears that the direct parametric continuation of bifurcation points exhibits a considerable efficiency and can be used as a fast tuning tool for the design of NES systems, since the optimal range of operating parameters can be determined with only one computation.

Nonlinear Jeffcott rotor

The second test case is the modified Jeffcott rotor presented in Fig. 10 which can interact at the disk location with a stator modelled as an additional stiffness [START_REF] Peletan | Quasi-periodic harmonic balance method for rubbing self-induced vibrations in rotor-stator dynamics[END_REF][START_REF] Jiang | Determination of the global responses characteristics of a piecewise smooth dynamical system with contact[END_REF]. The rotor is made of a weightless shaft carrying a disk with mass m and radius R disc at the middle of the span. The clearance between the rotor and the stator is denoted by h. The mass center of the rotor is located at a distance e from its geometrical center. The stator, which is rigidly fixed, has an elastic contact surface modelled as a symmetrical set of radial springs with isotropic stiffness k c . Let the investigated motion be located in the first natural frequency range. Therefore there is no gyroscopic effect to take into account in the equations of motion at the disk location described by

m ẍ + c ẋ + kx + k c (1 - h r )(x -µy sign(v rel )) = meω 2 cos ωt mÿ + cẏ + ky + k c (1 - h r )(µx sign(v rel ) + y) = meω 2 sin ωt (59)
with k the stiffness of the shaft, r = x 2 + y 2 the radial displacement, meω 2 the mass unbalance amplitude and v rel = ( x r ẏ -y r ẋ) + R disc ω the relative velocity between the rotor and the stator at the contact point. When r < 0, there is no rotor-stator contact and k c is set as nil.

All the calculations are carried out with the same set of parameters as in [START_REF] Peletan | Quasi-periodic harmonic balance method for rubbing self-induced vibrations in rotor-stator dynamics[END_REF] and [START_REF] Jiang | Determination of the global responses characteristics of a piecewise smooth dynamical system with contact[END_REF] Firstly, the calculation of the system response to the unbalance force is performed for a friction coefficient µ = 0.2 by HBM and continuation technique. Based on comparison with results obtained by time integration, N = 15 harmonics are used for HBM computations. This high number of harmonics is required for a good representation of the transition between the no-contact and the contact solution. Continuation, stability analysis and detection of bifurcations, as described in Sections 3 and 4, are performed at the same time. The maximum amplitude of the dimensionless radial displacement of the disk r/h with respect to the dimensionless excitation frequency ω = ω/ω 0 is plotted in Fig. 11, where solid lines stand for stable periodic solutions, dashed lines for unstable ones, dotted lines for quasi-periodic stable solutions, circle markers for LP and triangle markers for NS bifurcations. Periodic and quasiperiodic solutions are calculated by HBM and time integration respectively. There is no contact between the rotor and the stator until r/h = 1 ( ω 0.154). Then, the amplitude exceeds the initial clearance and the rotor rubs permanently against the stator in a so-called synchronous full annular rub motion. This periodic rub motion remains stable until ω 0.289 where a NS bifurcation is detected. For higher values of ω, the periodic rub motion is unstable and the only stable solution is the quasi-periodic partial rub motion characterized by rebounds and intermittent contact between the rotor and the stator. This quasi-periodic partial rub motion remains stable until ω 0.46. At this stage, there is a sudden jump in amplitude and a backward whirl motion is triggered, which may severely damage the rotor system in a very short time. The behaviours for µ = 0.05 and µ = 0.11 (Fig. 12) are qualitatively different. For µ = 0.05, the periodic synchronous full annular rub motion is stable until the end of the contact indicated by the limit point at ω 1 and is the only possible solution during the rotor-stator contact. For µ = 0.11, the response is also periodic after the contact initiation but it does not remain stable until the limit point is reached. Indeed, two NS bifurcations are found at ω = 0.847 and ω = 0.963. Between these two frequencies, the only stable motions are the quasi-periodic partial rub and the no-contact motions.

From this phenomenological analysis, it appears that the nature of the motion is highly dependent on the friction coefficient µ and on the excitation frequency ω. Thus, a parametric analysis defining the boundaries of the different dynamical regimes with respect to these two parameters is of great interest for designers and operators. Again, this parametric analysis can be performed by means of the direct bifurcation tracking presented in Section 5. To do so, a response curve is first computed for µ=0.2, until a NS bifurcation or a LP is detected. Then, µ is considered as a new unknown and this NS bifurcation or LP is used as a starting point and the bifurcation tracking is performed in both directions (increasing and decreasing values of µ). Doing so, the boundaries for the no-rub and the quasi-periodic motions are obtained with only one computation. These LP and NS tracking curves are represented as a 3D-plot in Fig. 13. A few response curves are also plotted to make the interpretation of the 3D-plot easier.

The projections of the tracking curves on plane ω -µ plane are plotted in Figs. 14 and15. For a given value of µ, these curves provide the frequencies corresponding to the occurrence of LP and NS bifurcations. The horizontal lines in Fig. 14 for µ = 0.05, µ = 0.11 and µ = 0.2 permit to compare the results of the bifurcation tracking with the response curves of Figs. 11 and12. Their intersections with the vertical lines and the two tracking curves precisely delimit the different dynamical regimes encountered along the response curves. Thus, the two LP and NS tracking curves divide the ω -µ plane into several regions corresponding to the different possible dynamical regimes of the rotor. These regions are identified in Fig. 15. The ends of the NS tracking curve give the two frequencies ω = 0.154 and ω = 0.854 corresponding to the initiation and end of the contact. The NS tracking curve also corresponds to the boundary of the quasi-periodic regime. The motion is periodic (synchronous full annular rub motion) below this curve and quasi-periodic (partial rub motion with rebounds) above it. Consequently, for µ < 0.108, the motion is purely periodic and the rotor undergoes only full annular rub during the contact with the stator. For 0.108 < µ < 0.36, the motion during contact can be periodic or quasi-periodic, depending on ω, and the range of periodic motion decreases with µ increasing. For µ > 0.36, the motion during contact is purely quasi-periodic. Finally, it can observed that the obtained NS tracking curve is in good agreement with the analytical results from Jiang [START_REF] Jiang | Determination of the global responses characteristics of a piecewise smooth dynamical system with contact[END_REF].

The computational cost for each numerical procedure is summarised in Table 1. The computation of a single response curve requires roughly 50-60 continuation steps. For µ=0.2, 56 adaptive continuation steps with an average number of 3 Newton-Raphson iterations per step are necessary to compute the whole response curve with an accuracy ε = 10 -6 and assess stability by Hill's method. The computational time for one step of response curve continuation is normalized to 1 and used as a reference time. It can be inferred from this table that one step of LP and NS tracking is respectively 1.6 and 2.3 times more costly than one step of response curve continuation. This increase is related to the larger size of the augmented systems to be solved. For this example, the LP and NS tracking require 54 and 79 continuation steps respectively. Therefore, the computational costs for the complete LP and NS trackings appear to be about 1.5 and 3 times higher than for a complete response curve. However, these costs do not include the computation of the starting bifurcation point, which requires the continuation of the response curve for µ=0.2 until the detection of the bifurcation points. For this application, 40 steps and 17 steps of continuation (with cost 1) are necessary to reach and detect the LP and NS bifurcation respectively. Thus, the actual costs of the complete LP and NS trackings are 87+40=127 and 182+17=199 respectively, i.e. they are equivalent to the computation of 2.3 and 3.5 complete response curves.

CONCLUSIONS

An efficient frequency-domain method for the fast parametric analysis of bifurcations of nonlinear dynamical systems has been presented. Instead of computing several response curves and detecting the bifurcation points, the proposed method consists in directly tracking the bifurcation points when a system parameter such as the excitation or nonlinearity level is varied.

To this end, a so-called extended system comprising the equation of motion and an additional equation characterizing the bifurcation of interest is solved by means of the Harmonic Balance Method coupled with an arc-length continuation technique. In particular, an original extended system for the detection and tracking of Neimark-Sacker (secondary Hopf) bifurcations has been proposed. Even if the additional equations characterizing the bifurcations increase the size of the problem, it is shown that the computational cost can be reduced with the use of a block elimination algorithm. Moreover, the additional equations have the same form and involve the same derivatives whatever the type of bifurcation, thus making the numerical implementation easier.

The examples considered in this paper have shown the interest of the proposed bifurcation tracking for the parametric analysis of a nonlinear vibration absorber (NES) and a nonlinear rotor system. Indeed, with only one computation a map showing the stability boundaries and changes of dynamical regime has been obtained. Moreover, the bifurcation tracking has been shown to be only 2 or 3 times more costly than a conventional response curve calculation. As a result, it can be used by engineers at the design stage to tune nonlinear systems efficiently.

The numerical developments have been fulfilled in both Matlab and Cast3M [START_REF]CAST3M, Finite Element software[END_REF] softwares, paving the way for applications of the method to the nonlinear dynamics of structures modelled with finite elements.

Figure 1 :

 1 Figure 1: Floquet exponents crossing the imaginary axis on the complex plane indicate a bifurcation

Figure 2 :

 2 Figure 2: Bifurcation points : (a) LP: Limit point (b) BP: Branch point (c) NS: Neimark-Sacker bifurcation ; Solid/dashed line: stable/unstable branch

Figure 3 :

 3 Figure 3: Vibration amplitudes of the main mass x 1 vs. excitation frequency for =0.1, λ=0.4, A=0.3 and various nonlinear stiffness coefficients k nl =1, 3, 5, 7, 9. Stable and unstable periodic responses.

Figure 5 :Figure 6 :

 56 Figure 5: Frequency responses and Neimark-Sacker bifurcation tracking for λ=0.4, A=0.3 and varying k nl .

Figure 7 :

 7 Figure 7: Average total system energy for =0.1, λ=0.4, k nl =5 and A=0.1, 0.2, 0.3, 0.4. Stable periodic and quasi-periodic solutions.

Figure 8 :

 8 Figure 8: Frequency responses and Neimark-Sacker bifurcation tracking for =0.1, λ=0.4, k nl =5 and varying A.

Figure 9 :

 9 Figure 9: Projection of the Neimark-Sacker bifurcations tracking for =0.1, λ=0.4, k nl =5 and varying A.

Figure 10 :

 10 Figure 10: Nonlinear Jeffcott rotor with stator contact

  : m = 1 kg, c = 5 N.s/m, k = 100 N/m, k c = 2500 N/m, h = 0.105 m, e = 0.1 m, R disc = 20h, ω 0 = √ k c /m = 50 rad/s, except that the friction coefficient µ is considered as the variable parameter.

Figure 11 : 11 Figure 12 :

 111112 Figure 11: Forced response curve of the modified Jeffcott rotor for a friction coefficient µ = 0.2.

  ti o n c o e ff ic ie n t µ

Figure 13 :Figure 14 :Figure 15 :

 131415 Figure 13: Limit points tracking and Neimark-Sacker bifurcations tracking of the Jeffcott rotor as a function of friction coefficient µ.

  4, A=0.3 and various nonlinear stiffness coefficients k nl =1, 3, 5, 7, 9. Stable and unstable periodic responses. Average total system energy for =0.1, λ=0.4, A=0.3 and k nl =1, 3, 5, 7, 9. Stable periodic and quasiperiodic solutions.

		0.14
		0.12		Periodic Quasi-periodic				k nl =1
	Total system energy	0.06 0.08 0.1 0.04							k nl =3 k nl =5
		0.02							k nl =9	k nl =7
		0
			0.8		0.85	0.9	0.95	1	1.05	1.1	1.15	1.2
									Frequency ω
	0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0 0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 Figure 4: 0 Total system energy	2 2 2	4 4 4 k n l	6 6 6	8 8 8	10 10 10	12 12 12	0.8	0.9	Fr eq ue nc y ω 1 NS tracking 1.1 Stable Unstable NS	1.2

Table 1 :

 1 Computational cost of the bifurcation tracking for the Jeffcott rotor

	Type of computation	Number of steps Relative average time per step Total time
	Response curve + stability (µ=0.2)	56	1	56
	LP tracking	54	1.6	87
	NS tracking	79	2.3	182

ACKNOWLEDGMENTS

The authors would like to acknowledge the financial support of CEA Saclay.