Baptiste Bergeot
email: baptiste.bergeot@insa-cvl.fr

Sergio Bellizzi
email: bellizzi@lma.cnrs-mrs.fr

Steady-state regimes prediction of a multi-degree-of-freedom unstable dynamical system coupled to a set of nonlinear energy sinks

Keywords: Multi-degree-of-freedom unstable system, set of nonlinear energy sinks, passive mitigation, relaxation oscillations, mitigation limit, asymptotic analysis

A general method to predict the steady-state regimes of a multi-degree-of-freedom unstable vibrating system (the primary system) coupled to several nonlinear energy sinks (NESs) is proposed. The method has three main steps. The first step consists in the diagonalization of the primary underline linear system using the so-called biorthogonal transformation. Within the assumption of a primary system with only one unstable mode the dynamics of the diagonalized system is reduced ignoring the stable modes and keeping only the unstable mode. The complexification method is applied in the second step with the aim of obtaining the slowflow of the reduced system. Then, the third step is an asymptotic analysis of the slow-flow based geometric singular perturbation theory. The analysis shows that the critical manifold of the system can be reduced to a one dimensional parametric curve evolving in a multidimensional space. The shape and the stability properties of the critical manifold and the stability properties of the fixed points of the slow-flow provide an analytical tool to predict the nature of the possible steady-state regimes of the system. Finally, two examples are considered to evaluate the effectiveness and advancement of the proposed method. The method is first applied to the prediction of the mitigation limit of a breaking system subject to friction-induced vibrations coupled to two NESs, and next an airfoil model undergoing an aeroelastic instability coupled to a NESs setup (from one to four) is discussed. Theoretical results are compared, for validation purposes, to direct numerical integration of the system. The comparisons show good agreement.

Introduction

In the framework of passive vibration control, it is now recognized that the nonlinear absorbers based on the concept of Targeted Energy Transfer (TET) also known as Nonliner Energy Sink (NES) are good candidate to consider especially when low frequency and high level are concerned. A basic NES generally consists of a light mass, an essentially nonlinear spring and a viscous linear damper. Due to the essentially nonlinear stiffness (with no linear part), a NES can tune to any frequency content displayed by the primary vibrating structure even if its mass is small. The TET concept results from nonlinear mode bifurcations and it can be described as a 1:1 resonance capture [START_REF] Gendelman | Energy Pumping in Nonlinear Mechanical Oscillators: Part I-Dynamics of the Underlying Hamiltonian Systems[END_REF][START_REF] Vakakis | Energy pumping in nonlinear mechanical oscillators: Part II -Resonance capture[END_REF]. A first review of the state-of-the-art can be found in [START_REF] Vakatis | Nonlinear Targeted Energy Transfer in Mechanical and Structural Systems[END_REF] and more recently in [START_REF] Lu | Nonlinear dissipative devices in structural vibration control: A review[END_REF] as a part of the state-of-the-art about technologies for nonlinear dissipative devices including for example (without being exhaustive) polynomial, non-polynomial, piecewise linear and non-smooth stiffness NES, and also vibro-impact NES, track NES, bistable NES and magnetic NES. The generation of new technologies is always in progress. Recent studies include for example a novel archetype of NES enhanced by an inerter for vibration reduction [START_REF] Zhang | Nonlinear energy sink with inerter[END_REF][START_REF] Zhang | An inertial nonlinear energy sink[END_REF], a acoustic NES [START_REF] Gourdon | Targeted energy transfer from one acoustical mode to an helmholtz resonator with nonlinear behavior[END_REF] and a virbo-acoustic NES with a controlled acoustic load [START_REF] Bryk | Experimental study of a hybrid electro-acoustic nonlinear membrane absorber[END_REF] for noise attenuation and a grounded NES for rotor system vibration suppression [START_REF] Yao | Using grounded nonlinear energy sinks to suppress lateral vibration in rotor systems[END_REF].

NESs are also used to control dynamic instabilities. The possible suppression of the limit cycles of a Van der Pol oscillator utilizing a NES was first demonstrated numerically in [START_REF] Lee | Suppression of limit cycle oscillations in the van der Pol oscillator by means of passive non-linear energy sinks[END_REF] and next considered theoretically in [START_REF] Gendelman | Bifurcations of self-excitation regimes in a Van der Pol oscillator with a nonlinear energy sink[END_REF]. A series of papers was also dedicated to aeroelastic mitigation [START_REF] Lee | Suppression aeroelastic instability using broadband passive targeted energy transfers, part 1: Theory[END_REF][START_REF] Vaurigaud | Passive control of aeroelastic instability in a long span bridge model prone to coupled flutter using targeted energy transfer[END_REF][START_REF] Gendelman | Asymptotic analysis of passive nonlinear suppression of aeroelastic instabilities of a rigid wing in subsonic flow[END_REF][START_REF] Luongo | Aeroelastic instability analysis of nes-controlled systems via a mixed multiple scale/harmonic balance method[END_REF], helicopter groud resonance instability mitigation [START_REF] Bergeot | Analysis of steady-state response regimes of a helicopter ground resonance model including a non-linear energy sink attachment[END_REF][START_REF] Bergeot | Passive suppression of helicopter ground resonance using nonlinear energy sinks attached on the helicopter blades[END_REF] and mode coupling instability mitigation in a friction system [START_REF] Bergeot | Mode coupling instability mitigation in friction systems by means of nonlinear energy sinks : numerical highlighting and local stability analysis[END_REF]. As a recent works, aeroelastic suppression of an airfoil with control surface using a nonlinear energy sink is considered in [START_REF] Guo | Aeroelastic suppression of an airfoil with control surface using nonlinear energy sink[END_REF]. Panel flutter suppression with NES is discussed in [START_REF] Pacheco | Panel flutter suppression with nonlinear energy sinks: Numerical modeling and analysis[END_REF]. Note that in these two recent works only complete Limit Cycle Oscillations (LCO) suppression or make the LCOs less intense are considered.

A single NES (coupled to the primary system) was the NES configuration the most often studied. Multi-Degree-Of-Freedom (multi-DOF) NES in series [START_REF] Tsakirtzis | Multi-frequency nonlinear energy transfer from linear oscillators to mdof essentially nonlinear attachments[END_REF][START_REF] Gourdon | Energy pumping for a larger span of energy[END_REF][START_REF] Grinberg | Response regimes in linear oscillator with 2dof nonlinear energy sink under periodic forcing[END_REF][START_REF] Taghipour | Steady state dynamics and robustness of a harmonically excited essentially nonlinear oscillator coupled with a two-DOF nonlinear energy sink[END_REF] or in parallel [START_REF] Vaurigaud | Passive control of aeroelastic instability in a long span bridge model prone to coupled flutter using targeted energy transfer[END_REF][START_REF] Zhang | Impulse-induced vibration suppression of an axially moving beam with parallel nonlinear energy sinks[END_REF][START_REF] Li | Dynamics of two vibro-impact nonlinear energy sinks in parallel under periodic and transient excitations[END_REF][START_REF] Boroson | Optimization under uncertainty of parallel lnonlinear energy sinks[END_REF][START_REF] Bergeot | Asymptotic analysis of passive mitigation of dynamic instability using a nonlinear energy sink network[END_REF] were also considered. The use of several NESs allow to increase the range of excitations (level and frequency) within which the NES can absorb energy efficiently and/or to extract vibration energy simultaneously from multiple linear modes of the primary system.

To capture the complexity of the dynamic of the coupled system (primary system + NES) due to the nonlinear behaviour of the NES, various analysis have been developed including direct time integration method [START_REF] Zhang | Impulse-induced vibration suppression of an axially moving beam with parallel nonlinear energy sinks[END_REF][START_REF] Yao | Using grounded nonlinear energy sinks to suppress lateral vibration in rotor systems[END_REF], harmonic balance method [START_REF] Zang | The evaluation of a nonlinear energy sink absorber based on the transmissibility[END_REF][START_REF] Zhang | An inertial nonlinear energy sink[END_REF] or complex-averaging method [START_REF] Zhang | Nonlinear energy sink with inerter[END_REF]. The complex-averaging method [START_REF] Manevitch | Complex representation of dynamics of coupled nonlinear oscillators[END_REF] gives also access to the slow-flow model of the nonlinear system. This slow-flow is governed by different time scales and analytical methods such as multiple scale analysis [START_REF] Nayfeh | Perturbation Methods[END_REF], or Geometric Singular Perturbation Theory (GSPT) [START_REF] Jones | Geometric singular perturbation theory[END_REF] can be used to extract the complex dynamics. In [START_REF] Starosvetsky | Strongly modulated response in forced 2dof oscillatory system with essential mass and potential asymmetry[END_REF], for example harmonic forcing was considered where response regimes are characterized in terms of periodic and relaxation oscillations of the slow-flow also named Strongly Modulated Responses (SMR). An alternative method known as mixed multiple scale/harmonic balance method has been proposed by Luongo and Zulli [START_REF] Luongo | Dynamic analysis of externally excited nes-controlled systems via a mixed multiple scale/harmonic balance algorithm[END_REF] leading to an equivalent result without the complexification-averaging being necessary. This approach was used by the same authors [START_REF] Luongo | Aeroelastic instability analysis of nes-controlled systems via a mixed multiple scale/harmonic balance method[END_REF] in case of NES control of aeroelastic instability.

To go further, a discussion on the relationship between the dimensionality of the critical manifold, the nature and the distribution of the fixed points of the slow dynamics and the observed response regimes is explored in [START_REF] Gendelman | Targeted energy transfer in systems with external and self-excitation[END_REF] where it is apparent that the analysis is difficult for systems with large dimension. Hence, it seems interesting to investigate the possibility of implementing a procedure able to deal with a multi-DOF system coupled to a large number of NESs.

In this work, auto-oscillating systems are considered and without loss of generality we focus on multidegree-of-freedom unstable vibrating systems undergoing cubic nonlinearities and having only one mode which can become unstable through Hopf bifurcation. Our objective is to predict the steady-state regimes (i.e complete suppression of instability, mitigation through periodic response or mitigation through SMR) when the (primary) system is coupled to a family of single NESs. The proposed procedure starts with the diagonalization of the primary underlying linear system written in the state-space form using the so-called biorthogonal transformation. Then, the dynamics of the diagonalized system is reduced ignoring the stable modes and keeping only the unstable mode of the primary system. The complexification-averaging method is next applied to the resulting coupled system, leading to slow-flow of the system. Based on the GSPT, it is shown that the critical manifold of the system can be reduced to a one dimensional parametric curve evolving in a multidimensional space. This form of critical manifold is similar to that obtained considering a network of parallel single NESs [START_REF] Bergeot | Asymptotic analysis of passive mitigation of dynamic instability using a nonlinear energy sink network[END_REF]. Knowing the stability properties of the critical manifold and the fixed point of the slow-flow (position and stability), the response regimes can be predicted. Finally as in [START_REF] Bergeot | Asymptotic analysis of passive mitigation of dynamic instability using a nonlinear energy sink network[END_REF] we introduce and predict the mitigation limit of the set of NESs, defined as the value of a chosen bifurcation parameter which separates harmless situations (corresponding to responses in which the set of NESs acts) from harmful situations (corresponding to responses in which the set of NESs does not act).

The system under study is introduced in the next section starting from the physical model, applying successively a rescaling step, a diagonalization step and a reduction step. In Sect. 3, the slow-flow dynamics is first obtained applying the complexification-averaging method and then analyzed using the geometric singular perturbation theory. Sects. 3.1 to 3.4 describe in detail the mathematical developments which allow to obtain the critical manifold and the fixed points (position and stability) of the slow-flow. From the knowledge of that, the Sect. 3.5 provides the procedure to predict the steady-state response regimes and consequently the mitigation limit. Finally, in Sect. [START_REF] Lu | Nonlinear dissipative devices in structural vibration control: A review[END_REF], two examples are considered. In Sect. 4.1, the method is first applied to the prediction of the mitigation limit of a breaking system subject to friction-induced vibrations coupled to two NESs. A airfoil model undergoing an aeroelastic instability coupled to a set of NESs is next considered in Sect. 4.2.

The system under study

The initial model

One considers in this work the following primary nonlinear system

̃ ̈ + ̃ ̇̃ + ̃ ̃ + ̃ (̃ ) = , ( 1
)
where ̃ = ̃ 1 , … , ̃ with (.) denotes the transpose operator, the dot represents time-differentiation and

̃ = ⎡ ⎢ ⎢ ⎣ ̃ 1,1 … ̃ 1, ⋮ ⋱ ⋮ ̃ ,1 … ̃ , ⎤ ⎥ ⎥ ⎦ , ̃ = ⎡ ⎢ ⎢ ⎣ ̃ 1,1 … ̃ 1, ⋮ ⋱ ⋮ ̃ ,1 … ̃ , ⎤ ⎥ ⎥ ⎦ and ̃ = ⎡ ⎢ ⎢ ⎣ ̃ 1,1 … ̃ 1, ⋮ ⋱ ⋮ ̃ ,1 … ̃ , ⎤ ⎥ ⎥ ⎦ (2a)
are constant matrices. The nonlinear vector-valued function ̃ is called nonlinearity of the primary system. We assume that each component of ̃ is a linear combination of monomial terms of order 3.

We also assume that Eq. ( 1) can undergo a single dynamic instability of the trivial solution and that the nonlinear function ̃ allows the existence of Limit Cycle Oscillations (LCOs) on which the system can saturate.

System (1) can undergo a single dynamic instability of the trivial solution. In order to protect the primary system from this instability, purely cubic ungrounded NESs with masses ̃ ℎ , damping coefficients ̃ ℎ and cubic stiffness ̃ ℎ ( = 1, … , ) are used. Taking into account the NESs displacements h ( ) ( = 1, … , ), the equations of the whole coupled system result on the following form

̃ ̈ + ̃ ̇̃ + ̃ ̃ + ̃ (̃ ) + ̃ ( ̃ ℎ ) ̇̃ -̇̃ + ( ̃ ℎ ) ̃ -̃ = (3a) ( ̃ ℎ ) ̈ - ( ̃ ℎ ) ̇̃ -̇̃ + ( ̃ ℎ ) ̃ -̃ = , ( 3b
)
where ̃ = ( h1 , … , h ) , and the constant matrices

̃ = ⎡ ⎢ ⎢ ⎣ ̃ 1,1 … ̃ 1, ⋮ ⋱ ⋮ ̃ ,1 … ̃ , ⎤ ⎥ ⎥ ⎦ and = ⎡ ⎢ ⎢ ⎣ 1,1 … 1, ⋮ ⋱ ⋮ ,1 … , ⎤ ⎥ ⎥ ⎦ (4)
are the influence coefficient matrices (depending on the position of the NESs). The nonlinear vector function , due to the NESs attachments, is given by

( ) = 3 1 , … , 3 . ( 5
)
Finally ( ) denotes the diagonal matrix with the main diagonal ( 1 , 2 , ⋯ , ). System (3) is the system under study (a multi-DOFs the primary system coupled with one-DOF NESs).

We now assume that the masses, ̃ ℎ , and the damping coefficients, ̃ ℎ , of the NESs are small and we introduce a small, dimensionless parameter (0 < ≪ 1) and the associated rescaled coefficients, ℎ and ℎ , as

̃ ℎ = ℎ and ̃ ℎ = ℎ , for = 1, … , . (6)
Then, rescaling the variables ̃ and ̃ through as

= ̃ √ and = ̃ √ (7)
and inserting the rescaled parameters ℎ and ℎ and the rescaled variables and into (3) lead to

̃ ̈ + ̃ ̇ + ̃ + ̃ ( ) + ̃ ( ) ̈ = (8a) ( ℎ ) ̈ - ( ℎ )( ̇ -̇ ) + ( ̃ ℎ ) ( -) = . ( 8b
)
Finally multiplying (8a) by ̃ -and (8b) by (1∕ ℎ ) we obtain

̈ + ̇ + + ( ) + ̈ = (9a) ̈ - ( ℎ )( ̇ -̇ ) + ( ℎ ) ( -) = , (9b) with = ̃ -̃ , = ̃ -̃ , = ̃ -̃ ( ), ℎ = ℎ ∕ ℎ , ℎ = ̃ ℎ ∕ ℎ and ( ) = ̃ -̃ ( ).

Model reduction

In this section, the system ( 9) is reduced taking into account that the primary system (1) can undergo only a single instability of the trivial solution.

First, to simplify asymptotic analysis in next section, it is convenient to introduce new coordinates as

= + (10a) = -, (10b)
giving reciprocally,

= ( + ) -1 ( - ) (11a) = ( + ) -1 ( + ), (11b)
where is the identity matrix of size . Using Eqs. [START_REF] Gendelman | Bifurcations of self-excitation regimes in a Van der Pol oscillator with a nonlinear energy sink[END_REF] with the inverse matrices ( + ) -1 and ( + ) -1 , approximated respectively by ( -) and ( -), Eq. ( 9) is transformed into the following form

̈ + ̇ + - ( ̇ + ̇ ) + ( + ) + ( ) = (12a) ̈ + ( ℎ ) ̇ - ̇ - + ( ℎ ) ( )+ ( ℎ ) ̇ + ( ̇ + ̇ ) + ( + ) + ( ℎ ) ( ) = , ( 12b
)
where only the first order term in has been retained in the nonlinearity of the primary underline linear system.

To capture the essential features of the single instability, model reduction is only performed on Eq. (12a), keeping the whole second equation. To achieve this, Eq. (12a), is written in state-space form as follows

̇ = + + + ̇ + ( ) (13)
where

= ( 1 , … , , ̇ 1 , … , ̇ ) , trivially redefined as ( ) = ( ), = - - , (14) 1
= , 2 = , 3 = and 4 = - . ( 15
)
The matrix characterizes the linear dynamics of the primary system. Because the matrix is not a symmetric matrix, its diagonalization must necessarily be carried out using biorthogonality property of the right eigenvectors ( = 1, … , 2 ) and the left eigenvectors ( = 1, … , 2 ) of hereafter defined. The method is briefly recalled in this section. The general procedure is provided for example in [START_REF] Meirovitch | Fundamentals of Vibrations[END_REF].

Let us consider the right and left eigenvector matrices

= * … * and = * … * (16)
and the associated eigenvalue diagonal matrix = ( 1 , * 1 , ⋯ , , * ) corresponding to the following right and left eigenvalue problems

= and = . ( 17
)
The right and left eigenvectors satisfy the biorthogonality properties stating that and are diagonal matrices. It is convenient to normalize the two sets of eigenvectors in order to obtain = , [START_REF] Bergeot | Mode coupling instability mitigation in friction systems by means of nonlinear energy sinks : numerical highlighting and local stability analysis[END_REF] that results in = .

(

) 19
The biorthogonal transformation consists in introducing the biorthogonal coordinates which are constituted of pairs of complex conjugates, and * ( = 1, … , ), and defined by the following relations

= ⇔ = , ( 20
)
where

= ( 1 , * 1 , … , , * ) . ( 21
)
Introducing Eq. [START_REF] Pacheco | Panel flutter suppression with nonlinear energy sinks: Numerical modeling and analysis[END_REF] in Eq. ( 13), the equations of motion take the form of the following system

̇ = + 1 + 2 + 3 ̇ + 4 ( ), (22)
and therefore Eqs. [START_REF] Lee | Suppression aeroelastic instability using broadband passive targeted energy transfers, part 1: Theory[END_REF] are equivalent to

̇ - - 1 + 2 + 3 ̇ - 4 ( ) = 0 (23a) ̈ + ( ℎ ) ̇ - ( + ) - ( + )+ ( ℎ ) ( ) + ( ℎ ) ̇ + ( + ) + ̇ + ( + ) + + ( ℎ ) ( ) = , ( 23b
)
where the matrix (respectively the vector ) has been split into a × -block matrix (respectively × 1block vector) as

= and = ( , ) (24)
respectively.

From now we assume that one (and only one) mode (without loss of generality the first mode) of the primary system become unstable through Hopf bifurcation, when the chosen bifurcation parameter, denoted , crosses a particular parameter value 0 called bifurcation point, i.e

• if < 0 , Re < 0 ∀ ∈ [1,
] and the trivial solution of the primary system is stable,

• if > 0 , Re 1 > 0, Re < 0 ∀ ∈ [2,
] and the trivial solution of the primary system is unstable.

Because we assume that only the first pair of eigenvalues ( 1 , * 1 ) may have positive real part, after some exponentially decaying transients the components and * ( = 2, … , ) become small. Therefore, all terms related to and * ( = 2, … , ) are omitted from further consideration and Eq. ( 22) is reduced to

̇ 1 = 1 1 + 1 1 ( 1 1 + * 1 * 1 ) + 2 + 3 ̇ + 1 4 1 1 + * 1 * 1 , (25)
in which we stated that = * = 0 ( = 2, … , ). Finally, grouping Eqs. ( 25) and (12b), we obtain the reduced system

̇ 1 -( + ) 1 -1 1 ( 1 1 + * 1 * 1 ) + 2 + 3 ̇ -1 4 1 1 + * 1 * 1 = 0 (26a) ̈ + ( ℎ ) ̇ - ( 1 1 + * 1 * 1 ) - ( 1 1 + * 1 * 1 ) + ( ℎ ) ( )+ ( ℎ ) ̇ + ( 1 1 + * 1 * 1 ) + ̇ + ( 1 1 + * 1 * 1 ) + + ( ℎ ) ( ) = (26b)
where the vector 1 has been split as 1 = ( 1 , 1 ) . Note that 1 has also been replaced by

Asymptotic analysis

We focus on the motion of the system in the vicinity of a 1 ∶ 1 resonance associated with a frequency close to the natural frequency ( see Eq. ( 27)) of the primary system.

The slow-flow

The slow-flow describes the evolution of the amplitude and phase of the system which evolve slowly compare to the oscillations at the frequency . To obtain it, we use the Complexification-Averaging method introduced by [START_REF] Manevitch | Complex representation of dynamics of coupled nonlinear oscillators[END_REF] and discussed in detail in [START_REF] Vakatis | Nonlinear Targeted Energy Transfer in Mechanical and Structural Systems[END_REF].

The first step is the complexification1 which consists in introducing the following complex vector

= ̇ + (28)
or equivalently

= - * 2 , ̇ = + * 2 and ̈ = ̇ - 2 + * (29)
and expressing the complex variables 1 and as

1 = and = ( 30
)
where the scalar and the vector are the complex (assumed slow modulated) amplitude of the fast component . The second step is the averaging which consists in, first, substituting Eqs. ( 29) and ( 30) into Eqs. [START_REF] Li | Dynamics of two vibro-impact nonlinear energy sinks in parallel under periodic and transient excitations[END_REF], and next, averaging the resulting equations over one period = 2 ∕ . We obtain the slow-flow described by the following differential equations

̇ = ( , ) , (31a) ̇ = ( , , ) (31b)
where the functions and are defined by

( , ) = + 1 1 1 + 1 4 ( )| | 2 + 1 2 1 2 + 1 2 1 3 (32a) = 0 + 1 | | 2 + 1 , (32b)
where the vector ( ) is obtained from the nonlinear function

( ) replacing each monomial 1 1 1 by * 1 1 1 + 1 * 1 1 + 1 1 *
1 (where 1 is the -th component of ) and

( , , ) = ( 1 + 1 ) - 1 2 + ℎ - 3 ℎ 4 3 | | 2 - 1 + 1 1 1 + 1 + 1 2 1 (33a) + 1 1 2 1 2 + 1 2 1 3 + 1 2 + 1 + ℎ + 3 ℎ 4 3 | | 2 , ( 33b
) = 0 ( , ) + 1 ( , ) . ( 33c
)
The scalars 0 and 1 , the vector 1 and the vector functions 0 and 1 have been introduced to simplify the analysis in next Section. Finally, it is convenient to rewrite Eqs. [START_REF] Nayfeh | Perturbation Methods[END_REF] considering real variables. To achieve this, polar coordinates are introduced as

= , (34a) = , ( 34b
)
where = 1 , … , and and (respectively and ) characterized the modulus and the argument of (respectively ). Substituting Eqs. [START_REF] Luongo | Dynamic analysis of externally excited nes-controlled systems via a mixed multiple scale/harmonic balance algorithm[END_REF] into Eqs. [START_REF] Nayfeh | Perturbation Methods[END_REF] and separating in real and imaginary parts, Eqs. [START_REF] Nayfeh | Perturbation Methods[END_REF] reduce to

̇ =  ( , , ) (35a)
̇ =  ( , , , ) (35b) ̇ =  ( , , , ) , ( 35c
)
where the argument differences = -have been considered and written in vector form as = 1 , … , . The functions  ,  and  are deduced from functions and in Eqs. [START_REF] Nayfeh | Perturbation Methods[END_REF] as

 = Re - , (36a)
 = Re - , (36b)
 = Im - -Im - (36c)
where = (1, … , 1) is a vector of size .

The Critical Manifold

The slow-flow described by Eqs. ( 35) has a slow-fast nature. Indeed, Eq. (35a) contains only ( ) terms and Eqs. (35b) and (35c) contain both (1) and ( ) terms. Consequently is the slow variable whereas the vectors and contain the fast variables and ( = 1, … , ). Equations ( 35) can be reformulated by switching from the fast time scale to the slow time scale = as

′ =  ( , , ) (37a)
′ =  ( , , , ) (37b)
′ =  ( , , , ) , (37c)
where ′ = . Solutions of the slow-fast system (35) (or ( 37)) can exhibit slow and fast epochs characterized by the speed at which the solution advances. Stating = 0, the following subsystems are derived from ( 35) and ( 37) respectively

̇ = 0 (38a) ̇ =  ( , , , 0) , (38b) ̇ =  ( , , , 0) , ( 38c
)
which is the fast subsystem, and

′ =  ( , , ) (39a)
0 =  ( , , , 0) , (39b) 0 =  ( , , , 0) , ( 39c
)
which is the slow subsystem.

The algebraic equations (39b) and (39c) of the slow subsystem define the so-called Critical Manifold [START_REF] Jones | Geometric singular perturbation theory[END_REF] as

∶= ( , , ) ∈ ℝ 2 +1 | | |  ( , , , 0) = 0 and  ( , , , 0) = 0 (40)
corresponding to the subspace where the slow motions take place and also to the fixed points of the fast subsystem [START_REF] Fenichel | Geometric singular perturbation theory for ordinary differential equations[END_REF]. Looking at Eqs. (33c) and ( 36), one can deduce that

 ( , , , 0) = Re ( -) 0 and  ( , , , 0) = Im ( - ) 0 , ( 41
)
and after some algebraic manipulation, Eqs (39b) and (39c) can be written in the reduced complex form as

- 1 + 1 - 1 2 + ℎ - 3 ℎ 4 3 2 = 0. (42)
Introducing the complex functions ( = 1, … , ) of a real variable as

( ) = 1 - ℎ - 3 ℎ 4 4 2 = ( ) + ( ), (43) where
= -2 1 + 1 (44)
with (.) stands for the -th coordinate of the vector (.), the Critical Manifold is characterized by

2 = ( ), ( = 1, … , ), (45a)
= -arg ( ) , ( = 1, … , ) (45b)
where the real functions ( = 1, … , ) of a real variable are defined as

( ) = 2 ( ) 2 + ( ) 2 . ( 46
)
It is interesting to note that the general form (45) of the Critical Manifold is the same as found in [START_REF] Bergeot | Asymptotic analysis of passive mitigation of dynamic instability using a nonlinear energy sink network[END_REF] in the case of a one-DOF unstable primary system coupled to several parallel NESs. Indeed, in the latter case, the equivalent functions can be obtained from (43) stating = 1. From this form, as suggested in [START_REF] Bergeot | Asymptotic analysis of passive mitigation of dynamic instability using a nonlinear energy sink network[END_REF], specific properties of the Critical Manifold can be deduced which can be advantageously used to characterize the response regimes of the system under study.

Critical Manifold properties and stability

Due to Eqs. (45), the critical manifold appears as a one dimensional parametric curve evolving in ℝ + +1 as a solution of the nonlinear equations (45a) with respect to the +1 unknown ( 1 , … , , ). A classical continuation method can be used to obtain . Equivalently can be obtained solving the following differential-algebraic equations

∑ =1 ( ) 2 + ( ) 2 = 1, for 1 < < 2 (47a) ( ) 2 - ( ( )) = 0, ( = 1, … , ) (47b)
(0) = 1, (0) = 0, ( = 1, … , ) (47c)
(0) = 0 , (0) = 0 , ( = 1, … , ) (47d)
where is arclength variable in the range 1 (< 0) to 2 (> 0) and ( 0 1 , … , 0 , 0 ) ∈ ℝ + +1 is a root of the algebraic equations (47b), chosen to be close to trial solution . An example of typical critical manifold for a system containing two NESs are plotted Fig. 1. starts at the origin by a straight line followed by a finite number of successive simple curved lines curving in opposite directions (with respect to ) and ended by an unbounded straight line. Each simple curve is characterized by a curving point defined by a horizontal tangency in at least one ( , )-plane (i.e ∕ = 0). The critical manifold can also be viewed as the intersection of the manifolds

= {( , ) ∈ ℝ + +1 | | | 2 = ( )}, ( = 1, ⋯ , ) (48)
in the ( 1 , … , , )-space.

A typical projection of onto the subspace ℝ + 2 corresponding to the variables ( , ) is depicted in Fig. 2. Its geometrical form is characterized by the properties of the function . From Eq. ( 46), it can be shown that, for each , if the following relation holds

ℎ < √ 3 , ( 49
)
then admits a local maxima max and a local minima min given as the positive roots of its derivative ∕ . An easy calculus shows that

max = 2 3 2 √ √ √ √ √ √ 2 - √ 1 -3 ℎ 2 ℎ and min = 2 3 2 √ √ √ √ √ √ 2 + √ 1 -3 ℎ 2 ℎ ( 50
)
with max < min . Each extrema, defined by ( p , p ) with p = √ ( p ) where the upper-script p stands for min or max (see Fig. 2 where condition (49) is satisfied) is associated to a curving point of .

As already mentioned, each point of is also a fixed point for the fast subsystem [START_REF] Fenichel | Geometric singular perturbation theory for ordinary differential equations[END_REF] and the stability can be determined by examining the sign of the eigenvalue real parts of the Jacobian matrix of the differential system (38b)(38c). Hence the stability range of is defined by the points 1 , … , , on that satisfy Stability condition can be also written equivalently in terms of arclength variable at the point 1 ( ) , … , ( ) , ( ) of as

∀ ∈ [1, ], > 0. (51)
∀ ∈ [1, ], ( ) ( ) -1 > 0. (52)
We assume from now that Eq. ( 49) is satisfied for the NESs. For each (1 ≤ ≤ ), the extrema points ( max , max ) and ( min , min ) characterize two bounds on where ceases to be hyperbolic2 connecting stable or attractive (continuous line) and unstable or repulsive (dashed line) parts of (see Fig. 2). It follows that ( max , max ) (respectively ( min , min )) defines a point ma, (respectively mi, ) on called fold point by ma,m = ( ma,

1 , … , ma, -1 , max , ma, +1 , … , ma, , max ) (53a) (respectively mi, = ( mi, 1 , … , mi, -1 , min , min, +1 , … , min, , min )) (53b)
with ma, (respectively mi, ) satisfied for = 1, … , , ≠ , ( ma, ) = ( max ) (respectively ( mi, ) = ( min )). The number of fold points is equal to 2 (2 by NES).

Moreover in each ( , )-plane, the dynamics of any point not on is entirely defined by the direction of the fast variable as given by Eq. (38a) (i.e. = Cte) and which is indicated by the horizontal arrows in Fig. 2. Hence at a fold point ma, (respectively mi, ) on , a jump occurs (fast dynamic) to a point on named arrival point and denoted ma, (respectively mi, ). During the jump, = max and for = 1, … , , ≠ the trajectory lies in the manifold leading to ma,

Stable Unstable • • • • • • • • 0.0 0.
= ( ma, 1 , … , ma, -1 , up , ma, +1 , … , ma, , max ) (54a)
(respectively mi, = ( mi,

1 , … , mi, -1 , do , mi, +1 , … , mi, , min )) (54b)
with the ordinate up (respectively do ) is obtained solving up = (respectively do = ) that leads to

up = 2 √ 2 3 2 √ √ √ √ √ √ 1 + √ 1 + 3 ℎ 2 ℎ (respectively do = 2 √ 2 3 2 √ √ √ √ √ √ 1 - √ 1 -3 ℎ 2 ℎ ). ( 55
)
It follows that the stability domain of is only characterized in the ( 1 , 2 , … , )-space (ℝ + ) by the subspace given in terms of Cartesian product as

= ∏ =1 , ( 56
) where = 0 max ∪ min + ∞ (57)
denotes the stability range associated to the -th NES (see Fig. 2).

Eq. ( 56) may be expanded as the union of 2 disjunct Cartesian products (with = 1, … , 2 ) as

= 2 ⋃ =1 , ( 58
)
where appears as the product of close or left-close intervals resulting in an isolated stable or attractive part of the the critical manifold as

= ∏ ∈ [0, max ] ∏ ∈ [ min , +∞), (59)
with is a part of the finite set {1, 2, … , } and its complementary. The first and the last Cartesian products are chosen as 1 = ∏ =1 [0, max ] and 2 = ∏ =1 [ min , +∞)] respectively. With this choice, the first branch of is defined on 1 whereas the last branch of is defined on 2 (see Fig. 1b).

Finally a more compact representation of the CM can be used obtained projecting in the ( 1 , 2 , … , )space resulting in a one dimensional parametric curve in ℝ + . The intersection of this curve with defines the stable zone of the CM which can be decomposed as successive branches, one branch by subset named hereafter . The first branch 1 starts at origin of ℝ + denotes 1 1 = and ends at the fold point 2 1 = ma, ̄ where ̄ = argmin 1≤ ≤ max (i.e ̄ is the subscript value at which max is minimal). The last branch 2 starts at the fold point 1 2 = mi, ̄ where ̄ = argmax 1≤ ≤ min (i.e ̄ is the subscript value at which min is maximal) and it is unbounded (i.e 2 2 = ∞). The other branches are ended by two fold points defined by 1 = mi, ̄ where ̄ = argmax ∈ min and 2 = ma, ̄ where ̄ = argmin ∈ max . Note that for each of the points 1 and 2 (except for 1 1 and 2 2 ) we can also associated the arrival points (following the jump) named 1 and 2 defined by 1 = ma, ̄ and 2 = mi, ̄ with the corresponding value of ̄ .

Such a multi-S-shape of the critical manifold provides a possibility for relaxation oscillations [START_REF] Grasman | Asymptotic Methods for Relaxation Oscillations and Applications[END_REF] of the slow-flow characterized by fast transitions (jumps) of the dynamics during each cycle (the possible fast transitions between the two stable branches of are denoted by arrows in Fig. 2). Because several NESs are considered, complex relaxation oscillations scenarios may be contemplated. Such relaxation oscillations of the slow-flow explain the existence of Strongly Modulated Responses [START_REF] Starosvetsky | Strongly modulated response in forced 2dof oscillatory system with essential mass and potential asymmetry[END_REF][START_REF] Gendelman | Asymptotic analysis of passive nonlinear suppression of aeroelastic instabilities of a rigid wing in subsonic flow[END_REF][START_REF] Gendelman | Bifurcations of self-excitation regimes in a Van der Pol oscillator with a nonlinear energy sink[END_REF] (SMRs) of the nonaveraged system [START_REF] Li | Dynamics of two vibro-impact nonlinear energy sinks in parallel under periodic and transient excitations[END_REF].

We proceed to a more detailed analysis of the possible steady-state regimes in next section computing fixed points of the slow-flow and their stability.

The fixed points of the slow-flow

Following the Geometric Singular Perturbation Theory (GSPT) [START_REF] Fenichel | Geometric singular perturbation theory for ordinary differential equations[END_REF][START_REF] Jones | Geometric singular perturbation theory[END_REF][START_REF] Desroches | Mixed-mode oscillations with multiple time scales[END_REF], the nontrivial fixed points of the slow-flow (31) can be approximated computing the fixed points of Eq. [START_REF] Desroches | Mixed-mode oscillations with multiple time scales[END_REF].

Substituting first Eq. (45a) with = 1 (any ∈ [1, ] can be chosen) into Eq. (39a) , we obtain

1 2 √ 1 1 1 1 ′ 1 =  √ 1 1 , 1 , … , , 1 , … , (60)
which can be rewritten, substituting Eq. (45b) into the resulting equation, as

1 1 ′ 1 = 1 1 , … , , (61)
( ) + ∑ =2 ( ) 1 (69a) = 1 ( ) + ∑ =2 ( ) 1 1 . ( 69b
)
Moreover, to be stable a fixed point must be on a stable branch of the critical manifold . Finally, a fixed point * = * 1 , … , * of Eq. ( 61) is stable if the two following conditions are satisfied 1.

1 | | | | = * < 0, 2. ∀ , 1 ≤ ≤ , * ∈ [0 ] ∪ [ + ∞).
A stable (respectively unstable) fixed point is denoted * = ( * ,1 , … , * , ) (respectively * = ( * ,1 , … , * , )). Finally, each fixed point (stable or unstable) as a point of is a part of one of branches .

Prediction of the steady-state regimes

As discussed for example in [START_REF] Lee | Suppression aeroelastic instability using broadband passive targeted energy transfers, part 1: Theory[END_REF][START_REF] Gendelman | Bifurcations of self-excitation regimes in a Van der Pol oscillator with a nonlinear energy sink[END_REF][START_REF] Gendelman | Asymptotic analysis of passive nonlinear suppression of aeroelastic instabilities of a rigid wing in subsonic flow[END_REF][START_REF] Bergeot | Mode coupling instability mitigation in friction systems by means of nonlinear energy sinks : numerical highlighting and local stability analysis[END_REF][START_REF] Bergeot | Asymptotic analysis of passive mitigation of dynamic instability using a nonlinear energy sink network[END_REF] when a NES is coupled to an unstable system, four types of steady-state regimes can be generated corresponding to a complete suppression of the instability, mitigation through Periodic Response (PR), mitigation through Strongly Modulated Response (SMR) or no mitigation. Complete suppression means that the trivial fixed point, common to both the non-averaged system and the slow-flow, is stable. Mitigation through PR means that a stable (nontrivial) fixed point of the slow-flow is reached leading, for non-averaged system, to a LCO with an amplitude smaller than the amplitude of the corresponding LCO undergone by the primary system without absorbers. Mitigation through SMR means that the slow-flow undergoes relaxation oscillations corresponding to a phase-amplitude modulated regime for the non-averaged system, the so-called "Strongly modulated response". Finally no mitigation means that the NESs are not able to mitigate the instability and the non-averaged system saturates on a LCO which has an amplitude close to that of the case without NES. Regarding the slow-flow, it reaches a stable fixed point with large amplitude with respect to the variable . This stable fixed point denoted * corresponds to the stable fixed point at which the amplitude with respect to the variable is maximized. This stable fixed point is positioned in the branch 2 .

As proposed in [START_REF] Bergeot | Asymptotic analysis of passive mitigation of dynamic instability using a nonlinear energy sink network[END_REF], these regimes can also be classified into two groups named harmless situation and harmful situation reflecting the situations in which the NESs act or not. Harmless situation includes complete suppression, mitigation through PR and mitigation through SMR regimes. Harmful situation corresponds to no mitigation regime. The main objective of NESs is to put the system in harmless situations.

The nature of the steady-state regime depends on some features of the slow-flow. Firstly, the initial conditions from which one can know where the fast dynamics leads the system trajectory on the critical manifold. Secondly, the slow subsystem [START_REF] Desroches | Mixed-mode oscillations with multiple time scales[END_REF], in particular its fixed points (position and stability), which guides the dynamics on at the slow time scale. In accordance with real word situations, the following prediction is limited to a set of initial conditions (0), 1 (0), … , (0) as a small perturbation of the trivial solution.

To predict the steady-state regimes, we have the CM , the stable domain of in terms of the Cartesian products and the associated branches in the ( 1 , … , )-space with the corresponding bounds 1 and 2 (which are fold points) and the associated arrival points 1 and 2 . We have also the stable fixed points * and unstable fixed points * which are positioned in the corresponding branches . Note that the bound point 1 1 is the trivial equilibrium of the system. Recalling the two subsystems [START_REF] Fenichel | Geometric singular perturbation theory for ordinary differential equations[END_REF] and [START_REF] Desroches | Mixed-mode oscillations with multiple time scales[END_REF], the trajectory of the system is a succession of slow (given by [START_REF] Fenichel | Geometric singular perturbation theory for ordinary differential equations[END_REF]) and fast (given by ( 39)) parts. The slow parts lead in whereas the fast parts are jump occurring out of between two points of .

From a given initial condition, after a first fast transient response (a jump), the trajectory reaches the critical manifold on the branch 1 = 1 at the point 1 .

We assume that we are at the end of the step of the procedure i.e on the branch = at the (arrival) point .

The step + 1 (from to + 1) reads as follows.

The trajectory follows = slowly starting form and the following cases are possibles:

Case 1. There are no fixed points in , then, depending of the sign of the function

1 , … , | | | = (see
Eq. ( 67)), 1 or 2 are reached, the trajectory leaves jumping to the arrival point +1 = 1 (or 2 ) on an other branch ̄ given +1 = ̄ and +1 = ̄ .

Case 2. The arrival point is between an unstable fixed point and 1 (resp. 2 ), then 1 (resp. 2 ) is reached, the trajectory leaves jumping to the arrival point +1 = 1 (resp. 2 ) on an other branch ̄ given +1 = ̄ and +1 = ̄ .

Case 3. The arrival point is between an stable fixed point and 1 (resp. 2 ), then the stable fixed point is reached and the procedure stops.

Case 4. The arrival point is between an stable fixed point and unstable fixed point, then the stable fixed point is reached and the procedure stops.

At the end of step + 1, if the trajectory reaches a new branch +1 at the point +1 a new step can be performed.

The procedure stops when one of the following conditions are meet resulting in harmless or harmful situations:

• The system is in a harmless situation if

• During a step, the trajectory stops on a stable fixed point which differs from * (∈ 2 ) resulting to a Mitigation through Periodic Response (Cases 3 or 4).

• At each step, the trajectory reaches a new branch. Due to the finite number of branch , the trajectory will reach the same branch in the same condition resulting to a Mitigation through Strongly Modulated Response (SMR) (Cases 1 or 2).

• The system is in a harmful situation if the trajectory stops on the stable fixed point * (∈ 2 ). Indeed, in this case, the reached fixed point has a large amplitude close to that of the case without NES.

Using the notion of harmless and harmful situation, the concept of mitigation limit of the system can be introduced. For a given parameter of the system named "bifurcation parameter", the mitigation limit is defined as the value of a chosen bifurcation parameter which separates harmful situation from harmless situation.

From a given set of parameters, the procedure described hereinbefore gives a theoretical prediction of the resulting steady-state regime and therefore allows to know if the system is a harmless or harmful situation. Consequently, the mitigation limit can be predicted theoretically as the first value of the chosen bifurcation parameter for which harmfull situation holds.

Applications

Two self-excited mechanical systems are used in this section to check the relevance of the analytical procedure proposed in Sect. 3.5. The first one is the well-known Hultén's model [START_REF] Hultén | Friction phenomena related to drum brake squeal instabilities[END_REF] which reproduces the typical dynamic behavior of friction systems. The second one is a 2-DOF nonlinear airfoil model. It has been already used by Lee et al in [START_REF] Lee | Suppression aeroelastic instability using broadband passive targeted energy transfers, part 1: Theory[END_REF].

Mode coupling instability mitigation in a friction system using two NESs

The system is shown in Fig. 3. It is composed of a 2-DOFs Hultèn's model (the primary system) coupled to two NESs in ungrounded configuration. This model has already been considered in [START_REF] Bergeot | Mode coupling instability mitigation in friction systems by means of nonlinear energy sinks : numerical highlighting and local stability analysis[END_REF] where it was shown that the use of NES appears to be an interesting way to control mode-coupling instability in braking systems. The objective here is to verify that the method presented here allows us to distinguish between strongly modulated responses and no mitigation responses, the method in [START_REF] Bergeot | Mode coupling instability mitigation in friction systems by means of nonlinear energy sinks : numerical highlighting and local stability analysis[END_REF] being not able to make this distinction in all the cases. Following Sect. 2.1, the equations of motion of the system take the form (see Appendix A)

′′ 1 + 1 ′ 1 + 1 -Ω 2 2 + 1 3 1 - 2 3 2 + ℎ 1 ℎ ′′ 1 = 0 (70a) ′′ 2 + 2 ′ 1 + Ω 2 2 + 1 + 1 3 1 + 2 3 2 + ℎ 2 ℎ ′′ 2 = 0 (70b) ℎ ′′ -ℎ ′ -ℎ ′ -ℎ -ℎ 3 = 0 ( = 1, 2) (70c)
where Eqs. (70) correspond to Eqs. [START_REF] Yao | Using grounded nonlinear energy sinks to suppress lateral vibration in rotor systems[END_REF]. The corresponding Eqs. [START_REF] Lee | Suppression aeroelastic instability using broadband passive targeted energy transfers, part 1: Theory[END_REF] and Eqs. [START_REF] Grinberg | Response regimes in linear oscillator with 2dof nonlinear energy sink under periodic forcing[END_REF] hereafter referred indistinctly as the Reference Model can be easily deduced as well as the slow-flow model corresponding to Eqs. [START_REF] Nayfeh | Perturbation Methods[END_REF].

The bifurcation parameter under consideration throughout this section is the friction coefficient . The following set of numerical values are used, for the other parameters:

1 = 5, 2 = 0, 1 = 0.01, 2 = 0.08 and Ω = 0.9 for the primary system, and ℎ 1 = 0.75, ℎ 2 = 0.2, ℎ 1 = 4 and ℎ 2 = 9 for the NESs. Moreover, ℎ 1 = ℎ 2 = 2 and several values for will be considered.

This set of numerical values leads to a Hopf bifurcation for the Hultèn's model without NESs at wo Hopf = 0.073. If < wo Hopf the system is stable and if > wo Hopf the system is unstable. In the latter case, the model variables increase and reach a LCO due to the nonlinear terms. The LCOs are what we call "harmful situations", i.e. the large amplitude vibrations we want to mitigate using the two NESs. Here we focus on the possible mitigation regimes which can appear when the primary system is unstable.

First of all, the procedure described in Sect. 3.5 for the prediction of the trajectory of the slow-flow in the ( 1 , 2 )-space is illustrated in Fig. 4 for = 2 ⋅ 10 -3 and four values of the bifurcation parameter (friction coefficient). The CM , computed from (47) and (52), is depicted with their stable parts (in gray), , i.e. the (with = 1, … , 2 = 4), and unstable parts (in black). The stable and unstable fixed points, computed from (63) and (69), are represented by magenta and blue points respectively. We can see in Figs. 4b to 4d a stable fixed point with a large value on the unbounded straight line (i.e. 2 = 4 ). The latter corresponds to a large amplitude LCO mentioned above, if it is reached, we are in a harmful situations.

Following the procedure (as described in Sect. 3.5) the theoretical trajectory of the slow-flow is determined and plotted with a dashed green line. This theoretical trajectory is compared to the trajectory obtained 40) together with theoretical and numerical trajectories of the slow-flow [START_REF] Nayfeh | Perturbation Methods[END_REF]. Selected parameter values with = 2 ⋅ 10 -3 . from the direct numerical integration of the slow-flow (plotted in red). The Fig. 4a shows an harmless situation through a periodic regime, i.e. a stable fixed point of the slow-flow (in magenta) which is not in 4 (here it is in 1 ) is reached. In this case, the theoretical and numerical trajectories of the slow are almost superimposed. In Figs. 4b and4c a stable fixed point is reached in 4 , the Fig. 4d corresponds to an harmful situation. One can see that for the last three situations (i.e Figs. 4b, 4c and4d), during the slow epochs on the theoretical and numerical trajectories of the slow-flow are again almost superimposed. However, during fast epochs (i.e. the jumps between each slow epoch) a difference is observed. In particular, unlike the theoretical trajectory, the numerical trajectory is not at each jump parallel to only one direction of the ( 1 , 2 )-space. This difference is due to the 0-order approximation of the proposed asymptotic approach. Indeed, figures similar to Fig. 4 using = 10 -4 would show a perfect agreement between theoretical and numerical trajectories.

• µ • µ c x1 c x2 k x2 , k N L x2 k x1 , k N L x1 x2 x1 m NES1 NES2 (a) h1 NES 1 mh1 ch1 k h1 , k N L h1 (b) h2 NES 2 mh2 ch2 k h2 , k N L h2 (c)
In addition, Fig. 5 shows the direct numerical integration of the Reference Model (solid blue line) and of the slow-flow (red line) for situations corresponding to Figs. 4c and4d. In Fig. 5a, we can see that only 1 ( ) is a strongly amplitude modulated signal (i.e. SMR) which corresponds effectively to relaxation oscillations of the slow-flow in the 1 -direction. A large amplitude periodic regime (which is a harmful situation) corresponding to a large stable fixed point of the slow-flow is observed for the reference model signals in Fig. 5b. Now, the relevance of the theoretical estimation of the trajectory with respect to the parameter is investigated. The comparison is first performed in term of the mitigation limit. For this purpose, the following quantities are defined:

• th ml : the theoretical value of the mitigation limit, • num ml,SF : the numerical estimation of the mitigation limit "measured" on the graphs of the maximum steady-state amplitude obtained from numerical simulations of the slow-flow,

• num ml,RM : the numerical estimation of the mitigation limit "measured" on the graphs of the maximum steady-state amplitude obtained from numerical simulations of the Reference Model.

Then the comparison is completed checking the differences in term of NES efficiency rate, denoted Λ (in %) and defined as the relative difference between the mitigation limit and the Hopf bifurcation point of the airfoil model without NES (71c)

This quantity gives indeed a direct estimation of the efficiency of the NESs: the larger it is, the more efficient the NESs network is. The results of the comparison are reported Tab. 1. One can observe that the mitigation limits obtained from numerical simulations of the slow-flow and of the reference model are very close. Moreover, the comparison between theoretical and numerical values shows a good agreement. Indeed, a difference of the order of magnitude of is observed. However, in term of NES efficiency rate, we can see that in efficient situations (i.e. = 2 ⋅ 10 -2 ) the theoretical prediction can lead to a relatively large overestimation (a relative error of 16.28 %). As previously, these overestimations are due to the 0-order approximation of the asymptotic approach presented in Sect. 3. This could be reduced by considering larger orders of approximation.

we want to mitigate using the three NESs. We focus again on the possible mitigation regimes which can occur when the primary system is unstable.

First, in the Fig. 7 the CM , computed from (47) and ( 52), is depicted in the ( 1 , 2 , 3 )-space with their stable parts (in gray), i.e. the (with = 1, … , 2 = 8), and unstable parts (in black). The stable and unstable fixed points, computed from (63) and (69), are again represented by magenta ans blue points respectively. We can see in each of the four figures a stable fixed point with a large value on the unbounded straight line (i.e. 2 = 8 ). As in Fig. 4, this fixed point corresponds to a large amplitude LCO and if it is reached, we are in a harmful situations.

The theoretical trajectory of the slow-flow is determined for the four values of the bifurcation parameters Θ (in green) and compared to the trajectory obtained from the numerical simulation of the slow-flow (in red). The Fig. 7a shows an harmless situation through a periodic regime, i.e. a stable fixed point of the slow-flow (in magenta) which is not in 8 (here it is in 1 ) is reached. In Figs. 7b and 7c harmless situations through SMRs are shown. Indeed, relaxation oscillations of the slow-flow are observed with respect to 1 in Fig. 7b and to 1 and 2 in Fig. 7c. Finally, because a stable fixed point is reached in 8 , the Fig. 7d corresponds to an harmful situation.

Similar observations as in Fig. 4 can be expressed: (1) during the slow epochs on the theoretical and numerical trajectories of the slow-flow are again almost superimposed and (2) due to the 0-order approximation of the asymptotic approach, differences are observed during fast epochs. Indeed, in the latter case, unlike the theoretical trajectory, the numerical trajectory is not at each jump parallel to only one direction of the ( 1 , 2 , 3 )-space. Again, figures similar to Fig. 7 using = 10 -4 would show a perfect agreement between theoretical and numerical trajectories.

The Fig. 8 shows the direct numerical integration of the Reference Model (blue lines) and of the slowflow (red lines) for situations corresponding to Figs. 7c and7d. In Fig. 8a, 1 ( ) and 2 ( ) are strongly amplitude modulated signals (i.e. SMRs) which corresponds to relaxation oscillations of the slow-flow in both 1 and 2 directions. A large (harmful) amplitude periodic regime corresponding to a large stable fixed point of the slow-flow is observed in Fig. 8b. Now the theoretical value of the mitigation limit Θ th ml is compared to the numerical estimation of the mitigation limit "measured" on the graphs of the maximum steady-state amplitude obtained from numerical simulations of the slow-flow Θ num ml,SF and from numerical simulations of the Reference model Θ num ml,RM . Again, the comparison is also performed in term of NES efficiency rate Λ (in %), now defined with respect to the parameter Θ Λ th ml = 100 × Θ th ml -Θ wo

) 73c
The comparison is presented in the Tab. 2 for both configurations 1 and 2. Moreover, numerical bifurcation diagrams for the variables 1 and 2 (calculated from [START_REF] Lee | Suppression of limit cycle oscillations in the van der Pol oscillator by means of passive non-linear energy sinks[END_REF]) of the reference model are shown in Fig. 9 (in red). The diagrams are compared to diagrams obtained for the primary system alone (in blue). Each of them is obtained from direct numerical integration of the corresponding system of equations and plotting the maximum steady-state amplitude, denoted ( = 1, 2), of each variable as functions of Θ. Once more, similar conclusions as in Sect. 4.1 can be made from Tab. 2. Indeed, Tab. 2a shows that the mitigation limits obtained from numerical simulations of the slow-flow and of the reference model are very close and a difference of the order of magnitude of is observed between theoretical and numerical values. In term of NES efficiency rate (see Tab. theoretical prediction can lead to a relatively large overestimation of the NES efficiency rate, for the same reasons as above.

As a final step, the influence of the number of NESs is examined. The results are shown in Tab. 3 in which the configuration 1 is considered for a number of NESs from = 1 to = 4. In the case of = 4 a fourth NES is added with the following parameters: ℎ 4 = 0.1, ℎ 4 = 6 and 4 = 0.6. We can see that the prediction is relatively robust with respect to number of NESs even if a global slight increase of the relative error between theoretical and numerical efficiency rates is observed when increases.

Conclusion

In this paper, dynamic instability mitigation of a multi-DOF system using Nonlinear Energy Sinks (NESs) is considered. The steady-state regimes of the coupled system is classified into two categories depending on whether the dynamic instability is mitigated or not and therefore separating harmless situations from harmful situations. The value of a chosen bifurcation parameter which separates harmless situations from harmful situations is called mitigation limit.

The prediction method of steady-state regimes has three main steps. The first step is the diagonalization of the primary system using the so-called biorthogonal transformation. Assuming a primary system with only one unstable mode the diagonalized system is reduced ignoring the stable modes and keeping only the unstable mode. The second step consists in applying the complexification-averaging method for the purposes of obtaining the slow-flow of the reduced system. It appears that the slow-flow is a slow-fast system. In light of that fact, in the third step, performing an asymptotic analysis using geometric singular perturbation theory of slow-fast systems, the critical manifold and the fixed points of the slow-flow are computed. The analysis shows that the critical manifold of the system can be reduced to a one dimensional parametric curve evolving in a multidimensional space. The response regimes of the system (and therefore the mitigation limit of the NES setup) is predicted by locating attracting parts, repelling parts and fold points of the critical manifold Table 3: Comparison between numerical and theoretical estimations of the mitigation limit and the NES efficiency rate for several number of NESs attached on the airfoil for = 10 -3 . Configuration 1 is considered for a number of NESs from = 1 to = 3. In the case of = 4, a fourth NES is added with the following parameters: ℎ 4 = 0.1, ℎ 1 = 6 and 3 = 0.7. Finally, the method is applied to the prediction of the mitigation limit of two self-excited mechanical systems coupled to a NES setup. The first one is the well-known Hultén's model which reproduces the typical dynamic behavior of friction systems (it can undergo mode coupling instabilities) and coupled to two ungrounded NESs. The second one is an airfoil model undergoing an aeroelastic instability and coupled to a multi-NES setup (up to four). Theoretical results are compared, for validation purposes, to direct numerical integration of the system. The comparison shows a good agreement and highlights the limits of the asymptotic approach.

Number of

As a perspective, a generalization of this analysis could be suitable in the case of a primary system undergoing a multi-instability, as can be seen in some finite element models of braking systems. Indeed, during the first step of the method, not only one but several unstable modes would be kept. However, as a result of the increase of the critical manifold dimension, the prediction of the steady-state regimes would require a more advanced mathematical procedure. The expression of the nonlinear vector function is given by

(̃ ) = ̃ 3 ̃ 3 , (B.5)
and again ̃ -̃ , with given by ( 5). Following Eqs. ( 6) and ( 7), we define

Figure 1 :

 1 Figure 1: (a) Example of typical critical manifold of a system containing two parallel NESs ( = 2) in the ( 1 , 2 , )-space (i.e. ℝ + 3 ). (in red) is the intersection between the blue surface = √ 1 ( 1 ) and the green surface = √ 2 ( 2 ). The parameters of the primary system are those used in Sect. 4.2 resulting in = 0.82. The NESs parameters are: ℎ 1 = 0.3, ℎ 2 = 0.2, ℎ 1 = ℎ 2 = 6, 1 = 0.028 -0.923 and 2 = 0.030 -0.819 (corresponding to the second and the third NESs in the example of Sect. 4). (b) Projection of the CM in the ( 1 , 2 )-plane.

Figure 2 :

 2 Figure 2: Part of the Critical Manifold (CM) corresponding to the -th NES, Eq (46). The parameters of the primary system are those used in Sect. 4.2 resulting in = 0.82. The NESs parameters are: ℎ = 0.2, ℎ = 6, = 0.030 -0.819 (corresponding to the third NES in the example of Sect. 4).

Figure 3 :Figure 4 :

 34 Figure 3: (a) Hultèn's model with two NES. (b) Zoom on the NES 1 . (c) Zoom on the NES 2 . CM (stable) CM (unstable) • Unstable fixed points • Stable fixed points Theoretical trajectory Trajectory from numerical integration of (31)

Figure 5 :

 5 Figure 5: Direct numerical integration of the Reference Model (solid blue line) and of the slow-flow (red dashed line). Parameters: the same as those used for (a) Fig. 4c and (b) Fig. 4d.

Figure 7 :

 7 Figure 7: Critical manifold in the ( 1 , 2 , 3 )-space defined through Eq. (40) together with theoretical and numerical trajectories of the slow-flow (31). Configuration 1 is considered with = 10 -3 .

Figure 8 :

 8 Figure 8: Direct numerical integration of the Reference Model (solid blue line) and of the slow-flow (red dashed line). Parameters: the same as those used for (a) Fig. 7c and (b) Fig. 7d.

Figure 9 :

 9 Figure 9: Comparison between numerical bifurcation diagrams of the airfoil model with and without NESs with respect to the reduced speed of the flow Θ. Configuration 1 is considered with (a) = 10 -4 , (b) = 10 -3 and (c) = 10 -2 .

2 ̃

 2 For convenience, Eq. (B.1) is written in dimensionless form̃ ′′ + ̃ ′′ + Ω 2 ̃ + ̃ 3 + ̃ ′ + ΛΘ Θ ̃ + ̃ ′ + ∑ ′′ + ̃ ′′ + 2 ̃ + ̃ 3 + ̃ ′ -ΛΘ Θ + ̃ ′ + ∑ + ̃ ℎ ̃ ′ + h′ -̃ ′ + ℎ ̃ + h -̃ 3 = 0 ( = 1, … , ) (B.2c)where the time has been replaced by the dimensionless time = (with= √ ∕ ) and time derivative ∕ is denoted ε ′ ε. The dimensionless displacements are defined by = ∕ , ℎ = ∕ . Moreover, = ∕ ( = 1, … , ), = ∕( ), Ω = ∕ (with = √ ∕ ) and = ∕ . = √ ∕( 2 )is the radius of gyration of the cross section of the wing. The dimensionless nonlinear stiffnesses and damping coefficients are = 2 2 ), ̃ ℎ = ̃ ∕( ). The mass and density ratios are defined by ̃ ℎ = ̃ ∕( 2 ) and = ∞ ∕(2 ) respectively. Finally, the bifurcation parameter under consideration is the reduced speed of the flow Θ = ∕( ). The system of equations (B.2) has the general form of Eq. (3

ℎ

  ̃ ℎ = ℎ , ̃ ℎ = ℎ ( = 1, … , ) ′′ + ℎ ′ + ℎ ′ -′ + ℎ + ℎ -3 = 0 ( = 1, … , ) (B.8c) with ℎ = ℎ ∕ ℎ and ℎ = ℎ ∕ ℎ ( = 1, … , ).

Table 1 :

 1 Comparison between numerical and theoretical estimations of the mitigation limit and the NES efficiency rate for several values of the little parameter .

	2 ⋅ 10 -4	2 ⋅ 10 -3	2 ⋅ 10 -2

Table 2 :

 2 Comparison between numerical and theoretical estimations of (a) the mitigation limit and (b) the NES efficiency rate for several values of the little parameter .

		(a)
		Mitigation limit	10 -4	10 -3	10 -2
	Configuration 1	Θ th ml Θ num ml,SF Θ num ml,RM	0.93448 0.93446 0.93444	0.9471 0.9460 0.9459	1.057 1.021 1.020
	Configuration 2	Θ th ml Θ num ml,SF Θ num ml,RM	0.93430 0.93428 0.93426	0.9453 0.9444 0.9442	1.038 1.013 1.013
		(b)
		NES efficiency rate	10 -4	10 -3	10 -2
	Configuration 1	Λ th ml Λ num ml,SF Λ num ml,RM	0.16 0.16 0.16	1.51 1.39 1.38	13.29 9.43 9.32
	Configuration 2	Λ th ml Λ num ml,SF Λ num ml,RM	0.14 0.14 0.14	1.32 1.22 1.2	11.25 8.57 8.57

  together with position and stability properties of the fixed points of the slow-flow. It is worth noting that the number of NESs does not appear as a theoretical restriction of the prediction method.

					NESs attached ( )
				1	2	3	4
	Mitigation limit	Θ th ml Θ num ml,SF Θ num ml,RM	0.9392 0.9389 0.9388	0.9440 0.9432 0.9430	0.9471 0.9460 0.9459	0.9486 0.9469 0.9467
	NES efficiency rate	Λ th ml Λ num ml,SF Λ num ml,RM	0.66 0.63 0.62	1.17 1.08 1.06	1.51 1.39 1.38	1.64 1.47 1.45
	Relative error	100 ×	Λ th ml -Λ num ml,RM ml,RM Λ num	6.45 %	10.38 %	9.42 %	13.1 %

= + ,[START_REF] Boroson | Optimization under uncertainty of parallel lnonlinear energy sinks[END_REF] with

= -1 assuming a weak instability i.e. that Re 1 is in the order of ( ) (Re 1 = with is in the order of (1)).

This step is not necessary for the variable 1 because it is already a complex variable.

is hyperbolic if all eigenvalues of the Jacobian matrix of the differential systems (38b)(38c) have nonzero real part.

where

The function 1 will be explicitly given here after. Note that 1 is a function of a real variable (here 1 ) because all variables (with ∈ [2, ]) are linked to 1 ( ( ) = 1 ( 1 ), see Eq. (45a)).

From Eq. (61), it is possible to detect fixed points as

and folded singularities as

Folded singularities correspond to situations for which fixed points and fold points coincide, they are hints of possible canard explosions [START_REF] Benoît | Chasse au canard ("duck hunting[END_REF]. The remaining of the section is dedicated to the analysis of the fixed points which, for 0 < ≪ 1, are assumed to be the fixed points of [START_REF] Gendelman | Targeted energy transfer in systems with external and self-excitation[END_REF].

Combining Eq. ( 62) with Eqs. (45), Eqs. (43) and Eq. ( 46), and using the notations introduced Eq. (32b), the function

with 0 = 0 + 0 , 1 = 1 + 1 and 1 = 1 + 1 where 1 is the -component of the vector

Hence, finding the fixed points of Eq. ( 35) goes back to find the roots of

which can be easily solved with a regular computer (at least until ≈ 5 or 6).

To check the stability of a fixed point (denoted * = * 1 , … , * ), Eq. ( 61) is written as

where

Flutter mitigation in an airfoil system using NESs

The system is shown in Fig. 6. It is composed by a 2-DOF nonlinear airfoil model coupled to NESs. Following Sect. 2.1, the equations of motion of the system take the form (see Appendix B)

where, similarly as in Sect. 4.1, Eqs. (72) correspond to Eqs. ( 9) and the corresponding Eqs. [START_REF] Lee | Suppression aeroelastic instability using broadband passive targeted energy transfers, part 1: Theory[END_REF] and Eqs. [START_REF] Grinberg | Response regimes in linear oscillator with 2dof nonlinear energy sink under periodic forcing[END_REF] hereafter referred indistinctly as the Reference Model can be easily deduced as well as the slow-flow model corresponding to Eqs. [START_REF] Nayfeh | Perturbation Methods[END_REF].

The bifurcation parameter under consideration is now the reduced speed of the flow Θ.

In the first instance, a model with three NESs ( = 3) is considered and we use the following parameters of the primary system:

= 1, = 1, = 0.01, = 0.01, = 0.2, = 0.5, Ω = 0.5, = 0, Λ = 2 and = 1∕ [START_REF] Lee | Suppression of limit cycle oscillations in the van der Pol oscillator by means of passive non-linear energy sinks[END_REF]. Two sets of NESs parameters are used, hereafter referred as configuration 1 and configuration 2. The configuration 1 is defined with ℎ 1 = 0.4, ℎ 2 = 0.3, ℎ 3 = 0.2, ℎ 1 = 7, ℎ 2 = 6, ℎ 3 = 6, 1 = -0.9, 2 = -0.8 and 3 = -0.7 which considers almost identical NESs. The configuration 2 is defined with ℎ 1 = 0.35, ℎ 2 = 0.3, ℎ 3 = 0.25, ℎ 1 = 7, ℎ 2 = 6, ℎ 3 = 1, 1 = -0.9, 2 = -0.85 and 3 = -0.8 where the third NES has a small cubic stiffness compared to those of other NESs. Moreover, for both configurations ℎ 1 = ℎ 2 = ℎ 3 = 1 and several values for will be considered.

With this set of numerical values, the airfoil model (i.e. the primary system) undergoes an Hopf bifurcation at Θ wo Hopf = 0.933 due to a mode coalescence phenomenon. If Θ < Θ wo Hopf the system is stable and if Θ > Θ wo Hopf the system is unstable. As previously, in unstable situations, the system reaches an harmful LCO Appendix A. The nonlinear Hultèn's model with two NESs

The system is similar to the one used in [START_REF] Bergeot | Mode coupling instability mitigation in friction systems by means of nonlinear energy sinks : numerical highlighting and local stability analysis[END_REF], it is composed of the nonlinear Hultèn's model [START_REF] Sinou | Mode coupling instability in friction-induced vibrations and its dependency on system parameters including damping[END_REF] coupled to two ungrounded NESs with masses ̄ ℎ , linear stiffnesses ℎ , damping coefficients ̃ ℎ and cubic stiffnesses ℎ ( = 1, 2, see Fig. 3). The motion equations of the system are

where the time derivative is denoted ε . ε. In this model the tangential force due to friction contact is assumed to be proportional to the normal force as given by Coulomb's law: = , where is the friction coefficient.

Introducing the following notation:

) (with = 1, 2) for the primary system parameters and and ̃ ℎ = ̄ ℎ ∕( 21 ), ̃ ℎ = ̃ ℎ ∕( 1 ) and = ∕( 2 1 ) (with = 1, 2) for the NESs parameters, Eqs. (A.1) become

The system of equations (A.2) has the general form of Eq. ( 3) with

The nonlinear vector function contains cubic nonlinearities

and the nonlinearities of the NESs have the required form, i.e.

̃ -̃ , with given by ( 5). According to the method proposed in Sect. 2.1, the small parameter is introduced to take into account small inertia and damping terms of the NESs

and the variables ̃ ( = 1, 2) and h ( = 1, 2) are rescaled through as

Consequently, the equations of motion have the following final form

Appendix B. The airfoil--NESs model

The two degrees of freedom airfoil model coupled to several NESs (see Fig. 6) and used in Sect. 4 is described in this appendix.

Assuming small angle, we derive the equations of motion of the airfoil--NESs model as

where again the time derivative is denoted ε . ε. = ∕2 is the semichord length. is the aerodynamic center, the elastic axis, the center of gravity of the airfoil. is the location aerodynamic center measured from (positive ahead of ). and are the mass of the airfoil and its mass moment of inertia with respect to .

= is the mass unbalance in the airfoil. and are the linear heave and pitch stiffnesses respectively whereas and are the cubic heave and pitch stiffnesses. and are the heave and pitch damping coefficients.

is the constant and uniform flow speed around the airfoil and = 1 2 ∞ 2 is the dynamic pressure where ∞ is the density of the flow.

is the planform of the airfoil, Λ is the lift curve slope and ( = 1, … , ) are the offset attachments of the NESs to the airfoil. Finally, ̃ , ̃ and ( = 1, … , ) are the masses, the damping coefficients and the cubic stiffnesses of the NESs respectively.