B Bergeot

S Berger

S Bellizzi

Mode coupling instability mitigation in friction systems by means of nonlinear energy sinks: Numerical highlighting and local stability analysis

Keywords: Friction-induced vibration, Passive control, Non linear energy sink, Relaxation oscillations, Strongly modulated response

In this paper, we study a problem of passive control of friction-induced vibrations due to mode coupling instability in breaking systems. To achieve that, the well-known two degrees of freedom Hultèn's model, which reproduces the typical dynamic behavior of friction systems, is coupled to two ungrounded Nonlinear Energy Sinks (NES). The NES involves an essential cubic restoring force and a linear damping force. First, using numerical simulations it is shown that the suppression or the mitigation of the instability is possible and four steady-state responses are highlighted: complete suppression, mitigation through periodic response, mitigation through strongly modulated response and no suppression of the mode coupling instability. Then the system is analyzed applying complexification-averaging method, the resulting slow-flow is finally analyzed using geometric singular perturbation theory. This analysis allows to explain the observed steady state response regimes and predict some of them. The boundary values of the friction coefficient for some of the transitions between these regimes are predicted. However, the appearance of a three-dimensional super-slow flow subsystem highlights the limitation of the local linear stability analysis of the slow-flow to predict all these boundaries.

Introduction

Self-excited systems play a key role in numerous industrial applications related to the fields of aeronautics, railways, and cars. Dry friction systems are good examples of these systems (Sinou et al. [2006a,b], [START_REF] Chevennement-Roux | Flexible Wiper System Dynamic Instabilities: Modelling and Experimental Validation[END_REF], Sinou and Jézéquel [2007], [START_REF] Hervé | Analysis of squeal noise and mode coupling instabilities including damping and gyroscopic effects[END_REF]). They develop dynamic instabilities related to the friction which are explained in major cases by two main families of mechanisms. The first family explains the instabilities by the variation of the friction coefficient with respect to the relative speed or by a higher static friction coefficient than the dynamic one. The stick-slip is a well-known phenomenon in this context (Van De Velde and De Baets [1998b,a]). The second family attributes the appearance of instabilities to the sprag-slip mechanism and more generally to the mode-coupling phenomenon. In this case, self-excited oscillations may occur even with a constant friction coefficient. In most brake models the instability is due to mode-coupling phenomenon [START_REF] Fritz | Investigation of the relationship between damping and mode-coupling patterns in case of brake squeal[END_REF], [START_REF] Oden | Models and computational methods for dynamic friction phenomena[END_REF]) which is studied in this paper. Moreover, it has been shown that the well-known two degrees of freedom Hultèn's model [START_REF] Hultén | Friction phenomena related to drum brake squeal instabilities[END_REF][START_REF] Hultén | Brake squeal -a self-exciting mechanism with constant friction[END_REF]) is sufficient to investigate the mode-coupling instability (D 'Souza and Dweib [1990], [START_REF] Eriksson | Friction behaviour and squeal generation of disc brakes at low speeds[END_REF], [START_REF] Hoffmann | Effects of damping on mode-coupling instability in friction induced oscillations[END_REF], [START_REF] Nechak | Non-intrusive generalized polynomial chaos for the robust stabil-ity analysis of uncertain nonlinear dynamic friction systems[END_REF]). In a nutshell, friction systems and especially the braking systems are subject to dynamic instabilities leading to limit cycle oscillations that may affect their efficiency and the user comfort. Moreover, it is very difficult to design completely stable systems, particularly because of the dispersion of friction laws. It is therefore necessary to attenuate these vibrations.

The concept Targeted Energy Transfer (TET) is a relatively new passive control method which consists in coupling an essentially nonlinear attachment also named Nonlinear Energy Sink (NES) to an existing primary system prone to unwanted vibrations. TET has been extensively studied numerically, theoretically and more rarely experimentally. The results prove that the NES is very efficient for vibration mitigation [START_REF] Vakatis | Nonlinear Targeted Energy Transfer in Mechanical and Structural Systems[END_REF]) and noise reduction [START_REF] Bellet | Experimental study of targeted energy transfer from an acoustic system to a nonlinear membrane absorber[END_REF]). Impulsive loading was theoretically analyzed for example by [START_REF] Vakakis | Energy pumping in nonlinear mechanical oscillators: Part II -Resonance capture[END_REF] where TET is investigated in terms of resonance capture. [START_REF] Starosvetsky | Strongly modulated response in forced 2dof oscillatory system with essential mass and potential asymmetry[END_REF] investigate harmonic forcing where response regimes are characterized in terms of periodic and strongly modulated responses using an asymptotic analysis (multi scale approach) of the averaged flow obtained using the complexification-averaging method (Manevitch [1999]). Ahmadabadi and [START_REF] Nili | Annihilation of high-amplitude periodic responses of a forced two degrees-of-freedom oscillatory system using nonlinear energy sink[END_REF] investigate the role of a single degree of freedom NES with nonlinear damping characteristics in annihilating undesired periodic response regimes and simultaneously preserving strongly modulated responses (SMR). The studied system consists in a harmonically excited three degrees of freedom system consisting of two linear coupled oscillators and an NES attached to it. [START_REF] Gourc | Chatter Control in Turning Process with a Nonlinear Energy Sink[END_REF] use a NES to reduce chatter vibration in turning process. An application of NES as a nonlinear vibration absorber in rotor dynamics can be found in [START_REF] Bab | Lateral vibration attenuation of a rotor under mass eccentricity force using non-linear energy sink[END_REF] where the efficiency of a collection of NES is analyzed for vibration mitigation of a rotating system under mass eccentricity force. [START_REF] Bab | Vibration mitigation of a rotating beam under external periodic force using a nonlinear energy sink (nes)[END_REF] investigated the performance of a NES to mitigate vibration of a rotating beam under an external forced. We can also cite [START_REF] Farid | Tuned pendulum as nonlinear energy sink for broad energy range[END_REF] which study the applicability of common pendulum as the NES for mitigation of impulsive excitations. The authors present a theoretical analysis of the damped targeted energy transfer into the pendulum NES from the primary mass with an account of corrections caused by the effect of gravity.

NES are also used to control dynamic instabilities. The possible suppression of the limit cycle oscillations of a Van der Pol oscillator utilizing a NES is demonstrated numerically in [START_REF] Lee | Suppression of limit cycle oscillations in the van der Pol oscillator by means of passive nonlinear energy sinks[END_REF]. In Gendelman and [START_REF] Gendelman | Bifurcations of selfexcitation regimes in a Van der Pol oscillator with a nonlinear energy sink[END_REF] (resp. [START_REF] Domany | Dynamic responses and mitigation of limit cycle oscillations in Van der Pol-Duffing oscillator with nonlinear energy sink[END_REF]), the self-excitation response regimes of a Van der Pol (resp. Van der Pol-Duffing) oscillator with a NES are investigated. An asymptotic analysis of the system related to slow/super-slow decomposition of the averaged flow reveals periodic responses, global bifurcations of different types and basins of attraction of various self-excitation regimes. A series of papers by Lee et al. [2007a,b], Gendelman et al. [2010] demonstrated that a NES coupled to a rigid wing in subsonic flow can partially or even completely suppress aeroelastic instability. In Lee et al. [2007a], the suppression mechanisms are investigated numerically. Several aspects of the suppression mechanisms are validated experimentally in [START_REF] Lee | Suppression aeroelastic instability using broadband passive targeted energy transfers, part 2: Experiments[END_REF]. Moreover, an asymptotic analysis is reported in Gendelman et al. [2010] demonstrating the existence of the three passive suppression mechanisms based on TET. Suppression of aeroelastic instability of a general nonlinear multi degree of freedom system has also be considered in [START_REF] Luongo | Aeroelastic instability analysis of nes-controlled systems via a mixed multiple scale/harmonic balance method[END_REF]. A theoretical/numerical analysis of the capacity of a NES to control helicopter ground resonance instability (which is a mode-coupling instability) has been performed by Bergeot et al. [2016a,b]. More generally, the discussion on relationship between dimensionality of the super-slow manifold, structure of the fixed points and the observed response regimes is explored in review paper by [START_REF] Gendelman | Targeted energy transfer in systems with external and self-excitation[END_REF].

In this context, the use of NES appears to be an interesting way to control mode-coupling instability in braking systems. The goal of the paper is therefore to study the effect of coupling two NES to the two degree of freedom model defined by Hultèn. The originality of this paper focuses on two things: (1) to our knowledge, there are no previous studies on this subject for the dry friction system and (2) in the context of the analysis of this kind of systems, the appearance of a three-dimensional super-slow flow subsystem (the vocabulary will be clarified in the paper) highlights the limitation of the local linear stability analysis of the slow-flow to predict all the steady-state response regimes of the system.

The paper is organized as follows. In Sect. 2, the system under study is presented. It consists in a Hultén's model coupled to two ungrounded NES. In Sect. 3, the local stability analysis of the trivial solution for the with-NES and without-NES systems is presented. Using numerical simulations, the Sect. 4 presents some steady-state response regimes which result from the NES attachments. We count four regimes classified into two categories depending on the fact that the trivial solution of the coupled system is stable or not. The simple local stability analysis per- formed in Sect. 3 is sufficient to find if the NES are able to to suppress completely the dynamic instability. By contrast, trying to describe and predict the other mechanisms (i.e. when the trivial equilibrium remains unstable despite the presence of the NES), a more technical mathematical development is required. It is performed in Sect. 5 following an analytical procedure based on complexification-averaging method together with geometric singular perturbation theory. Positive results and limitations of this analytical work are discussed in Sect. 6. Finally, additional analysis using numerical simulations and a benchmark of the theoretical results obtained in Sect. 6 are performed in Sect. 7.

System under study

The primary system

In this paper a simple self-excited system proposed by [START_REF] Hultén | Friction phenomena related to drum brake squeal instabilities[END_REF][START_REF] Hultén | Brake squeal -a self-exciting mechanism with constant friction[END_REF] is used. Hultèn's model reproduces the typical dynamic behavior of friction systems. Therefore, it is sufficient to investigate how passive control of friction-induced vibration due to mode-coupling phenomenon by means of NES can be performed. The simplicity of this model allows to develop analytical expressions in order to better understand the role of the NES attachment.

This model is composed of a mass m held against a moving band; the contact between the mass and the band is modeled by two plates supported by two different springs (see Fig. 1). For the sake of simplicity, it is usually assumed that the mass and band surfaces are always in contact. This assumption may be due to a preload applied to the system. The contact can be expressed by two cubic stiffnesses, see for example Sinou and Jézéquel [2007]. Damping is integrated as shown in Fig. 1. The friction coefficient at contact is assumed to be constant and the band moves at a constant velocity. Then it is assumed that the direction of friction force does not change because the relative velocity between the band speed and dx 1 /dt or dx 2 /dt is assumed to be positive. All these assumptions are taken into account in order to study a simple non-linear theoretical 2 degrees-of-freedom system with friction such that the effects of damping on mode coupling instability and the associated analytical developments may be easily investigated. The tangential force F T due to friction contact is assumed to be proportional to the normal force F N as given by Coulomb's law: F T = γF N , where γ is the friction coefficient. Assuming the normal force F N is linearly related to the displacement of the mass normal to the contact surface, the resulting equations of motion can be expressed as

m d 2 x 1 dt 2 + c 1 dx 1 dt + k 1 x 1 -γk 2 x 2 + k N L 1 x 3 1 -γk N L 2 x 3 2 = 0 (1a) m d 2 x 2 dt 2 + c 2 dx 2 dt + k 2 x 2 + γk 1 x 1 + γk N L 1 x 3 1 + k N L 2 x 3 2 = 0. (1b)

Mechanical model with Nonlinear Energy Sinks

Two identical NES with masses m h , linear stiffnesses k h , damping coefficients c h and a cubic stiffnesses k N L h , are attached on the system in an ungrounded configuration (see Fig. 2). Because a NES is an essentially nonlinear oscillator the linear stiffness k h is assumed to be very smaller than cubic stiffness k N L h . This assumption is in agreement with experimental data (see for example [START_REF] Bellet | Experimental study of targeted energy transfer from an acoustic system to a nonlinear membrane absorber[END_REF], [START_REF] Gourdon | Nonlinear energy pumping under transient forcing with strongly nonlinear coupling: Theoretical and experimental results[END_REF], Kerschen et al. [2007]).

Taking into account the NES displacements h 1 (t) and h 2 (t), the equations of motion (1) become

m d 2 x 1 dt 2 + c 1 dx 1 dt + k 1 x 1 -γk 2 x 2 + k N L 1 x 3 1 -γk N L 2 x 3 2 + (a) (b) (c)
c h dx 1 dt - dh 1 dt + k h (x 1 -h 1 ) + k N L h (x 1 -h 1 ) 3 = 0 (2a) m h d 2 h 1 dt 2 + c h dh 1 dt - dx 1 dt + k h (h 1 -x 1 ) + k N L h (h 1 -x 1 ) 3 = 0 (2b) m d 2 x 2 dt 2 + c 2 dx 2 dt + k 2 x 2 + γk 1 x 1 + γk N L 1 x 3 1 + k N L 2 x 3 2 + c h dx 2 dt - dh 2 dt + k h (x 2 -h 2 ) + k N L h (x 2 -h 2 ) 3 = 0 (2c) m h d 2 h 2 dt 2 + c h dh 2 dt - dx 2 dt + k h (h 2 -x 2 ) + k N L h (h 2 -x 2 ) 3 = 0. ( 2d
)
Introducing the following notation ηi = c i / √ mk i ,

ω i = k i /m, ϕ i = k N L i /m (with i = 1, 2), = m h /m, ξ h = k h /m, μ = c h / √ mk 1 and ϕ h = k N L h /m, Eqs. (2) become d 2 x 1 dt 2 + η1 ω 1 dx 1 dt + ω 2 1 x 1 -γω 2 2 x 2 + ϕ 1 x 3 1 -γϕ 2 x 3 2 + μω 1 dx 1 dt - dh 1 dt + ξ h (x 1 -h 1 ) + ϕ h (x 1 -h 1 ) 3 = 0 (3a) d 2 h 1 dt 2 + μω 1 dh 1 dt - dx 1 dt + ξ h (h 1 -x 1 ) + ϕ h (h 1 -x 1 ) 3 = 0 (3b) d 2 x 2 dt 2 + η2 ω 2 dx 2 dt + ω 2 2 x 2 + γω 2 1 x 1 + γϕ 1 x 3 1 + ϕ 2 x 3 2 + μω 1 dx 2 dt - dh 2 dt + ξ h (x 2 -h 2 ) + ϕ h (x 2 -h 2 ) 3 = 0 (3c) d 2 h 2 dt 2 + μω 1 dh 2 dt - dx 2 dt + ξ h (h 2 -x 2 ) + ϕ h (h 2 -x 2 ) 3 = 0, (3d)
with 0 < 1, assuming that the mass of the NES is small with respect to the mass of the primary system. In theoretical and experimental works devoted to the systems with NES, the mass ratio is adopted to stay in a range 0.01-0.1 and this convention will be followed in current work. As it will be demonstrated below, relative smallness of is crucial, if one develops the analytic approach to the problem.

In the remaining of this section few notations are introduced in order to obtain a system written into a form which facilitate the theoretical study performed in Sect. 5.

First, changing the time scale from t to t * = ω 1 t and noting " • " the derivative with respect to time t * , Eqs. (3) take the following form

ẍ1 + η1 ẋ1 + x 1 -γ(1 -ã) 2 x 2 + φ1 x 3 1 -γ φ2 x 3 2 + μ( ẋ1 -ḣ1 )+ α1 (x 1 -h 1 ) + α (x 1 -h 1 ) 3 = 0 (4a) ḧ1 + μ( ḣ1 -ẋ1 ) + α1 (h 1 -x 1 ) + α (h 1 -x 1 ) 3 = 0 (4b) ẍ2 + η2 (1 -ã) ẋ2 + (1 -ã) 2 x 2 + γx 1 + γ φ1 x 3 1 + φ2 x 3 2 + μ( ẋ2 -ḣ2 )+ +α 1 (x 2 -h 2 ) + α (x 2 -h 2 ) 3 = 0 (4c) ḧ2 + μ( ḣ2 -ẋ2 ) + α1 (h 2 -x 2 ) + α (h 2 -x 2 ) 3 = 0, (4d)
with

ω 2 /ω 1 = 1 -ã, φ2 = ϕ 2 /ω 2 1 , φ2 = ϕ 2 /ω 2 1 , α1 = ξ h /ω 2 1 and α = ϕ h /ω 2 1 .
We assume that the parameters η1 , η2 , γ, φ1 , φ2 , μ, α and ã are of order (i.e η1 , η2 , γ, φ1 , φ2 , μ, α, ã ∼ O( ), with 0 < 1). Moreover, because the linear stiffness of the NES is supposed to be very smaller than the cubic stiffness, it is stated that α 1 /α ∼ O( 2 ).

In order to perform asymptotic analysis in next sections the parameters of the system are rescaled taking into account previous assumptions

η 1 = η1 ; ϕ 1 = φ1 ; γ = γ ; (5a) µ = μ ; η 2 = η2 ; ϕ 2 = φ2 ; (5b) a = ã ; α = α ; α 1 = α1 3 , (5c)
with η 1 , η 2 , γ, ϕ 1 , ϕ 2 , µ, a, α, α 1 ∼ O(1).

Using rescaled parameters (5), Eqs. (4) become

ẍ1 + η 1 ẋ1 + x 1 -γ(1 -a) 2 x 2 + ϕ 1 x 3 1 -2 γϕ 2 x 3 2 + µ( ẋ1 -ḣ1 )+ 3 α 1 (x 1 -h 1 ) + α (x 1 -h 1 ) 3 = 0 (6a) ḧ1 + µ( ḣ1 -ẋ1 ) + 2 α 1 (h 1 -x 1 ) + α (h 1 -x 1 ) 3 = 0 (6b) ẍ2 + η 2 (1 -a) ẋ2 + (1 -a) 2 x 2 + γx 1 + 2 γϕ 1 x 3 1 + ϕ 2 x 3 2 + µ( ẋ2 -ḣ2 )+ 3 α 1 (x 2 -h 2 ) + α (x 2 -h 2 ) 3 = 0 (6c) ḧ2 + µ( ḣ2 -ẋ2 ) + 2 α 1 (h 2 -x 2 ) + α (h 2 -x 2 ) 3 = 0. (6d)
System of Eqs. ( 6) is the Rescaled Hulten's Model including NES (RHM+NES). The remaining of the paper is devoted to the analysis of its steady-state regimes.

Linear stability of the trivial solution

Using the notation introduced in Sect. 2.2, Eq. ( 1) reduces to the Rescaled Primary System (RPS)

ẍ1 + η 1 ẋ1 + x 1 -γ(1 -a) 2 x 2 + ϕ 1 x 3 1 -2 γϕ 2 x 3 2 = 0 (7a) ẍ2 + η 2 (1 -a) ẋ2 + (1 -a) 2 x 2 + γx 1 + 2 γϕ 1 x 3 1 + ϕ 2 x 3 2 = 0. (7b)
We focus the analysis on the capacity of the NES attachments to suppress or mitigate vibrations when the primary system is unstable. That is why in this section the local stability of the trivial equilibrium of the RPS ( 7) is computed and compared to the local stability of the trivial equilibrium of the RHM+NES (6).

Local stability of the trivial equilibrium is found by looking the sign of the eigenvalues real parts of the Jacobian matrices of the vector functions F 1 and F 2 evaluated at the trivial equilibrium, denoted respectively DF 1 (0) and DF 2 (0). The vector functions F 1 and F 2 characterize respectively the RPS (7) and the RHM+NES (6) when they are formally written in state-space form

Ẋ = F 1 (X), with X = [x 1 x 2 ẋ1 ẋ2 ] t , ( 8
) and Ẋ = F 2 (X), with X = x 1 x 2 h 1 h 2 ẋ1 ẋ2 ḣ1 ḣ2 t .
(9) The evolution of the real and imaginary parts of the resulting eigenvalues of the Jacobian DF 1 (0) with respect to the rescaled friction coefficient γ is plotted in Fig. 3 and compared to that of the Jacobian matrix DF 2 (0) in Fig. 4 for the following values of the parameters

a = 1, η 1 = 0.4, η 2 = 1.2, ( 10a
) = 0.01, ϕ 1 = 5, ϕ 2 = 0, (10b) µ = 0.4, α = 7, α 1 = 7. ( 10c
)
For the RPS, there are four eigenvalues (two pair of complex conjugate) denoted λ wo i (with i ∈ [1, 4]), they are represented in Fig. 3, only two of these eigenvalues (λ 1 and λ 3 ) are plotted in Fig. 4. Due to the presence of the NES, four additional eigenvalues are observed compared to the system without NES. Because of the weak coupling between the NES and the primary system, these additional eigenvalues are close to the eigenvalues of the uncoupled linearized equation of motion of the NES. Each linearized equation of motion has two eigenvalues given by the following expression:

1 2 -µ ± µ 2 -4α 1 2 . The RHM+NES have therefore eight eigenvalues denoted λ w i (with i ∈ [1, 8]).
We can notice that the trivial equilibrium is hyperbolic for both the RPS (7) and the RHM+NES (6) regardless the value of γ, i.e. all the eigenvalues of DF 1 (0) and DF 2 (0) have non zero real parts (except at the Hopf bifurcation points hereafter defined). In this case, The Hartman-Grobman theorem (e.g., [START_REF] Wiggins | Introduction to Applied Nonlinear Dynamical Systems and Chaos[END_REF], Chap. 3) states that in the vicinity of such a hyperbolic equilibrium point, the nonlinear systems ( 8) and ( 9) have the same qualitative stability as does the corresponding linear systems. This guarantees the validity of the present local stability analysis 1 .

For both RPS and RHM+NES, as usual the Hopf bifurcation points is defined as the particular value of γ for which at least the real part of one of the eigenvalues switches from negative to positive while the real part of the other eigenvalues remains negative. At the Hopf bifurcation points the trivial equilibrium switches from stable to unstable. We notice in Fig. 4 that the presence of the NES shifts the bifurcation to a larger value of γ. This is the linear effect of the NES due the additional damping µ.

For the RPS the bifurcation point is denoted γ b and it is used as the origin for following theoretical study. Therefore, in the remaining of the paper the bifurcation parameter under consideration is the detuning term σ defined as follows

γ = γ b + σ (11)
The other parameters are fixed (see Eq. ( 10)).

For the set of parameters (10), γ b = 1.12 and the bifurcation point of the RHM+NES is equal to γ = 1.43, corresponding to

σ = 0.31. ( 12
)
The parameters (10) are chosen to illustrate with only one set of parameters the potential of the NES to modify the response regimes after the bifurcation 1 In the present paper all fixed points of all system of differential equations are checked to be hyperbolic. point. These modified steady-state regimes are presented in the following section.

Possible steady-state response regimes

The aim of this section is first to present the main steady-state response regimes which may result from the NES attachments and their relevance. For that, the time series x 1 (t * ) and x 2 (t * ), resulting from the numerical integration of the RHM+NES, Eqs. ( 6) and of the RPS, Eqs. (7), are compared in Fig. 5. In both cases, same initial conditions are used, chosen as small perturbation of the trivial solution: ḣ2 (0) = 0.

x 1 (0) = 0.05, x 2 (0) = h 1 (0) = h 2 (0) = ẋ1 (0) = ẋ2 (0) = ḣ1 (0) =
Observing the displacements x 1 (t * ) and x 2 (t * ) (solid red line in Fig. 5) of the RHM+NES, four main types of response regimes which may be generated when a NES is attached on the system are highlighted selecting different values of the parameter σ. They are classified into two categories depending on the fact that the trivial solution of the RHM+NES is stable or not:

• The trivial solution of the RHM+NES is stable:

-Complete suppression (see. Fig. 5(a)). In this case, the additional damping due to the NES attachment stabilizes the system and the mode-coupling instability is completely suppressed.

• The trivial solution of the RHM+NES is unstable:

-Mitigation through Periodic Response (PR) (see. Fig. 5(b)). In this case, the steadystate response regime is periodic with frequency close to 12 .

-Mitigation through Strongly Modulated Response (SMR) (see. Fig. 5(c)). In this case, the steady-state response regime is a quasiperiodic regime which exhibits a "fast" component with frequency close to 1 and a "slow" component corresponding to the envelope of the signal. The term "Strongly modulated response" has been introduced by [START_REF] Starosvetsky | Strongly modulated response in forced 2dof oscillatory system with essential mass and potential asymmetry[END_REF] for the study of a harmonically forced linear system coupled to a NES.

-No mitigation (see. Fig. 5(d)). The NES is not able to mitigate the instability. Indeed, a periodic regime is observed with an am-plitude close to the amplitude of the system without NES.

These four responses are also observed by Lee et al. [2007a] and study theoretically by Gendelman et al. [2010] in the context of the mitigation of aeroelastic instabilities of a rigid wing in subsonic flow by means of a NES. Furthermore, Bergeot et al. [2016a] observed these responses studying control of helicopter ground resonance instability attaching a NES on the fuselage of the helicopter.

Using numerical simulations, the capacity of NES to suppress or mitigate vibrations due dynamic instability in friction systems has been highlighted in this section. Nevertheless, numerical simulations are not sufficient to have a good understanding of the mitigation mechanisms. That is why the RHM+NES is analyzed in the following section.

Asymptotic analysis of the Hulten's Model including NES (RHM+NES)

The analysis presented in this section is first based on Complexification-Averaging method (CA-X) introduced by Manevitch [1999] and discussed in detail by [START_REF] Vakatis | Nonlinear Targeted Energy Transfer in Mechanical and Structural Systems[END_REF]. The CA-X leads to the determination of the slow-flow of the system. This slow-flow is then analyzed using Geometric Singular Perturbation Theory (GSPT) [START_REF] Fenichel | Geometric singular perturbation theory for ordinary differential equations[END_REF], Jones [1995], [START_REF] Desroches | Mixed-mode oscillations with multiple time scales[END_REF]).

The slow-flow

First, to simplify the following calculations, it is convenient to introduce barycentric coordinates u i (t) and v i (t) (with i = 1, 2)

u 1 = x 1 + h 1 , v 1 = x 1 -h 1 , ( 13
)
u 2 = x 2 + h 2 , v 2 = x 2 -h 2 , ( 14
)
and reciprocally,

x 1 = u 1 + v 1 + 1 , h 1 = u 1 -v 1 + 1 , ( 15
)
x 2 = u 2 + v 2 + 1 , h 1 = u 2 -v 2 + 1 . ( 16
)
Using Eqs. (13-16), Eqs. ( 6) are written as follows

ü1 + η 1 ( u1 + v1 ) + 1 + u 1 + v 1 + 1 - (a -1) 2 (γ b + σ)(u 2 + v 2 ) + 1 + ϕ 1 (u 1 + v 1 ) 3 ( + 1) 3 - 2 ϕ 2 (γ b + σ)(u 2 + v 2 ) 3 ( + 1) 3 = 0 (17a) v1 + 2 α 1 (1 + ) v 1 + η 1 ( u1 + v1 ) + 1 + u 1 + v 1 + 1 - (a -1) 2 (γ b + σ)(u 2 + v 2 ) + 1 + µ (1 + ) v1 + α (1 + ) v 3 1 + ϕ 1 (u 1 + v 1 ) 3 ( + 1) 3 - 2 ϕ 2 (γ b + σ)(u 2 + v 2 ) 3 ( + 1) 3 = 0 (17b) ü2 - η 2 (a -1) ( u2 + v2 ) + 1 + (a -1) 2 (u 2 + v 2 ) + 1 + (γ b + σ)(u 1 + v 1 ) + 1 + ϕ 2 (u 2 + v 2 ) 3 ( + 1) 3 + 2 ϕ 1 (γ b + σ)(u 1 + v 1 ) 3 ( + 1) 3 = 0 (17c) v2 + 2 α 1 (1 + ) v 2 - η 2 (a -1) ( u2 + v2 ) + 1 + (a -1) 2 (u 2 + v 2 ) + 1 + (γ b + σ)(u 1 + v 1 ) + 1 + µ (1 + ) v2 + α (1 + ) v 3 2 + ϕ 2 (u 2 + v 2 ) 3 ( + 1) 3 + 2 ϕ 1 (γ b + σ)(u 1 + v 1 ) 3 ( + 1) 3 = 0. (17d)
Secondly, the complexification consists in introducing the following change of variable

ψ 1 = u1 + jω st u 1 , ψ 2 = v1 + jω st v 1 , ( 18a
)
ψ 3 = u2 + jω st u 2 , ψ 4 = v2 + jω st v 2 . ( 18b
)
where j 2 = -1 and ω st is the frequency for which the observed steady-state responses are assumed to oscillate. Here, the observation of the responses made in Sect. 4 lead us to state ω st = 1.

Then, the variable u 1 , v 1 , u 2 , v 2 and their first and second derivatives with respect to time t * are expressed in term of the new variables ψ i (with i ∈ [1, 4]) as:

u 1 = ψ 1 -ψ 1 2j , u 2 = ψ 3 -ψ 3 2j , ( 19a
) u1 = ψ 1 + ψ 1 2 , u2 = ψ 3 + ψ 3 2 , ( 19b
) ü1 = ψ1 - j 2 ψ 1 + ψ 1 , ü2 = ψ3 - j 2 ψ 3 + ψ 3 , ( 19c
)
v 1 = ψ 2 -ψ 2 2j , v 2 = ψ 4 -ψ 4 2j , ( 19d
) v1 = ψ 2 + ψ 2 2 , v2 = ψ 4 + ψ 4 2 , ( 19e
) v1 = ψ2 - j 2 ψ 2 + ψ 2 , v2 = ψ4 - j 2 ψ 4 + ψ 4 , (19f)
where ψ i is the complex conjugate of ψ i . Numerical results shown in Sect. 4 motivate us to assume that the variable u 1 , v 1 , u 2 and v 2 may be broken down into fast and slow components. For that, the following representation is introduced

ψ 1 = φ 1 e jt * , ψ 2 = φ 2 e jt * , ( 20a
)
ψ 3 = φ 3 e jt * , ψ 4 = φ 4 e jt * , ( 20b
)
where φ i (with i ∈ [1, 4]) is the complex slow modulated amplitude of the fast component e jt * . Substituting Eqs. ( 19) into Eqs. ( 17) an equivalent complex system of differential equations is obtained. Then, using Eq. ( 20) in this complex system and performing an averaging over one period equal to 2π yield to a system of equations describing the behavior of the slow complex amplitudes φ i . Finally, since 0 < 1, these equations are expanded in a first-order Taylor series around = 0 giving φ1 = -

1 2 j φ 1 -3ϕ 1 |φ 1 | 2 4 -jη 1 + 1 + φ 3 (γ b + σ) -φ 2 (21a) φ2 = 1 2 j 3αφ 2 |φ 2 | 2 4 + φ 1 -φ 2 (1 -jµ) + 1 2 3jαφ 2 |φ 2 | 2 + 3jφ 1 ϕ 1 |φ 1 | 2 4 + φ 1 (η 1 -j) -φ 2 (µ -1) -jφ 3 (γ b + σ) (21b) φ3 = 1 2 j φ 3 3ϕ 2 |φ 3 | 2 4 + jη -1 -2a + φ 1 (γ b + σ) + φ 4 (21c) φ4 = 1 2 j 3αφ 4 |φ 4 | 2 4 + φ 3 -φ 4 (1 -jµ) + 1 2 j 3αφ 4 |φ 4 | 2 + 3φ 3 ϕ 2 |φ 3 | 2 4 + φ 1 (γ b + σ) -φ 3 (1 + 2a -jη 2 ) + φ 4 (1 + jµ) . (21d)
Eqs. ( 21) describe Complex Form of the Slow-Flow (CFSF) of the system (17) and it does not depend on the linear stiffness part of the NES.

The trivial fixed point is common to both the nonaveraged system (6) (or ( 17)) and the slow-flow ( 21). The stability of this trivial fixed point is calculated in Sect. 3 directly on the RHM+NES. On the other hand, the nontrivial fixed points of the slow-flow (21) (defined as φi = 0 for i ∈ [1, 4]) only characterizes periodic solutions of Eqs. (6) (or ( 17)) if the frequency of the periodic solutions is exactly equal to 1, the frequency used to defined the complex variables (18). However, computing the real form of the slow-flow by using the polar coordinates n i (t) and θ i (t) (with i ∈ [1, 4]), defined by

φ i (t) = n i (t)e jθi(t) , ( 22
)
and considering not the arguments θ i (t) directly but the argument differences δ 2i = θ 2 (t) -θ i (t) (the master component can be chosen arbitrary, φ 2 (t) is chosen for convenience), the periodic solutions of the system of Eqs. (21) (and consequently of the RHM+NES (6)) may be defined as the nontrivial fixed points of the system of the Real Form of the Slow-Flow (RFSF) which may be formally written as follows

Ẋ = F 3 (X), with X = [n 1 n 2 n 3 n 4 δ 21 δ 23 δ 24 ] t . ( 23
)
See Appendix A for more details about the link between the periodic solutions of the RHM+NES (on the form ( 17)) and the fixed points of the RFSF (23).

To summarize, the prediction of a situation in which a fixed points of the RFSF ( 23) is reached allows to predict the existence of stable periodic responses of RHM+NES. However, the numerical solution of F 3 (X) = 0, i.e. the fixed points, cannot be accessed with a regular computer. Therefore, in the following section, an asymptotic analysis of the slow-flow is developed which permits to easily compute the fixed points and to analyze the responses regime when fixed points are unstable.

The Critical Manifold

Slow and Super-slow subsystems. The slowflow described by Eqs. ( 21) can be written as follows:

φ1 = f 1 (φ 1 , φ 2 , φ 3 , φ 4 ) (24a) φ2 = f 2 (φ 1 , φ 2 , φ 3 , φ 4 , ) (24b) φ3 = f 3 (φ 1 , φ 2 , φ 3 , φ 4 ) (24c) φ4 = f 4 (φ 1 , φ 2 , φ 3 , φ 4 , ) , (24d)
which highlights the "slow/fast" nature of the system.

Here terminology introduced by Gendelman and Bar [2010] is preferred, i.e. the terms fast and slow are replaced by slow and super-slow respectively, whereby the term fast denotes the time scale determined by fast oscillations of the primary system with frequency 1. Therefore, system (24) consists of two super-slow complex variables φ 1 and φ 3 and two slow variables φ 2 and φ 4 . Eqs. ( 24) can be reformulated by switching from the slow time scale t * to the super-slow time scale τ = t * as

Φ 1 = f 1 (Φ 1 , Φ 2 , Φ 3 , Φ 4 ) (25a) Φ 2 = f 2 (Φ 1 , Φ 2 , Φ 3 , Φ 4 , ) (25b)
Φ 3 = f 3 (Φ 1 , Φ 2 , Φ 3 , Φ 4 ) (25c) Φ 4 = f 4 (Φ 1 , Φ 2 , Φ 3 , Φ 4 , ) , (25d)
where = d dτ and Φ i (τ ) = φ i t = τ . Solutions of the super-slow/slow system (24) (or (25)) can exhibit slow and super-slow epochs characterized by the speed at which the solution advances.

Stating = 0, the following subsystems are derived from ( 24) and ( 25) respectively:

φ1 = 0 (26a) φ2 = f 2 (φ 1 , φ 2 , φ 3 , φ 4 , 0) (26b) φ3 = 0 (26c) φ4 = f 4 (φ 1 , φ 2 , φ 3 , φ 4 , 0) , ( 26d
)
which is the slow subsystem, and

Φ 1 = f 1 (Φ 1 , Φ 2 , Φ 3 , Φ 4 ) (27a) 0 = f 2 (Φ 1 , Φ 2 , Φ 3 , Φ 4 , 0) (27b) Φ 3 = f 3 (Φ 1 , Φ 2 , Φ 3 , Φ 4 ) (27c) 0 = f 4 (Φ 1 , Φ 2 , Φ 3 , Φ 4 , 0) , (27d)
which is the super-slow subsystem.

In the following sections the Geometric Singular Perturbation Theory (GSPT) is used to describe the dynamics of the full system (24) (and ( 25)) for 0 < 1 from the analysis of the slow and superslow subsystems ( 26) and ( 27) (which are defined for = 0). More precisely, the following result of the GSPT is used: if 0 < 1, the dynamics of the full system (24) (or ( 25)) during slow (resp. superslow) epoch is given by the dynamic of the slow (resp. super-slow) subsystem (26) (resp. ( 27)).

Stability of the Critical Manifold.

The algebraic equations ( 27b) and (27d) of the super-slow subsystem define the so-called Critical Manifold (CM) (Jones [1995])

CM := (z 1 , z 2 , z 3 , z 4 ) ∈ C 4 f 2 (z 1 , z 2 , z 3 , z 4 , 0) = 0, f 4 (z 1 , z 2 , z 3 , z 4 , 0) = 0 . (28)
Looking at Eqs. (21b) and (21d) one see that the functions f 2 (Φ 1 , Φ 2 , Φ 3 , Φ 4 , 0) and f 4 (Φ 1 , Φ 2 , Φ 3 , Φ 4 , 0) have the same form. Therefore, from Eqs. (27b) and (27d) the CM can take the following equivalent form

Φ 1 (τ ) = Φ 2 (τ )F (|Φ 2 (τ )|) (29a) Φ 3 (τ ) = Φ 4 (τ )F (|Φ 4 (τ )|) , ( 29b
)
where the complex function F is defined by

F (X) = F R (X) + jF I (X) = 1 - 3α 4 X 2 -jµ. ( 30
)
Eq. ( 30) shows that the function F , which defines the CM, depends only on the NES parameters α and µ. The fact that Eq. ( 29a) and (29b) are defined with the same function is the direct consequence of using two identical NES.

It is convenient to characterize the CM in R. To achieve this, again polar coordinates are introduced Φ i (τ ) = N i (τ )e jΘj (τ ) for i = 1, 2, 3, 4, (31) and the module and the argument of (29) are successively computed, that lies to

N 2 1 = N 2 2 F R (N 2 ) 2 + F I (N 2 ) 2 = H(N 2 ) (32a) Θ 1 = Θ 2 + arg (F (N 2 )) (32b) N 2 3 = N 2 4 F R (N 4 ) 2 + F I (N 4 ) 2 = H(N 4 ) (32c) Θ 3 = Θ 4 + arg (F (N 4 )) . ( 32d
)
The local extrema of the real function H(X) are given by the positive roots of its derivative H (X). An easy calculus shows that the local extrema occur at

X M = 2 3 2 -1 -3µ 2 α (33) X m = 2 3 2 + 1 -3µ 2 α , ( 34
)
if the following relation holds

µ < 1 √ 3 , ( 35
)
and in this case X M < X m . The chosen value of µ (see Eq. ( 10)) satisfies the previous relation. Because φ2 and φ4 in ( 26) are reduced to zero in the CM, each point of the CM is a fixed point for the slow subsystem which consists in two independent systems of equations

φ1 = 0 (36a) φ2 = j 2 φ 1 -φ 2 F (|φ 2 |) , ( 36b
)
and,

φ3 = 0 (37a) φ4 = j 2 φ 3 -φ 4 F (|φ 4 |) . ( 37b
)
To determine the stability of the CM, one must know if the CM attracts or repels the slow dynamics. For that, using the polar coordinates φ i = n i e jθi , Eqs. ( 36) and ( 37) are rewritten as follows

ṅ1 = 0 (38a) ṅ2 = 1 2 (n 1 sin δ 21 + n 2 F I (n 2 )) (38b) δ21 = 1 2 n 1 n 2 cos δ 21 -F R (n 2 ) , ( 38c
)
and,

ṅ3 = 0 (39a) ṅ4 = 1 2 (n 3 sin δ 43 + n 2 F I (n 4 )) (39b) δ43 = 1 2 n 3 n 4 cos δ 43 -F R (n 4 ) (39c)
where the argument differences δ 21 = θ 2 -θ 1 and δ 43 = θ 4 -θ 3 have been introduced. Stability range of the CM is then determined by examining the sign of the eigenvalues real parts of the Jacobian matrix of the two identical differential systems ( 38) and ( 39) on the CM. It can be shown that the condition of stability of the CM is equivalent to

H (X) > 0, (40)
and the stability range of the CM is characterized by the points (X m , Y m ) and (X M , Y M ) where Y m = H (X m ) and Y M = H (X M ). Such points on which the CM ceases to be hyperbolic3 connecting attractive and repulsive parts of the the CM are called folded singularities. We point that the folded singularities exist whatever the values of the argument differences δ 21 = θ 2 -θ 1 and δ 43 = θ 4 -θ 3 .

A typical Critical Manifold and its stability range are depicted (see Fig. 6) in which X d and X u are solutions of

H (X m ) = H (X d ) ⇒ X d = 2 √ 2 3 1 -1 -3µ 2 α , (41)
and, respectively.

H (X M ) = H (X u ) ⇒ X u = 2 √ 2 3 1 + 1 -3µ 2 α , ( 42
)

Fixed points and folded singularities of the slow-flow

Introducing the CM Eq. ( 29) in Eqs. ( 27a) and ( 27c), the super-slow subsystem is written only with respect to the variables Φ 2 and Φ 4

∂ [Φ 2 F (|Φ 2 |)] ∂τ = f 1 Φ 2 F (|Φ 2 |) , Φ 2 , Φ 4 F (|Φ 4 |) , Φ 4 (43a) ∂ [Φ 4 F (|Φ 4 |)] ∂τ = f 3 Φ 2 F (|Φ 2 |) , Φ 2 , Φ 4 F (|Φ 4 |) , Φ 4 . ( 43b
)
Using the polar coordinates (31) and separating real and imaginary parts, Eqs. (43) take the following forms

F R (N 2 ) ∂N 2 ∂τ -N 2 ∂Θ 2 ∂τ F I (N 2 ) + N 2 ∂N 2 ∂τ F R (N 2 ) = f 1,R (N 2 , N 4 , ∆ 24 ) (44a) N 2 ∂Θ 2 ∂τ F R (N 2 ) + ∂N 2 ∂τ F I (N 2 ) + N 2 ∂N 2 ∂τ F I (N 2 ) = f 1,I (N 2 , N 4 , ∆ 24 ) (44b) F R (N 4 ) ∂N 4 ∂τ -N 4 ∂Θ 4 ∂τ F I (N 4 ) + N 4 ∂N 4 ∂τ F R (N 4 ) = f 3,R (N 2 , N 4 , ∆ 24 ) (44c) N 4 ∂Θ 4 ∂τ F R (N 4 ) + ∂N 4 ∂τ F I (N 4 ) + N 4 ∂N 4 ∂τ F I (N 4 ) = f 3,I (N 2 , N 4 , ∆ 24 ) (44d)
where the expressions of the functions f 1,R , f 1,I , f 3,R and f 3,I are respectively given by

f 1,R (N 2 , N 4 , ∆ 24 ) = Re f 1 N 2 e jΘ2 F (N 2 ) , N 2 e jΘ2 , N 4 e jΘ4 F (N 4 ) , N 4 e jΘ4
× e -jΘ2 , (45)

f 1,I (N 2 , N 4 , ∆ 24 ) = Im f 1 N 2 e jΘ2 F (N 2 ) , N 2 e jΘ2 , N 4 e jΘ4 F (N 4 ) , N 4 e jΘ4
× e -jΘ2 , (46)

f 3,R (N 2 , N 4 , ∆ 24 ) = Re f 3 N 2 e jΘ2 F (N 2 ) , N 2 e jΘ2 , N 4 e jΘ4 F (N 4 ) , N 4 e jΘ4
× e -jΘ4 , (47)

f 3,I (N 2 , N 4 , ∆ 24 ) = Im f 3 N 2 e jΘ2 F (N 2 ) , N 2 e jΘ2 , N 4 e jΘ4 F (N 4 ) , N 4 e jΘ4
× e -jΘ4 , (48) involving the argument difference ∆ 24 = Θ 2 -Θ 4 . Combining Eqs. (44a-44d), system of Eqs. ( 44) can be reduced to the following forms

g(N 2 ) ∂N 2 ∂τ = f N2 (N 2 , N 4 , ∆ 24 ) (49a) g(N 2 ) ∂Θ 2 ∂τ = f Θ2 (N 2 , N 4 , ∆ 24 ) (49b) g(N 4 ) ∂N 4 ∂τ = f N4 (N 2 , N 4 , ∆ 24 ) (49c) g(N 4 ) ∂Θ 4 ∂τ = f Θ4 (N 2 , N 4 , ∆ 24 ) , ( 49d
)
and finally

g(N 2 ) ∂N 2 ∂τ = f N2 (N 2 , N 4 , ∆ 24 ) (50a) g(N 4 ) ∂N 4 ∂τ = f N4 (N 2 , N 4 , ∆ 24 ) (50b) g(N 2 )g(N 4 ) ∂∆ 24 ∂τ = f ∆24 (N 2 , N 4 , ∆ 24 ) , ( 50c
)
where

f ∆24 (N 2 , N 4 , ∆ 24 ) = g(N 4 )f Θ2 (N 2 , N 4 , ∆ 24 ) -g(N 2 )f Θ4 (N 2 , N 4 , ∆ 24 ) , ( 51
)
and

g(X) = H (X) 2 . ( 52
)
From Eqs. ( 50), it is possible to detect fixed points and folded singularities.

f N2 (N e 2 , N e 4 , ∆ e 24 ) = 0 (53a) f N4 (N e 2 , N e 4 , ∆ e 24 ) = 0 (53b) f ∆24 (N e 2 , N e 4 , ∆ e 24 ) = 0 (53c) g(N e 2 ) = 0 (53d) g(N e 4 ) = 0. ( 53e
)
Since 0 < 1, fixed points computed from Eqs. (53) corresponds to fixed points of the real RFSF (23). The expression of these fixed points are obtained from Eq. ( 32)

n e 1 = H(N e 2 ) ( 54a
)

n e 2 = N e 2 (54b
)

n e 3 = H(N e 4 ) (54c
)

n e 4 = N e 4 (54d
)

δ e 21 = -arg (F (N 2 )) (54e
)

δ e 23 = ∆ e 24 -arg (F (N 4 )) (54f) δ e 24 = ∆ e 24 . ( 54g
)
Stability of the fixed points can now be found by looking the sign of the eigenvalues real parts of the Jacobian matrix of the vector function F 3 (see Eq. ( 23)) evaluated at

X e = [n e
1 n e 2 n e 3 n e 4 δ e 21 δ e 23 δ e 24 ] t .

(55)

Folded singularities. The system of Eqs. ( 50) is not defined on the folded singularities defined in Sect. 5.2 because the system is singular on these points. Introducing, the desingularized super-slow subsystem, which is obtained using the time rescaling τ → τ g(N 2 )g(N 4 )

∂N 2 ∂τ = g(N 4 )f N2 (N 2 , N 4 , ∆ 24 ) (56a) ∂N 4 ∂τ = g(N 2 )f N4 (N 2 , N 4 , ∆ 24 ) (56b) ∂∆ 24 ∂τ = f ∆24 (N 2 , N 4 , ∆ 24 ) , (56c)
one can notice that the folded singularities (here denoted {N s 2 , N s 4 , ∆ s 24 }) can be also defined as the fixed points of the desingularized super-slow subsystem (56) (see for example [START_REF] Desroches | Mixed-mode oscillations with multiple time scales[END_REF]). Since g X m/M = 0 (where X m/M denotes indistinctly X M or X m ) and looking at the definition of the function f ∆24 (see Eq. 51), folded singularities are defined from (56) by

{N s 2 , N s 4 , ∆ s 24 } = X m/M , X m/M , ∀∆ 24 . ( 57
)
6 Exploitation of theoretical results and discussion

Explanation of the observed responses

The existence of an S-shaped CM (i.e. the presence of folded singularities on which the stability of the CM changes) allows to explain the three steady-state regimes of the RHM+NES (6) observed when its trivial solution is unstable, namely: Mitigation through Periodic Response, Mitigation through Strongly Modulated Response (SMR) and No mitigation.

A fixed point of the RedM is reached. These situations corresponds to a periodic solutions of the RHM+NES which may correspond to either mitigation through Periodic Response or No mitigation of the instability (see Sect. 4). Indeed, the bifurcation diagram of the super-slow subsystem (50) with respect the variables N 1 (see Fig. The comparison with the bifurcation diagram of the RPS (7)4 (see also Fig. 7) allows to obtain to following conclusion: if a stable fixed point is reached in the domain 1, mitigation through Periodic Response are observed and if a stable fixed point is reached in the domain 2, this time, no suppression of the instability is observed.

Relaxation oscillations. The S-shape of the CM suggests also the possible existence of relaxation oscillations [START_REF] Grasman | Asymptotic Methods for Relaxation Oscillations and Applications[END_REF]): after reaching of fold singularity (X M , Y M ) (in the (N 2 , N 1 )-plane or in the (N 4 , N 3 )-plane), the system jumps to (X u , Y M ) , which is followed by a super-slow evolution of the trajectory of the system (in the stable domain of the CM) until it reaches (X m , Y m ). After another jump and a super-slow evolution (the stable domain of the CM), the trajectory returns to (X M , Y M ) (see Fig. 6). Such scenario of relaxation oscillations for the slow-flow can explain the existence of Strongly Modulated Responses √ 3, the S-shape nature of the CM is lost and therefore relaxation oscillations are not possible.

Prediction of the steady-state response regimes and limitation of the local linear stability analysis of the slow-flow

The aim of this section is to predict, for a given set of parameters, the nature of the resulting steady-state response regimes of the RHM+NES (6). The prediction is performed checking first the local stability property of the trivial equilibrium point (see Sect. 3), this allowed to predict the appearance of Complete suppression of the instability. Then, if the trivial solution is unstable, we intent to discriminate the three other steady-state response regimes analyzing the local stability of the fixed points of the slow-flow (see Sect. 5).

To illustrate the prediction method and highlight its limitation, the study is restricted to a one parameter dependent system, namely the detuning parameter σ. To simplify the presentation of the method (but without limiting its generality), the set of parameters (10) has been chosen to allow a continuous evolution of the nature of the steady-state response regimes when the parameter σ increases (see Fig. 8). In this case, the goal is to find the values σ 1 , σ 2 and σ 3 which separate the regions of existence of the four steady-state response regimes.

From the simple local stability analysis of the trivial equilibrium performed in Sect. 3, it is possible to obtain a theoretical value of σ 1 (denoted σ th 1 ), this is the difference between the Hopf bifurcation point of the system without and with NES (see Eq. ( 12)): σ th 1 = 0.31. ). From this analysis, a theoretical value σ th 2 of σ 2 and a lower σ th 3,l theoretical values of σ 3 can be deduced. Here, lower values mean that it is sure to obtain SMR mechanism if σ < σ th 3,l . However, if σ > σ th 3,l , it is not possible to conclude if SMR or no mitigation mechanisms are encountered. The theoretical value σ th 2 is obtained as the first intersection between the bifurcation diagram of the super-slow subsystem (50) with respect the variables N 2 or N 4 and the line N 2 = N 4 = X M . In Fig. 9(a), one can notice that the bifurcation diagram with respect to the variables N 2 is the first to cross the line N 2 = N 4 = X M at σ = 0.73 which is the theoretical value of σ th 2 . At this point, relaxation oscillations of the slow-flow are observed. Before the appearance of the second branch of stable fixed points at σ ≈ 0.93, relaxation oscillations phenomenon is the only possible regime, this defines the lower theoretical value σ th 3,l .

The fact that only a lower theoretical value of the σ 3 can be obtained highlights the limitation of the method to predict unequivocally all the steadystate response regimes, this is due to the dimensionality of the super-slow flow subsystem. Indeed, in works in which the super-slow flow subsystem can be reduced to a one-dimensional system [START_REF] Domany | Dynamic responses and mitigation of limit cycle oscillations in Van der Pol-Duffing oscillator with nonlinear energy sink[END_REF], Gendelman et al. [2010], Gendelman and [START_REF] Gendelman | Bifurcations of selfexcitation regimes in a Van der Pol oscillator with a nonlinear energy sink[END_REF]), global structure of possible response regimes can be deduced directly from the local stability analysis of the slow-flow. Here the super-slow subsystem is a three-dimensional system. In this case, determination of the basin of attraction of each solution (fixed points and relaxation oscillations) would be needed to conclude about the nature of the steadystate regimes, this may be subject of future work.

Additional analysis using numerical simulations and benchmark of theoretical results

In this section numerical simulations are first performed in order to explore more accurately the possible steady-state response regimes. Secondly, these numerical simulations allow to study the validity of the theoretical results presented in Sect. 6.2. Figs. 10 to 14 present the same type of numerical simulations obtained using the set of parameters (10) and for five different values of σ. In each figure, it is depicted:

• Figs. (a), (b), (c) and (d) compare numerical simulations of the RHM+NES written using barycentric coordinates (see Eq. ( 17)) and numerical simulations of the slow-flow (21).

• Figs. (e) and (f) compare the Critical Manifold (N 1 = H(N 2 ) and (N 3 = H(N 4 ), see Eq. ( 32)) with numerical simulations of the slowflow.

• Figs. (g) and (h) compare the bifurcation diagram obtained from Eq. ( 53) with the trajectory of the slow-flow in the 3D N 2 , N 4 , ∆ 24 -space. Two views are used to make the reading of the 3D graph easier.

In Fig. 10 the chosen value of σ = 0.2 is smaller than σ th 1 = 0.31 and Complete suppression is actually observed. In Fig. 11 the chosen value of σ = 0.6 is between σ th 1 = 0.31 and σ th 2 = 0.73 and Mitigation through Periodic Response is actually observed. In Fig. 12 the chosen value of σ = 0.9 is between σ th 2 = 0.73 and σ th 3,l = 0.93 and Mitigation through Strongly modulated Modulated Response is actually observed. We can notice in this example that relaxation oscillations of the slow-flow which explain the strongly modulated response of the RHM+NES are observed only in the N 1 , N 2 -plane and only small oscillation are observed in the N 3 , N 4 -plane (see Figs. 12(e) and 12(f)). We call this SMR regime: semi-SMR regime. The value of σ must be larger to observe relaxation oscillations in both N 1 , N 2 -plane and N 3 , N 4 -plane (see Fig. 13 obtained using σ = 1.3). We call this SMR regime: full-SMR regime. Finally, in Fig. 14, the results are obtained with σ = 2 and No mitigation is observed.

To observe more precisely the boundaries between the different observed regimes, the amplitudes of the steady-state responses of the numerical simulations of the RHM+NES written using barycentric coordinates (see Eq. ( 17)) and numerical simulations of the slowflow ( 21) are plotted and compared to the theoretical bifurcation diagram in Fig. 15. Following conclusion can be made from the observation of these figures:

Comparison between theoretical and numerical boundaries values of σ. The better figure to perform this comparison is Fig. 15(b). We can observe that the transition from Complete suppression to Mitigation through Periodic Response and the transition from Mitigation through Periodic Response to Mitigation through Strongly Modulated Response is well predicted by σ th 1 = 0.31 and σ th 2 = 0.73 respectively. Indeed, one can see that the steady-state amplitudes of the numerical simulations (for both RHM+NES and the slow-flow) follow the amplitudes predicted by the bifurcation diagram until σ = σ th 2 . At this value a jump of the steady-state amplitudes of the numerical simulations is observed corresponding to the transition from periodic regimes to SMR.

Not predicted boundaries values of

Prediction of the steady-state amplitudes

The theoretical bifurcation diagram is obtained from the super-slow subsystem (see Eq.( 53)). We can first see the good agreement between the theoretical bifurcation diagram and the steady-state amplitudes measures on numerical simulation of the slow-slow when the steady-state regime is a periodic regime, i.e. when σ < σ th 2 and σ > σ num

3

. Moreover the theoretical values H(X u ) = Y M and X u (see Fig. 6) give an approximated prediction of maximum amplitudes when full-SMR occurs, i.e. for 1.08 < σ < σ num 3 .

Comparison between RHM+NES and the slow-flow. The comparison between numerical simulations of RHM+NES and the slow-flow is important because the capacity of the slow-flow to reproduce the behavior of the RHM+NES reflect the quality of all mathematical developments which derive from (even those not presented in this paper). Observing Fig. 15, the comparison shows that the slow-flow well predicts the boundary values of σ. The slow-flow well predicts also the steady-state amplitudes except for periodic regimes in domain 2 (see Fig. 7) corresponding to No mitigation of the dynamic instability. Indeed, one can observe a significant underestimation of the steadystate amplitude (here especially for the variable N 2 , see Fig. 15(b)). However, in the context of the control of the instability the most pertinent values to predict are the boundaries, in particular the boundary σ 3 corresponding to the transition from Mitigation through Strongly Modulated Response to No mitigation of the instability. The fact that the slow-slow can describe this transition motive us to continue this research in order to find a theoretical value of σ 3 .

Conclusion

We studied the capacity of Nonlinear Energy Sinks (NES) to mitigate vibrations due to mode-coupling instability in braking systems. To achieve that, a simple two degree of freedom model which can reproduce mode-coupling instability (i.e. the well-known Hultèn's model) was coupled to two ungrounded NES.

To analyze the steady-state response regimes, the system is partitioned in slow-fast dynamics using complexification-averaging approach. The presence a small dimensionless parameter related to the mass of the NES in the slow-flow system implies that it involves two "slow" complex variables and two "superslow" complex variables. The "super-slow/slow" nature of the system allowed us to use multiple scale approach to analyze it. In particular, the Critical Manifold of the slow-flow was determined. Its S-shape (i.e. involving two folded singularities) and the associated stability properties provide an analytical tool to explain the existence of three regimes: periodic response regimes, strongly modulated responses regimes and no mitigation regimes that appear when the trivial solution is unstable. A complete suppression regimes is also observed and it is studied directly on the full system.

The boundary values of the friction coefficient corresponding to the transition from complete suppression regime to periodic response regimes and from periodic response regimes to strongly modulated responses are predicted analytically. However the prediction of the boundary value between strongly modulated responses and no mitigation responses is not performed, this highlights that global structure of possible response regimes can not be deduced from local stability analysis of a super-slow flow subsystem with dimension larger than one. The prediction of this boundary value could be important in the context of engineering applications. That is why advanced mathematical procedure will be developed to predict this boundary.

Using Eqs. ( 18) and ( 20), the steady-state regime described by ( 58) is written in terms of complex slow modulated amplitude φ i (with i ∈ [1, 4]) φ 1 = ( u1 + ju 1 ) e -jt = U 1 1 + e 2 e j(et+θu 1 ) + t.h.f, (59a) where "t.h.f" means term of higher frequency. Using the polar coordinates n i (t) and θ i (t) (with i ∈ [1, 4]) defined by Eq. ( 22), Eqs. (59) become

φ 2 = ( v1 + jv 1 ) e -jt = V 1 1 + e 2 e j
n 1 = U 1 1 + e 2 ≈ U 1 (60a) n 2 = V 1 1 + e 2 ≈ V 1 (60b) n 3 = U 2 1 + e 2 ≈ U 2 (60c
)

n 4 = V 2 1 + e 2 ≈ V 2 (60d)
δ 21 = θ 2 -θ 1 = θ v1 -θ u1 (60e)

δ 23 = θ 2 -θ 3 = θ v1 -θ u2 (60f) δ 24 = θ 2 -θ 4 = θ v1 -θ v2 . ( 60g
)
Comparing Eqs. ( 59) and (60) we understand why the argument differences must be used. Indeed, even if e 1, the error e caused a linear growth of the arguments of the complex amplitude φ i , which are therefore not stationary. 17)) (blue squares) and from numerical simulations of the slow-flow (21) (red circles). Set of parameters ( 10) is used.

Figure 1 :

 1 Figure 1: Mechanical system without NES.

Figure 2 :

 2 Figure 2: (a) Mechanical system with NES. (b) Zoom on the NES1. (c) Zoom on the NES2.

Figure 3 :

 3 Figure 3: Evolution of real and imaginary parts of the eigenvalues of the Jacobian matrix of the vector function F1 evaluated at the trivial equilibrium. Parameters used: a = 1, η1 = 0.4, η2 = 1.2 and = 0.01. γ b is the Hopf bifurcation point.

Figure 4 :

 4 Figure 4: Evolution of real and imaginary parts of the eigenvalues of the Jacobian matrices of the vector functions F1 (solid black lines) and F2 (dashed gray lines) evaluated at the trivial equilibrium. Parameters used: a = 1, η1 = 0.4, η2 = 1.2, = 0.01, α1 = 7 and µ = 0.03.

Figure 5 :

 5 Figure 5: Comparison between time series x1(t * ) and x2(t * ) resulting from the numerical integration of the RHM+NES (Eqs. (6)) (solid black line) and x1(t * ) and x2(t * ) resulting of the numerical integration of the RPS (Eqs. (7)) (solid gray line). Set of parameters (10) is used with (a) σ = 0.2, (b) σ = 0.6, (c) σ = 1.3 and (d) σ = 2. On the left, the position of each simulation in the graph representing the real parts of the eigenvalues with respect to γ is shown.

Figure 6 :

 6 Figure 6: Critical Manifold (CM). Following parameters are used: α = 7 and µ = 0.4.

Fixed

  points. The (regular) fixed points of Eqs. (50), {N e 2 , N e 4 , ∆ e 24 }, are defined by

  7(a)) and N 3 (see Fig. 7(b)) highlights two ranges of σ in which stable fixed points are encountered (see Fig. 7(a)): • Domain 1: 0.31 σ 0.72, • Domain 2: σ 0.93.

Figure 7 :

 7 Figure 7: Comparison between the bifurcation diagrams of the slow-flow of the system without NES (RPS (7)) and with NES which is obtained from the analysis of the super-slow subsystem (50); with respect to the variable N1 (a) and to the variable N3 (b). Set of parameters (10) is used.

Figure 8 :

 8 Figure 8: Outline schematic showing the definition of the boundary values σ1, σ2 and σ3 which separate the regions of existence of the four steady-state response regimes.

Figure 9 :

 9 Figure 9: Bifurcation diagram of the slow-flow of the system with NES obtained from the analysis of the super-slow subsystem (50); with respect to the variable N2 and N4 (a) and to the variable ∆24 (b). Set of parameters (10) is used.

Fig. 9

 9 shows the bifurcation diagram deduced from local stability analysis of the super-slow subsystem (50) performed in Sect. 5 (Fig. 9(a) with respect to the variables N 2 et N 4 and Fig. 9(b) with respect to the argument difference ∆ 24

  σ. The transition from Mitigation through Strongly Modulated Response to No mitigation and from semi-SMR regime to full-SMR regime have not theoretical values. However, they can be found on numerical simulations. Observing for example Fig. 15(b), one can see that the transition from Mitigation through Strongly Modulated Response to No mitigation of the dynamic instability appears at σ = 1.63 (denoted σ num 3 on the graph). The transition from semi-SMR regime to full-SMR regime is observed on Fig. 15(d), it appears at σ = 1.08.

Figure 10 :

 10 Figure 10: (a), (b), (c) and (d) compare numerical simulations of the RHM+NES written using barycentric coordinates (see Eq. (17)) (solid grey line) and numerical simulations of the slow-flow (21) (solid black line). (e) and (f) compare the Critical Manifold N1 = H(N2) and N3 = H(N4) (dashed grey line) (see Eq. (32)) with numerical simulations of the slow-flow (solid black line). (g) and (h) compare the bifurcation diagram obtained from Eq. (53) (black and grey points for stable and unstable branches respectively) with the trajectory of the slow-flow (solid red line) in the 3D N2, N4, ∆25-space, The blue point shows the initial condition and two views are used to make the reading of the 3D graph easier. Set of parameters (10) is used and σ = 0.2.

Figure 11 :

 11 Figure 11: Same caption as for Fig. 10. Set of parameters (10) is used and σ = 0.6.

Figure 12 :

 12 Figure 12: Same caption as for Fig. 10. Set of parameters (10) is used and σ = 0.9.

Figure 13 :

 13 Figure 13: Same caption as for Fig. 10. Set of parameters (10) is used and σ = 1.3.

Figure 14 :

 14 Figure 14: Same caption as for Fig. 10. Set of parameters (10) is used and σ = 2.

Figure 15 :

 15 Figure 15: Comparison between steady-state amplitudes obtained from numerical simulations of the RHM+NES written using barycentric coordinates (see Eq. (17)) (blue squares) and from numerical simulations of the slow-flow (21) (red circles). Set of parameters (10) is used.

This can be shown for example by computing the power spectrum of the steady part of the signal.

The CM is hyperbolic if all eigenvalues of the Jacobian matrices of the differential systems (38) and (39) have nonzero real part.

The bifurcation diagram of the RPS (7) is also obtained computing the stability of the non trivial fixed points of its slow-flow expanded in fist-order Taylor series around = 0.

A Link between the periodic solutions of the RHM+NES and the fixed points of the RFSF

To explain the link between the periodic solutions of the RHM+NES (on the form of ( 17)) and the fixed points of the RFSF ( 23), let's consider a periodic steady-state regime of the system (17) defined by

where

, θ v2 and e are real constant. e characterizes the error made using 1 as fast component frequency in the complexificationaveraging method with the assumption that e 1.


Dominique Eyheramendy
email: dominique.eyheramendy@ec-marseille.fr

David Loureiro
email: david.loureiro@ens-lyon.fr

Fabienne Oudin-Dardun
email: fabienne.oudin@univ-lyon1.fr

An integrated object-oriented approach for parallel CFD

Keywords: Finite elements, object-oriented programming, domain decomposition, multi-threaded computing, distributed computing

In this paper, we present a global computer science approach to deal with parallel computations. The proposed approach consists in managing at the same level either multithreading or distributed strategies, whatever the computation may be. The integration of the concept is held in a Java framework which proposes both, a pure object-oriented paradigm and, convenient libraries to deal with threads management and communications schemes. The approach is illustrated on a domain decomposition method for a Navier-Stokes flow.

Introduction

In this paper, we present the key ideas of a Java application to a finite element code. This approach is based on the reusability of code and the portability of application in the case of parallel application. The key idea of the developments is that the parallel algorithms are programmed within the application using a single programming concept and a single language. In classical approaches, the programming language (C, C++, Fortran) is associated to additional libraries such as PVM or MPI. Existing libraries in the Java environment are used here. Thus, the maintenance of the application is made easier, which is important from an industrial point of view. The proposed approach offers the advantage of a high level of reusability of the different parts of the code whatever the computational strategy is. The consequence of it is a high reliability of the software for domain decomposition methods. This reliability is obtained through a classical object-oriented paradigm which allows the programmer to separate the management of finite elements data, the solution algorithms, the management of different processes and the communication schemes.

In section 2, we give a state of the art for Java applications in Computational Mechanics that most often includes parallel computations. In section 3, we described two different implementations of an overlapping Schwarz domain decomposition method within the same finite element code. In section 4, we give a numerical application that has been run using the same code on different systems.

Java in Computational Mechanics

Until now, most of the developments in Java computational mechanics have been considered by the computational science community, in general concerned by parallel computations. In the domain of numerical computations, Java retained some attention for its networking capabilities and its Internet easy portability. E.g. in [START_REF] Nuggehally | An internet-based computing plateform for the boundary element method[END_REF], a trivial application based on a boundary element method has been developed. Similar developments may be found directly on the Internet, including basic finite element applications, such as in [START_REF] Nikishkov | The development of a Java engineering application for higher-order asymptotic analysis of crack-tip fields[END_REF] for an application of fracture mechanics. In the computational mechanics community, Java is often considered as a simple tool to produce applications on the Internet and/or to effectuate computation on the network. For example, Java kind technologies are often used to couple and manage traditional codes written in C/C++/Fortran. This permits the developers to use ancient codes or part of code in coupled applications. A consequence of it is to keep a real computational efficiency. E.g. in [START_REF] Miller | Localized tensor-based solvers for interactive finite element applications using C++ and Java[END_REF], an interactive finite element application based on a coupled C++/Java is described. Comparative tests with FORTRAN and C are conducted on small problems using direct solvers based on tensor computations; this aims at illustrating the high efficiency computational potential of Java. In [START_REF] Marchand | A program for generating unstructured triangular meshes with a Java Graphics User interface[END_REF], the development of GUIs is put in prominent position on an unstructured mesh generator. In [START_REF] Padial-Collins | Parallel operation of CartaBlanca on shared and distributed memory computers[END_REF], the Java code CartaBlanca is presented. It is an environment for distributed computations of complex multiphase flows. Based on a finite volume approach, a solution scheme based on a Newton-Krylov algorithm is described. The code exhibits good performances. A similar environment has been developed to simulate electromagnetisms problems in [START_REF] Baduel | A parallel object-oriented application for 3-D electromagnetism[END_REF]. Both applications show the high potential of the approach to design more complex and general computational tools in mechanics for example. These developments exhibit the networking facilities provided in Java. A large number of publications shows the interest of Java and its efficiency: direct solution of linear systems [START_REF] Nikishkov | Comparison of C and Java performance in finite element computations[END_REF], FFT and iterative and direct linear systems solvers on Euler type flows [START_REF] Bull | Benchmarking Java against C and Fortran for Scientific Applications[END_REF], solution of Navier-Stokes flows [START_REF] Huser | A test suite for high-performance parallel Java[END_REF] and [START_REF] Riley | High-performance Java codes for computational fluid dynamics[END_REF]. Again, in all these papers, multiprocesses' management and networking capabilities of Java are put forward. In [START_REF] Bull | Benchmarking Java against C and Fortran for Scientific Applications[END_REF], [START_REF] Huser | A test suite for high-performance parallel Java[END_REF] and [START_REF] Eyheramendy | Object-oriented parallel CFD with JAVA[END_REF], performances of Java are tested on simple matrix/vector products. Compared to C/C++ code, only 20 to 30% of efficiency is lost. More recently, the description of a finite element code in Java was proposed in [START_REF] Eyheramendy | Advanced object models for mathematical consistency enforcement in scientific computing[END_REF]. The proposed design remains rather similar to the existing ones based on classical object-oriented approaches. Nevertheless, this paper shows that it is possible to develop a global code based on a pure Java approach. The latter can be applied to design large complex applications in computational mechanics taking into account complexity in modern computational mechanics: multiscale, multiphysics and multiprocesses applications. Firstly, enhanced and homogeneous data organization schemes to deal with complexity are proposed in Java. Secondly, developing a global application in a uniform environment including all the libraries needed is somehow an attractive idea for scientists and industrials. It is worth noting that similar approaches exist in C# (e.g. see [START_REF] Mackie | Object-oriented design of pre-conditionned iterative equation solvers using[END_REF], [START_REF] Mackie | Lessons learnt from using .NET for distributed finite element analysis[END_REF]). This brief analysis shows that even if performances do not achieve the one of C/C++/Fortran but get close to a certain extent, Java or similar approaches can be used in computational mechanics including CFD. These strategies offer the main advantage to design the simulation tools such a way that the technical strategies and the numerical algorithms are integrated in a seamless way.

Schwarz multiplicative multithreaded and distributed applications

The Schwarz multiplicative algorithm is described in [START_REF] Eyheramendy | Object-oriented parallel CFD with JAVA[END_REF]. In Fig. 1, the same decomposed domain is shown to be solved in a multi-threaded application or in a distributed application. Solving both problems is made through the same basic code for the physics and numerical algorithm. A description of the code goes beyond the scope of this paper and can be found in [START_REF] Eyheramendy | High abstraction level frameworks for the next decade in computational mechanics[END_REF], [START_REF] Eyheramendy | Advanced object-oriented techniques for coupled multiphysics[END_REF] and [START_REF] Eyheramendy | Advanced object models for mathematical consistency enforcement in scientific computing[END_REF]. The description of the multithreaded application is given in [START_REF] Eyheramendy | Advanced object-oriented techniques for coupled multiphysics[END_REF]. As shown in Fig. 1, the general structure of the code organization is the same. In the distributed version, the global Schwarz algorithm is managed through a thread located on the computer 1. The algorithm lies in two points: the management of the domains located on alternative computers and the communication schemes between these domains located on various computers. The global management is held by the way of distributed objects. The Java RMI package is used for that purpose. It permits the programmer to keep a natural and homogeneous organization of classes to deal with the domain decomposition algorithm. The consequence of it is to keep exactly the same structure as in both algorithms, managing real objects in the one case, and distributed objects in the other case. The later remains seamless for the programmer. The Schwarz algorithm is given in both cases as follows (programming details are omitted in the following code): In this piece of code, the equivalent domain is called using the code: listOfServer-Domain.get(indice0a) but the global scheme remains the same. In the first approach, the domain is located in the same memory space, as in the second approach, it is located on a separate computer. The lack of efficiency in communication schemes is the main drawback of the use of RMI distributed objects. The latter is not a problem for the management of the global Schwarz algorithm but produces a bottleneck when exchanging data. This is the reason why communications are programmed using the Socket Java class which a point to point communication scheme. This scheme is implemented at the level of boundary conditions object where the information provided by the neighbor domain id needed. Once again, only a local modification is needed within the global framework. For illustration purposes, both schemes are applied to a Navier-Stokes flow. For programming details in Java, the reader may refer to [START_REF] Flanagan | Java in a Nutshell[END_REF].

Numerical application

To illustrate both approaches, we study a Galerkin Navier-Stokes formulation stabilized by adding least-squares type terms. The equations of the problem and the formulation are presented in Fig. 2. Linearization of the problem is introduced in a Newton like scheme in the Schwarz multiplicative scheme. A direct linear system solver based on a Crout decomposition is used to solve the linear system at each iteration. Note, the same code is compiled once and run on all systems. The multithreaded version of code is run on a SGI Altix 16-processors Itanium2 1.3Ghz, 32 Go of RAM. The distributed version is run on a set of simple cluster of PC linked by a classical network. The code used for both applications is exactly the same. The global solution algorithm and the boundary conditions in charge of the communication scheme are of course not the same. The formulation is applied to the computation of a flow through a set of cylinders shown in Fig. 3. The domain is a periodic layer of cylinders. Numerical results are given in Fig. 3. First, the pressure contour on a typical cell is plotted. Secondly, the mean velocity computed overall the domain is given with respect to the gradient of pressure over the cell in the direction of the main flow. The latter represent the homogenized flow through the cylinders. For low velocities, the relation between the mean velocity and the gradient of pressure is linear. From a global point of view this can be assimilated to a Darcy's flow. But as far as the mean velocity is increasing the linearity disappears. The advection term is no more negligible and the global flow cannot be assimilated to Darcy's flow. It shows the influence of the advection term on the homogenized flow, that cannot any more considered as a linear Darcy's flow. Filtration laws could be established such a way. The same computation has been held using both approaches. From a practical point of view, we show that the same code can be run on heterogeneous systems. As both systems are different, we cannot compare the efficiency of both algorithms. This goes beyond the scope of our qualitative test. We show here that we can easily switch from a computer to another one without any problem using exactly the same code and without compiling it. This feature can be very interesting from an industrial point of view.

Conclusion

In this paper, we have presented two computational approaches based on the same numerical algorithm, a Schwarz overlapping domain decomposition method. The first approach is a multithreaded approach of the algorithm; the second one is a distributed version of it. In the first one, which is to be run on a share memory system, no communication between the domains is needed. The second one is based on two different communication schemes. The first one permits a master to manage the global Schwarz algorithm based on a classical object-oriented approach for distributed objects. Communications for data exchanges are implemented using classical sockets. This point to point communication scheme allows us to achieve efficiency. The implementation is held in Java which ensures the code to be run directly on all the systems. The same code is run for both applications. The advantages of such an approach are:

1. from a programmer point of view: a single language and a single approach for all the algorithms which means simplicity and reliability due to the fact that the finite element core remains identical. 2. From a user point of view: the same code can be used on heterogeneous systems depending on the availability of different systems.

We advocate that such an approach may simplify a lot the use of complex systems for single applications. This opens new tracks in the design of code that can be used in the context of either a shared memory system or a distributed memory system. In this context, mixing both approaches within a single application, taking advantage of the heterogeneous computers systems available at a given time is made possible in a simple way.

Fig. 1 .

 1 Fig. 1. Multi-threaded application versus distributed application.

Fig. 2 .

 2 Fig. 2. Navier-Stokes model. Initial-boundary value problem and stabilized finite elements formulation.

Fig. 3 .

 3 Fig. 3. Flow through cylinders. Computational domain, decomposed domain, numerical results: detail of iso-pressures, curve velocity/gradient of pressure (Experiment and simulation).