A Loyer

J.-J Sinou

O Chiello

X Lorang

Study of nonlinear behaviors and modal reductions for friction destabilized systems. Application to an elastic layer

Keywords: Squeal, non-linear dynamics, stability analysis, frictional contact, model reduction

As noise reduction tends to be part of environmental directives, predicting squeal noise generated by disc brakes is an important industrial issue. It involves both the transient and stationary non-linear dynamics of selfexcited systems with frictional contact. Time simulation of the phenomenon is an attractive option for reducing experiment costs. However, since such computations using full finite element models of industrial disc brake systems is time-consuming, model reduction has to be performed. In this paper, both the transient and stationary non-linear behaviors of the friction destabilized system and the effect of dynamical reduction on the non-linear response of a simple friction destabilized system are carried out. The first part provides a description of the general modeling retained for friction destabilized systems. Then, discretization and solving processes for the stability analysis and the temporal evolution are presented. The third part presents an analysis of a sliding elastic layer for different operating conditions, in order to better understand the non-linear behaviour of such systems. Finally, spatial model reduction is performed with different kinds of reduction bases in order to analyze the different effects of modal reductions. This clearly shows the necessity of including static modes in the reduction basis and that non-linear interactions between unstable modes are very difficult to represent with reduced bases. Finally, the proposed model and the associated studies are intended to be the benchmark cases for future comparison.

Introduction

Brake squeal is a complex phenomenon and a serious industrial issue. Many studies have been carried out in the automotive [START_REF] Ouyang | Numerical analysis of automotive disc brake squeal: a review[END_REF][START_REF] Kinkaid | Automotive disc brake squeal[END_REF][START_REF] Fritz | Investigation of the relationship between damping and mode-coupling patterns in case of brake squeal[END_REF][START_REF] Wagner | Minimal models for disk brake squeal[END_REF][START_REF] Herve | Analysis of squeal noise and mode coupling instabilities including damping and gyroscopic effects[END_REF] and aeronautics [START_REF] Gordon | Nonlinear model for aircraft brake squeal analysis: stability analysis and parametric studies[END_REF][START_REF] Chevillot | The destabilization paradox applied to friction-induced vibrations in an aircraft braking system[END_REF][START_REF] Sinou | Friction-induced vibration for an aircraft brake system-part 1: Experimental approach and stability analysis[END_REF][START_REF] Chevillot | Nonlinear transient vibrations and coexistences of multi-instabilities induced by friction in an aircraft braking system[END_REF] fields. However, until recently, little had been done regarding railway brake systems [START_REF] Lorang | Tgv disc brake squeal[END_REF][START_REF] Brizard | Performances of some reduced bases for the stability analysis of a disc/pads system in sliding contact[END_REF]. Squeal noise is closely related to friction and can be explained by two main kinds of approaches. The first is based on tribological considerations such as decreasing the friction coefficient with sliding speed [START_REF] Ibrahim | Friction-induced vibration, chatter, squeal, and chaos-part i: Mechanics of contact and friction[END_REF][START_REF] Ouyang | Friction-induced parametric resonances in discs: Effect of a negative friction-velocity relationship[END_REF]. The second is based on structural instability due to geometrical coupling of degrees of freedom [START_REF] Spurr | A theory of brake squeal[END_REF] or modes [START_REF] Fritz | Investigation of the relationship between damping and mode-coupling patterns in case of brake squeal[END_REF][START_REF] Jarvis | Vibration induced by dry friction[END_REF][START_REF] Coudeyras | Periodic and quasi-periodic solutions for multi-instabilities involved in brake squeal[END_REF] in the system through the frictional contact. In this paper, the instability of the operating point due to mode coupling is considered. Thus, frictional contact will be described by using a Coulomb law with a constant coefficient.

In order to overcome the computational cost of a transient dynamics simulation on a friction destabilized system, one solution is to perform model reduction. However, the impact of reduction on non-linearities is not very well known and thus great care is needed when using such a method. The study has two main objectives. This paper first focuses on the behaviour of a simple self excited system under different operating conditions and frictional states. Both transient-state and the stationary-state behaviours are studied with different system parameters and initial conditions. Secondly, different reduction bases for approximating the non-linear selfexcited response of the system are presented. The impact of the reduction on the non-linear self-excited solution is studied.

This paper is divided into six sections. The first presents the general model retained for the frictional contact problem in view to studying self-excited systems with a frictional interface. The second deals with the solvers used for both the stability analysis and the transient non-linear dynamics. The third section is dedicated to the presentation of the different tools used for analysing the transient and stationary self-excited non-linear solutions. Then the stability of the system retained is examined for different friction coefficients, after which the transient and stationary non-linear behaviours of the system for different initial conditions and operating points are studied without reduction. Finally, the efficiency of the different modal reduction bases is demonstrated for different operating points and physical parameters.

Modelling

Contact laws and description of the non-linear interface

Friction induced instabilities arise when two bodies slide against each other. Tribological reasons can be invoked in order to explain the occurrence of instability (e.g. variation of the friction coefficient with sliding speed [START_REF] Ibrahim | Friction-induced vibration, chatter, squeal, and chaos-part i: Mechanics of contact and friction[END_REF]) as can be structural causes (e.g. coupling modes due to friction forces [START_REF] Coudeyras | Periodic and quasi-periodic solutions for multi-instabilities involved in brake squeal[END_REF]). The structural approach is adopted here since it allows simpler modeling of the contact interface with unimpaired representation of structural dynamics phenomena. Thus a Coulomb law with a constant friction coefficient µ is used:

   r t ≤ -µr n r t = -µr n ⇒ ∃ξ ∈ R + , ut -v g = -ξr t r t < -µr n ⇒ ut -v g = 0 (1)
where r is the contact reaction, u is the displacement field, v g is the Eulerian sliding speed, dot defines the time derivative, n and t subscripts stand for the normal and tangential projection of a field on the contact interface respectively. This law states that the tangential component of the contact reaction cannot be greater than a limit that is proportional to the normal component. Moreover, if this limit is reached, there is a sliding phenomenon in the opposite direction of the tangential friction force, otherwise there is no sliding. To deal with the unilateral contact, a non-regularized Signorini law was chosen:

   u n -g ≤ 0 r n ≤ 0 (u n -g) r n = 0 (2)
where g is the initial gap at the contact interface. The advantage of the Signorini law is that it does not require the introduction of a coefficient such as contact stiffness that would require measurement. Indeed using a contact stiffness might be interesting for accurate modeling at the frictional interface if valuable experimental tests were available. However, the use of a non-regularized Signorini law is sufficient for the main purposes of this paper which illustrate different global nonlinear behavior and the effect of model reduction for the elastic layer with frictional interface. The Coulomb law with a constant coefficient also means that only one coefficient has to be measured rather than characterizing a complex law. If the two bodies have a uniform rectilinear relative movement in translation, then it is easy to find the equation in an Eulerian framework. If the movement is not a uniform translation, neglecting the convective terms allows obtaining simple equations in which the sliding speed only appears in the contact laws. This formulation is still valid for dealing with contact between two bodies by considering the relative displacement, velocity and contact reaction. The deformable body in contact is considered to be a linear visco-elastic continuous medium under the small perturbation hypothesis, such that the only non-linearity considered in the system is the frictional contact interface (i.e. contact, no-contact, sliding and adhesion states). As illustrated in Figure 1(a) the system is in a domain Ω with boundary ∂Ω admitting {∂Ω u , ∂Ω f , S c } as partition, subjected to a surface load f and a volume load f v . ∂Ω u is the part of ∂Ω where displacement u d is prescribed, ∂Ω f is the part of ∂Ω where surface load f s is applied and S c is the potential contact interface with initial gap δ. The dynamics of the system is then given by the principle of virtual power:

                     find u ∈ U and r such as ∀ u⋆ ∈ U 0 Ω ρ u⋆ üdΩ + Ω ε⋆ (B ε + Aε) dΩ = ∂Ω f u⋆ • f s dS + Ω u⋆ • f v dΩ +
where B and A are respectively the viscous stress tensor and the elastic stress tensor. In this equation, both the kinematical fields and the contact reactions are unknown. Any linear law could be considered in the previous modeling. However, in the practical example discussed in Section 5, an isotropic elastic linear material with Rayleigh damping will be used. In order to simplify modeling, the sliding speed is taken to be an Eulerian sliding speed, meaning that it only appears in the friction equation by modifying the relative speed between the bodies. This allows considering the fixed configuration of the bodies. Finally, the contact and friction laws can be rewritten in terms of projections on the negative real set (proj R -) and on the Coulomb cone (proj Kµ ) [START_REF] Jean | The non-smooth contact dynamics method[END_REF]. This reformulation is only used to facilitate the numerical implementation in the treatment of the contact state. By doing this, the frictional contact laws (1) and [START_REF] Kinkaid | Automotive disc brake squeal[END_REF] give the following equations:

r n = proj R -(r n -ρ u n (u n -g)) , ∀ρ u n > 0 where proj R -(x) = min(x, 0) (4)
r t = proj Kµ (r -ρ t ( ut -v g )) , ∀ρ t > 0 with proj Kµ (x) = x t if x t /|x n | ≤ µ µ |xn| xt x t otherwise (5)
where ρ u n and ρ t are two arbitrary positive scalars called normal displacement augmentation parameter and tangential augmentation parameter respectively.

Stability analysis of the self-excited system with a frictional interface

If the boundary conditions are constant in time, the problem defined in Equation (3) has a quasi static solution (u e , r e ) that verifies: It may occur that this quasi static operating point is unstable. In order to perform a stability analysis of this solution, it is necessary to linearize the equations governing the non-linear dynamical behaviour around the equilibrium point. Restriction on the perturbation is imposed such that it keeps the contact state unchanged (i.e. no normal displacement is allowed on the effective contact zone). Perturbation of the quasi static solution is sought in the form (u, r) = (u e + û, r e + r) where the . quantities represent the perturbation. After linearization around the equilibrium point, the variational equation satisfied by (û, r) is given by:

                       Ω ε⋆ Aε e dΩ = ∂Ω f u⋆ • f s dS + Ω u⋆ • f v dΩ + Sc u⋆ • r e dS, ∀ u⋆ ∈ U 0 u e n ≤
                         find û ∈ U and , r such as ∀ u⋆ ∈ U 0 Ω u⋆ ρ üdΩ + Ω ε⋆ B ε+ Aε dΩ =
By searching harmonic solutions (û, r) = e λt ( û0 , r0 ) to Equation ( 6), the variational equation becomes:

           find û0 ∈ U, r0 and λ ∈ C such as ∀ u⋆ ∈ U 0 λ 2 M (û 0 , u⋆ ) + λ Ĉ (û 0 , u⋆ ) + K (û 0 , u⋆ ) = S e c r0n u⋆ • (n -µt) dS û0 • n = 0 on S e c (7)
where M, C, K are respectively the mass, damping and stiffness operators of the structure without contact and

Ĉ (û 0 , u ⋆ ) = C (û 0 , u ⋆ ) -S e c µ r e n vg û0b u ⋆ b dS.
The additional part in the damping operator Ĉ is due to the linearization of the sliding direction for 3D-cases that gives the term along b in Equation [START_REF] Gordon | Nonlinear model for aircraft brake squeal analysis: stability analysis and parametric studies[END_REF]. The last condition in Equation (7) (i.e. û0 • n = 0) can be imposed by searching û0 in W 0 and the right-hand side of the equation vanishes for u⋆ taken in W µ . This leads to the following reformulation:

       find û0 ∈ W 0 and , λ ∈ C such that ∀ u⋆ ∈ W µ λ 2 M (û 0 , u⋆ ) + λ Ĉ (û 0 , u⋆ ) + K (û 0 , u⋆ ) = 0 W 0 = {u, u • n = 0 on S e c } and W µ = {u, u • (n -µt) = 0 on S e c } (8)
The latter equation is a non symmetric quadratic eigenvalue problem (QEP). Finite element discretization can be used to obtain non-symmetric mass, stiffness and damping matrices M, K and C from the classical finite element matrices M, K and C with the use of T N and T µ matrices, which are respectively the bases of the fields orthogonal to the normal to the contact and orthogonal to the direction n -µt on the contact interface. For 3D-cases, taking into account the extra damping term due to the linearization of the sliding direction for computation Ĉ gives:

Ĉ = C -µP B T DP B with D ii = (r e n ) i (v g ) i (9)
where (v g ) i is the norm of the sliding speed v g at the i th node, and P B the projection matrix on the second tangential direction b at contact nodes.

The quadratic eigenvalue problem to be solved is then given by:

λ 2 M + λ C + K u 0 = 0 with M, C, K = T µT M, Ĉ, K T N ( 10
)
where T N and T µ are respectively the basis of the finite element approximated W 0 and W µ . Classical linearization of this quadratic eigenvalue problem (QEP) leads to the generalized eigenvalue problem (GEP):

Λ Ã + B Φ 0 = 0 with à = C M M 0 and B = K 0 0 -M (11)
where Φ 0 corresponds to mode shapes in state space. Taking the half first or last rows of Φ 0 gives the mode shapes in physical space with the bilateral contact constraint, and these can be transferred to physical space without constraint by multiplication by T N . Solving this problem is achieved by using the Residual Iteration Method [START_REF] Bobillot | Iterative computation of modal sensitivities[END_REF]. Among the eigenvalue solutions to the problem, those with a positive real part correspond to unstable modes. The divergence rate of a mode is defined as 2 ℜ(λ) ℑ(λ) , which corresponds to a negative damping rate (with ℜ(λ) and ℑ(λ) are the real and imaginary parts of the associated unstable mode). Perturbation of equilibrium leads to increase vibrations on the unstable modes. The growth of these vibrations is limited by non-linear events. In the current model, since the only non-linearity accounted for is frictional contact, adhesion and separation are the only phenomena capable of limiting increase.

Solving

Fixed point statement

For the static and dynamic non-linear problem, the contact laws can be interpreted as a fixed point statement:

r = f (r, ρ u n , ρ t ) (12)
where the function f is defined by:

f : (r, ρ u n , ρ t ) → proj Kµ r t +n proj R -r n -ρ u n (u n (r) -g 0 ) -ρ t ( ut (r) -v g ) +n proj R -r n -ρ u n (u n (r) -g 0 ) (13
) This fixed point statement defines an iterative solver that applies f to a given initial guess for the contact reaction r and iterates until a stop criterion is reached. Application of f requires solving the principle of virtual power given in Equation (3) with the current estimation of r. This type of algorithm requires that function f is contractive. In order to ensure contraction, the ρ u n and ρ t parameters of function f must be optimized. Using classical finite element discretization of the problem with linear elements on the potential contact zone leads to the following equation:

Mü + C u + Ku = f + r c (14)
with (u, u, r c ) verifying Equations ( 4) and ( 5) at each mesh node. f and r c are the generalized force and contact reaction respectively and M, K and C are the classical mass, stiffness and damping matrices of the structure. The contact laws are verified at each node with the discretized displacement and the generalized contact reaction, meaning that the fixed point statement is true when applied at each node. This allows defining two algorithms for solving both the static and the dynamic problems.

Transient non-linear temporal response computation

For the computation of the transient solution a first order θ-method time integration scheme is used with a time step τ . The resulting integrator matrix is denoted J. θ is taken as equal to 0.5 in order to avoid numerical damping. Avoiding numerical damping is especially important since it is well known that damping has critical effects on self-excited systems. Equation ( 4) is then reformulated in:

r i+1 n = proj R -r i+1 n -ρ n ġi+1 n + bi τ , ∀ρ n > 0 with bi = u i n -g 0 + τ (1 -θ) ui n ( 15
)
so that the only kinematic unknown dealt with is the velocity field. Term b is an equivalent gap computed with a modified integration scheme [START_REF] Jean | The non-smooth contact dynamics method[END_REF] in order to impose inelastic shock at the expense of additional damping. Detailed information on the integration scheme, computation of the equivalent gap and solving algorithm can be found in Appendix A. Evaluation of the function f given in Equation ( 13) requires two projections: one for the normal component and one for the tangential one. This involves two numerical augmentation parameters ρ n and ρ t . Optimization of these parameters can be considered as one of the most important step for efficiency of the numerical scheme. It is performed by minimizing the Lipschitz constant of the application. For the sake of simplicity, it is assumed that minimizing the Lipschitz constant for each projection is equivalent to minimizing it to chain both parameters. Usually, the latter are chosen as a global scalar related to the smallest eigenvalue of the integrator matrix [START_REF] Khenous | Hybrid discretization of the signorini problem with coulomb friction. theoretical aspects and comparison of some numerical solvers[END_REF]; however, here they are taken as diagonal operators, in order to allow finer optimization. Optimization of ρ n and ρ t thus leads to the use of diagonal approximations of the integrator matrix J on the normal and tangential contact degrees of freedom. In this study extraction of the diagonal of the condensation of the integrator matrix on the contact degrees of freedom was used. To facilitate the reader's comprehension, the static solver can be derived from this algorithm by using the stiffness matrix instead of the integrator scheme matrix and computing the tangential reaction in the sliding direction while imposing the sliding condition. By using these solvers, it is possible to compute both the static equilibrium of sliding systems with the associated stability analysis and perform the temporal simulation of the system. The stability analysis can also be used to validate the results of the temporal solver with those of the linearized system around an equilibrium point before the occurrence of strong non-linear events. It also offers the possibility of initiating a system on one or more unstable modes, in order to reduce the time in which the final non-linear state of the self-excited system is reached.

Analysis of the solution

In this section, the different tools that will be used to provide a good description of non-linear contributions during the non-stationary transient response of the frictional system are presented. Firstly, the time-frequency features and a tool based on the Continuous Wavelet Transform (CWT) are discussed. Secondly we propose the use of projection on the complex modes and peak-to-peak or peak-to-mode M-MAC to perform the analysis of the evolutions of the mode shapes and their contributions.

The continuous wavelet transform

For the time frequency analysis, calculation of the frequency components of the whole signal is sufficient for stationary signals. It can be easily obtained with a fast Fourier transform (FFT). In the case of non-stationary signals, the time-dependence of the frequency components is essential for better understanding of transient nonlinear behaviour. One classical method for this type of analysis is the short-time Fourier transform (STFT), first proposed by Gabor in [START_REF] Gabor | Theory of communication[END_REF], which is appropriate for weak non-stationary signals. For stronger non-stationary signals, the use of continuous wavelet transform (CWT) is a good alternative to the STFT [START_REF] Morlet | Wave propagation and sampling theory-part ii: Sampling theory and complex waves[END_REF]. In the following section, only the main principles are stated.

The first step of a wavelet analysis is the choice of a mother wavelet, φ(t), which is a very fast decaying function that fulfills the admissibility conditions defined in [START_REF] Goupillaud | Cycle-octave and related transforms in seismic signal analysis[END_REF]. A family of analysing functions is created from the mother wavelet by translation and dilatation according to Equation [START_REF] Herve | Extension of the destabilization paradox to limit cycle amplitudes for a nonlinear self-excited system subject to gyroscopic and circulatory actions[END_REF] where a is the scale -or dilatationfactor and b is the time translation factor. Multiplication by 1 √ a ensures energy normalization across the different scales. The CWT is then performed by projecting the signal on the wavelet family obtained with Equation ( 16), producing the CWT coefficients given by Equation [START_REF] Coudeyras | Periodic and quasi-periodic solutions for multi-instabilities involved in brake squeal[END_REF] where ⋆ denotes the complex conjugate.

φ a,b (t) = 1 √ a φ( t -b a ) ( 16
) CWT s φ (a, b) = ∞ -∞ s(t)φ ⋆ a,b (t)dt (17)
The wavelet coefficients represent a measurement of the correlation of signal s(t) with the family of analysis. Thus CWT s φ (a, b) expresses the local information concerning s(t) at scale a and time b. A time-scale representation of the CWT is obtained by plotting the coefficients over time and scale. In practice, the scale is linked to the inverse of a frequency, thus a time-frequency representation can be achieved by plotting over time and inverse scale. The exact correspondence between scale and frequency depends on the mother wavelet. The wavelet power spectrum defined as | CWT s φ (f, b)| 2 is usually preferred to the wavelet time-frequency plot due to its analogy with the Fourier power spectrum. Many different types of mother wavelets adapted to different specific cases have been developed [START_REF] Mesa | Adapted wavelets for pattern detection[END_REF][START_REF] Du | Crackle detection and classification based on matched wavelet analysis[END_REF]. In the following, one of the most widely used mother wavelets, the Morlet wavelet, has been chosen. It can be expressed as:

φ(t) = e jk0t e -(1/2)t 2 F T → φ(ω) = √ 2πe -(1/2)(ω-k0) 2 (18)
For the Morlet wavelet, the relation between scale and frequency is given by the central frequency ω c = k 0 /a at which φ(aω) is maximal. Factor k 0 is called the wavenumber. It is important to note that the Morlet wavelet does not respect the admissibility condition, since φ(0) = 0, which means that the inverse transform is not properly defined. Due to the proximity of the Morlet wavelet with the exponential function, the Morlet CWT can be seen as an STFT with a short window at high frequency and a long window at low frequency. This adaptive window size partially overcomes the time-frequency resolution limitation of the STFT.

Complex modal projection

Another tool used in order to analyse the transient non-linear response is the complex modal projection. The evolution of the system at the beginning of the time interval is given by the complex mode shapes and eigenvalues computed around the equilibrium point (i.e. the modes obtained via the stability analysis). Then, the evolution of the participations of these modes, especially the unstable ones, during the transient solution can be behaviourly useful for studying the non-linear behaviour of the self-excited system and the associated physical phenomenon more effectively.

Since the complex modes obtained via the stability analysis have bilateral contact constraints, the displacement and velocity fields are projected on the subspace of the fields with zero normal displacement by using T N and T J matrices, (i.e. the bases of the fields with no normal displacement on the contact zone and of the normal displacement of the contact zone). The matrices T N and T J verify T N T T J = 0. Thus the displacement u can be written as u = T N u w + T J u k where u w represents the component with no normal displacement and u k the normal gap. The component u w can then be extracted from the whole field using:

u w = H -1 T N T u = T N T T N -1 T N T u (19)
Considering a temporal evolution of the system given by (u, u), the corresponding state space perturbation evolution α is given by:

α = H -1 T N (u -e u) H -1 T N u = u w uw (20)
By introducing the adjoint eigenvalue problem (i.e. by using ÃT and BT matrices defined in Equation ( 11)) we obtain adjoint mode shapes Ψ 0 . The solution Φ 0 of Equation [START_REF] Brizard | Performances of some reduced bases for the stability analysis of a disc/pads system in sliding contact[END_REF] and the adjoint mode shapes Ψ 0 verify the bi-orthogonality relation:

∀i = j Ψ 0i T ÃΦ 0j = 0 Ψ 0i T ÃΦ 0i = 0 (21)
α admits a unique decomposition on the mode shapes:

α = n k=1 β k (t)Φ 0k (22)
For the reader comprehension, β k (t) can be expressed by γ k e (r k +jω k )t where r k and ω k are close to the real and imaginary parts of the eigenvalue of the complex mode (at the beginning of the time interval). During the onset and the evolutions of non-linear events, r k (t) tends to 0 (or to a null-mean function) and ω k (t) can be different from the initial frequency and oscillate. It can be noted that even if the mode shapes are complex, the displacements are real. Indeed, if one complex mode is present in the decomposition, its complex conjugate is also present. So, the complete result is therefore real. Thus when seeking the decomposition of α given by Equation [START_REF] Morlet | Wave propagation and sampling theory-part ii: Sampling theory and complex waves[END_REF], left multiplication by the adjoint modes allows the computation of the modal participation β j with the following formula:

β j = Ψ 0j T Ãα Ψ 0j T ÃΦ 0j (23)
If both direct and adjoint modes are normalized with respect to mechanical energy, the representation of β j versus time gives information on the energy contribution of the j th mode in the temporal solution. This evolution is given in a plot presenting log(|β j |) versus time. It is important to note that since the complex modes are coupled, the summation of the different modal energies is greater than the total energy of the system.

u w (t) = n k=1 a k γ k e λ k t = n k=1 a k γ k e r k t e jω k t (24)
where a k is the first half of row Φ k and γ k depends on the initial conditions. It can be noted that the complex modal projection gives additional information with physical and modal meaning compared to CWT or FFT analysis that are more signal processing tools. Indeed, one of the most useful piece of information in the complex modal projection is to detect the contribution of one or more complex modes computed through the stability analysis in the transient and final non-linear solution of the system This can be an efficient tool to compare and correlate the non-linear response with the results of the stability analysis. Finally, another advantage of this tool is to detect the onset of the non-linear stationary signal (i.e. periodic or quasi-periodic response of the system) through evolutions of the complex modal projection: indeed, when the limit cycle is reached, contributions of all the complex modal projections are constant or oscillate around an average value.

Peak-to-mode and peak-to-peak M-MAC

In order to compare the shapes of the stabilized solution at various spectral peaks with shapes of a solution computed on another reduction basis or with the complex mode shapes obtained via the stability analysis, an

Parameter value E 5 10

9 M P a h 0.05 m L 0.15 m mesh 50 × 20 elements ν 0.2 ρ 2500 kgm -1 α 7.8 × 10 -1 β 2 × 10 -7
Table 1: System parameters extension of the classical Mass normalized Modal Assurance Criterion (M-MAC) is proposed. This criterion will be referred to as peak-to-peak M-MAC when it compares spectral peaks of solution computed on different bases and as peak-to-mode M-MAC when used to compare peaks with mode shapes. The M-MAC and this extension is presented to aid the reader's comprehension. The expression of the M-MAC matrix used in the present study is

M-MAC (A, B) = |A ⋆T MB| (25)
with a normalization of each vector of A and B versus the mass matrix of the system. This M-MAC is a direct extension of the complex mode NCO reviewed by Morales in [START_REF] Morales | Comments on the mac and the nco, and a linear modal correlation coefficient[END_REF]. For peak-to-peak M-MAC, A and B are extractions of the spectral peak shapes of a reference solution and of the solution computed on a reduced basis.

For peak-to-mode M-MAC, B is the matrix of the complex mode shapes of the system.

Description of the system and stability analysis

In this section, the self-excited system is defined and the parameters and operating conditions are presented. Secondly, the stability analysis is performed for different friction coefficients.

Description of the system

The system considered here is a homogeneous 2D elastic layer of length L and height h, sliding on a rigid plane, as described in Figure 1(b). The layer is discretized by using a 50 by 20 linear element FE mesh. Its upper boundary is fixed and the two lateral sides have a periodicity constraint. The layer has an imposed displacement δ towards the rigid plan. The friction interface S c between the layer and the plane is characterized by the constant friction coefficient µ and the Eulerian velocity of the plane v g . The damping in the structure is taken as a Rayleigh damping proportional to the mass and stiffness operators such that the damping operator is given by C = αK + βM. The physical values of the system parameters (L, h, the density ρ, the Young modulus E, the Poisson coefficient ν, the Rayleigh coefficients α and β) are given in Table 1. The operating parameters are the sliding speed of the plane v g , the friction coefficient µ and the imposed displacement δ.

Stability analysis

First of all, the stability analysis is performed for different friction coefficients from µ = 0 to µ = 1. In Figure 2, it appears that the system is stable for µ = 0. Then, increasing the friction coefficient increases the number of unstable modes. For example, for µ = 0.19 (µ = 0.35 and µ = 0.5 respectively), five (nine and twelve respectively) unstable modes appear. Moreover, it can be observed that the divergence rate of the unstable modes increases as the friction coefficient increases. It should be noted that the contact state of the nodes of S c in the quasi static solution is independent from the three operating parameters defined previously. This implies that the modes and eigenvalues calculated via the stability analysis depend only on the evolution of the friction coefficient µ. Figure 3 shows the mode shapes of the twelve unstable modes at µ = 0.5. The complex mode shapes presented here are actually traveling waves. Most unstable modes are left traveling (LT) waves (i.e. opposed to the sliding speed), however some of them are right traveling (RT) waves.

(c) (d) (e) (f) (g) (h) (i) (j) (k) (l)

Transient analysis of system behaviour

As previously explained in [START_REF] Sinou | Transient non-linear dynamic analysis of automotive disc brake squeal -on the need to consider both stability and non-linear analysis[END_REF], different stationary or quasi-periodic solutions (i.e. limit cycles) can be observed due to a small evolution of the initial conditions. Therefore in this section, we investigate the transient behavior of the system with fixed parameters and various functioning and initial conditions.

Presentation of case studies

All the cases that are presented in the following sections have been chosen to give an overview of the different transient and stationary non-linear behaviors that can be obtained for different physical parameters and operating conditions. The aim of study is to obtain better understanding of the physical phenomena and show the possible links between a stability analysis and the temporal solution.

For the reader comprehension, the non-linear behavior of system with only one unstable mode will not be presented in order to condense the results presented in this paper and only focus on the most original results. Indeed, in this case (one unstable mode), the non-linear response has a simple signature with only the participation of the unstable mode and its harmonics. The interested reader can refer to the papers [START_REF] Chevillot | Nonlinear transient vibrations and coexistences of multi-instabilities induced by friction in an aircraft braking system[END_REF][START_REF] Sinou | Transient non-linear dynamic analysis of automotive disc brake squeal -on the need to consider both stability and non-linear analysis[END_REF] for more details.

In order to cover a reasonable range of situations, two values for the friction parameters were chosen (µ = 0.19 and µ = 0.35). They correspond to cases with 5 and 9 unstable modes that are presented in the Sections 6.2 and 6.3, respectively. They have been selected to undertake different non-linear transient and stationary dynamical behaviors that can be more complex with contributions from all the unstable modes, their harmonics and combinations. For both cases with 5 and 9 unstable modes, the two other parameters δ and v g are used to choose between separation and adhesion as first occurring non-linear events (i.e. the first nonlinear phenomenon that appears in the transient response). It can be observed that this terminology only means that either adhesion or separation occurs before the other in the transient solution. However there is no possibility to ensure that after this first event, the other type of non-linear phenomenon (separation or adhesion) does not occur. Thus two main functioning conditions are defined: the first is referred to as adhesion configuration (with δ = 1 10 -5 m and v g = 0.1 ms -1 ) and the second as separation configuration (with δ = 2.5 10 -6 m and v g = 0.5 ms -1 ). The description of the different case are given in Table 2.

Moreover, five different initial conditions were tested for each case. They are defined as follows: a nonzero initial condition is given for

• only the most unstable mode • only the least unstable mode • the most unstable and a mildly unstable mode have no zero • the least unstable and a mildly unstable mode • all the unstable modes where the term "most" (or least) unstable mode" stands for the mode with the most (or least) important divergence rate. These initial conditions are summarized in Table 2 by indicating the selected modes for initial conditions.

In the following parts of the paper, the following notations will be used for the test cases: case X IT defines the case where X is the number of unstable modes, I are the modes present in the initial condition (i.e. then modes with a nonzero initial condition) and T is the type of first occurring non-linear event (A for adhesion or S for separation). For example, case 5 9S is a case with 5 unstable modes and an initial condition only on mode 9 with separation as a first occurring non-linear event, and case 5 9×33A is a case with 5 unstable modes with modes 9 and 33 as the initial condition and adhesion as the first occurring non-linear event. For a condition on all unstable modes, I is denoted all.

Since there is a periodicity condition, all nodes of the contact interface are kinematically equivalent. Thus in the following figures that represent displacements or velocities on the contact interface, the observation node will not be specified. Figures 4(a) and 4(b) illustrate the evolutions of tangential velocity and displacement at the frictional interface. One can see that a first "cycle" occurs very early. The phase plot of this first cycle is represented in dashed black lines in Figure 4(c). This cycle has a fundamental frequency (f 1 = 16.3 kHz) very close to the frequency of mode 9, as shown in Figure 4(d). It can also be seen that the harmonic components of orders 2, 3, 4 and 5 are visible but with rather slight contributions. However, a transition can be observed at t = 0.09 s, with a surge in tangential velocity (see Figure 4(a)) and a significant drop in the mean tangential displacement (see Figure 4(b)). This transition appears to be very fast with a brutal variation of the global equilibrium position of the system (i.e. variation of the mean value of displacement as indicated on Figure 4(b)). This reduction of the average tangential displacement can be related to a reduction of the apparent global friction coefficient due to local adhesion. Analysis of the phase plot of tangential component (see Figure 4(c)) also clearly shows the two different cycles.

Case

As indicated in Figure 4(d), the fundamental frequency of the second cycle (f 0 = 14.4 kHz) is very close to the frequency of mode 6 (i.e. 14.25 kHz). The harmonic components of orders 2, 3 and 4 are also present. However, the previous participation of mode 9 in fundamental and harmonic frequencies is not visible. Finally, combinations of the fundamental frequencies f 0 and f 1 appear clearly during the transition state (t ∈ [0.09; 0.11] s). Complex modal projection (see Figure 4(e)) then provides a clear view of the evolution of the non-linear transient and stationary solutions. It can be seen that during the transition phase, the contribution of mode 9 (denoted M9 in Figure 4(e)) decreases whereas the contribution of mode 6 (denoted M6 in Figure 4(e)) increases. Moreover, the observation made by Lorang [START_REF] Lorang | Stability and transient analysis in the modeling of railway disc brake squeal[END_REF] that the predominant unstable mode in the final cycle is that capable of exhibiting the higher energy contribution is verified. Finally, as illustrated in Figure 4(f), the main peak in the FFT of the final state is very close to mode 6.

In this case, one of the shortcomings of time integration for squeal analysis is illustrated: it is difficult to ensure that the final state is reached. Integration could have been stopped during the first cycle. Overcoming this issue can be achieved by using methods such as constrained harmonic balance method [START_REF] Herve | Extension of the destabilization paradox to limit cycle amplitudes for a nonlinear self-excited system subject to gyroscopic and circulatory actions[END_REF][START_REF] Coudeyras | Periodic and quasi-periodic solutions for multi-instabilities involved in brake squeal[END_REF][START_REF] Coudeyras | A new treatment for predicting the self-excited vibrations of nonlinear systems with frictional interfaces: The constrained harmonic balance method, with application to disc brake squeal[END_REF] which are able to predict a periodic or quasi periodic final state, but they do not provide information on the transient behavior of the non-linear solution. Transient analysis allows visualization of the different phases of the solution, sometimes showing that the amplitude during the transient phase is higher than in the stabilized final state. Moreover, using complex modal projection to analyze the time integration results allows effective detection of the final state: as indicated in Figure 4(e), evolutions of the complex modal projection are observed during the transient non-linear response of the system (see for example the growth rates for modes 3 and 6). When the limit cycle is reached, contributions of all the complex modal projections are constant or oscillate around an average value. Finally, it is observed that the growth rate of all the instabilities of modes 3, 6 and 11 are approximately the same. This fact illustrates that evolutions of the real part for the unstable complex non-linear modes (i.e. the real part of eigenvalues of the Jacobian matrix for a given time) are approximately the same. This may be consistent with the values of the divergence rate of the corresponding unstable modes that are almost identical via the stability analysis (see Figure 2(b)).

Case 5 9S

This case represents the case with 5 unstable modes (i.e. µ = 0.19) and an initial condition only on mode 9 with separation as a first occurring non-linear event.

Figures 5(a), 5(b), 5(c)) and 5(d) represent the velocities and phases plots for the normal and tangential components respectively. In Figures 5(a) and 5(b), a succession of three phases is shown: the first two correspond to an increase of the separation phenomenon, with an increase of the normal component on the contact zone. The third and last phase is characterized by saturation and stabilization through adhesion: the normal velocity falls to a very low level in comparison with the first two phases when the tangential component reaches its maximum level. Here again, due to the adhesion, the reduction of the global friction coefficient leads to a reduction of the average tangential displacement in the third phase. This succession of phases is also visible on the phase plots (see Figures 5(c)) and 5(d)) where the amplitude of the final tangential component cycle is greater than that of the first cycle. Moreover, the final cycle for the normal component has a very small amplitude.

As indicated in Figure 5(e), the fundamental frequency of the final stationary solution, which is equal to f 2 = 10.45 kHz, is very close to the frequency of mode 3 (i.e. 10.52 kHz). The higher harmonics, especially order 3, make a significant contribution. The evolutions of the participations of the unstable mode frequencies, their harmonics and combinations during the transient and stationary non-linear behavior are also observed. As explained previously for case 1 (Case 5 9A), the first and third phases have a simple non-linear signature with the major contribution of one unstable mode frequency (modes 9 and 3 respectively) and its harmonics. During the second phase, the non-linear behavior is more complex with the contribution of three unstable mode frequencies (modes 3, 6 and 9), their harmonics and combinations (see for example f 0 -f 1 at 3.79 kHz). Then, on the modal projection illustrated in Figure 5(f), quite complex competition in the second phase is observed: mode 9 is saturated while mode 6 reaches its saturation around 0.015 s. This saturation almost stabilizes the system. It can also be noted that mode 3 increases in the second phase and finally reaches saturation through adhesion in the third phase. Finally, Figure 5(g) indicates the modal identification for the main peak of the non-linear stationary solution. In this case, the shape appears to be very similar to that of mode 3. All these physical comments are in perfect correlation with the previous observations made for the CWT (in Figure 5(e)).

Influence of initial condition

Finally, the influence of different initial conditions (previously described in Section 6.1 and in Table 2) is tested for all cases. In both adhesion and separation cases, the limit cycle obtained was not dependent on the initial conditions for the five initial conditions tested and both adhesion and separation cases. Figures 6(a) and 6(b) illustrate the convergence of the modal projections on mode 3 for example. It also shows that the time before the cycle is reached greatly depends on the initial condition, but that some parts of the transient signal can be recognized as a "translation in time" from another initial condition. This case represents the case with 9 unstable modes (i.e. µ = 0.35) and an initial condition only on mode 9 with adhesion as a first occurring non-linear event.

This case exhibits one of the simplest non-linear behaviors for the system with 9 unstable modes. The initial condition on mode 9 leads to the saturation of this mode and this state gives the limit cycle. Figure 7(a) shows that there is a transient surge both in displacement when adhesion saturation occurs around t = 5 10 -4 s. After this surge, amplitude continues to grow until limit cycle is reached. The phase plot on Figure 7(b) clearly shows that the limit cycle is mono-harmonic. This fact is confirmed by the CWT analysis on Figure 7(c). However, the complex modal projection given in Figure 7(d) shows significant participation of mode 33. Moreover, no participation of frequencies of other unstable modes in the transient phase is observed. To facilitate comprehension, it should be noted that mode 33 has the double frequency of mode 9 and thus its participation may be a response to forced excitation by harmonic of order 2 of mode 9. Peak to mode M-MAC (see Figure 7(e)) shows that the main peak (at f 0 = 16.16 kHz) of the response is spatially close to mode 9. The second main peak of the stationary non-linear response (at 2f 0 ) is spatially very close to mode 33. The limit cycle obtained with 9 9A is also obtained with all the tested initial conditions 9 XA except 9 allA.

Case 9 allA

This case represents the case with 9 unstable modes (i.e. µ = 0.35) and an initial condition on all unstable modes with adhesion as a first occurring non-linear event.

This is a special case in the 9 XA series. Among the five initial conditions tested, it is the only one having a different final state. Figure 8(a) illustrates the tangential displacement at the frictional interface. A complex transient non-linear behavior is observed which is characterized by beats produced by the presence of different fundamental frequencies. Figure 8(b) illustrates the phase plot for the whole transient response. The last time steps in the phase space indicates two loops around two different positions, when the trace of the phase plot between t = 5 ms and 10 ms is composed of one single loop slowly shifting to the right. Indeed, as shown in Figure 8(c), the final limit cycle is quasi-periodic with two fundamental frequencies (f 0 = 14.3 kHz and f 1 = 10 kHz) and the presence of their harmonics (2f 0 and 3f 0 ) and the combinations of frequencies (f 0 + f 1 and f 0 -f 1 ). The two fundamental frequencies f 0 and f 1 are close to the frequencies of modes 6 and 3 respectively. Moreover, CWT analysis (Figure 8(c)) shows a progressive enrichment of the spectrum with harmonics of f 0 and f 1 with no significant participation of orders higher than 2. The appearance of the frequency of mode 9 (f 2 = 16.3 kHz) is also observed during the first transient limit cycle. However, its participation vanishes rapidly before the final limit cycle. Complex modal projection (see Figure 8(d)) confirms the main participations of modes 6, 3, 9 and 11. In the final cycle, mode 9 falls below mode 11 which retains almost constant amplitude. It is also interesting to note that except for mode 6, the amplitudes of most of the modes are not constant in the final limit cycle but oscillate. In particular, the oscillations of the modal projection on mode 3 are very large. It proves that this mode changed notably compared with its linear version. Finally, Figure 8(e) shows the peak-to-mode M-MAC criterion for the final limit cycle at frequencies f 0 and f 1 . This modal identification shows that the peaks at f 0 and f 1 are spatially very close to modes 6 and 3 respectively. Both peaks also have projections on certain stable complex modes.

Case 9 33S

This case represents the case with 9 unstable modes (i.e. µ = 0.35) and an initial condition only on mode 33 with separation as a first occurring non-linear event.

It gives a final state with both adhesion and separation with an amplitude of velocity of the same order for the tangential and the normal velocity. It can be noted that the final state obtained here is obtained for all the initial conditions tested 9 XS except 9 3x33S. This comment will be discussed in more detail in the next section.

Tangential displacements are given in Figure 9(a). The temporal evolution of displacements suggests the existence of three states: the first occurs just after the linear part and finishes at t = 0.004 s. Then the second occurs from t = 0.005 s to t = 0.01 s, and which the final limit cycle is reached at t = 0.015 s. An increase of the tangential displacement is observed when separation occurs. This is due to an increase of the global normal reaction due to the shocks occurring with separation. On the phase plots (see Figure 9(b)), the first small transient cycle (indicated in dashed black line) is visible. Then, the final state corresponds to a quasi stationary cycle. This fact is also observed in the CWT analysis in Figure 9(c): there are two fundamental frequencies (f 0 = 15.7 kHz and f 1 = 27.6 kHz). The first is close to the frequency of mode 9 but significantly lower, the second corresponds to the frequency of mode 24. It can also be seen that the final quasi-stationary response is very complex with the contribution of harmonics and combinations of the fundamental frequencies f 0 and f 1 (2f 0 , 3f 0 , 4f 0 , 2f 0 , 3f 0 -f 1 , f 1 -f 0 , 2f 0 -f 1 , 2f 1 -3f 0 ,...). Then, it can be seen clearly that the fundamental frequency f 1 and combinations of the two fundamental frequencies f 0 and f 1 appear only during the second state (for t > 0.01 s). All these previous comments are also in correlation with the modal projection given in Figure 9(d): the final state is mainly composed of modes 9 and 33, with significant participations of modes 24 and 6. The three different phases can be easily identified with saturation of a specific mode at each state. For example the first state corresponds to the saturation of mode 33. Then mode 9 saturates at a higher energy level, whereas modes 24 and 6 continue to increase. Finally, modes 24 and 6 saturate at a lower level than mode 9. It can be observed that the saturation of modes 24 and 6 induces a reduction of the energy level of modes 9 and 33. Finally, Figure 9(e) gives three M-MAC criteria computed on the peaks at f 0 , f 1 and 2f 0 respectively. It clearly appears that they are spatially similar to mode 9, mode 24 and mode 33 respectively. The main peaks seem to be essentially composed of one unstable mode.

Influence of initial conditions

Finally, the influence of different initial conditions is tested for all cases with 9 unstable modes defined in Table 2. As a reminder, these conditions have been described in Section 6.1.

In both adhesion and separation cases, the limit cycle obtained depends on the initial conditions for the five initial conditions tested and both adhesion and separation cases. Figure 10 represents the evolution of the modal amplitude of the most significant modes for the different initial conditions. It is clear that there are at least two attractors for each functioning condition for the system with 9 unstable modes. One attractor seems to be governed by the most unstable mode (i.e. mode 9, see for example the previous case 9 33S). The other one seems to be closer to the observation made with the 5 unstable modes system where modes 3 and 6 have strong coupled participations, with mode 6 being the main mode for the adhesion cases and mode 3 for the separation cases. To illustrate this and allow a better understanding for the reader, Figures 11(a) and 11(b) give the tangential phase plot and the peak to mode M-MAC at the main peak f 0 for the case 9 3×33S which is the only case of the 9 XS series that gives a different final state. In this case, the final state is mainly composed of frequency f 0 = 10.5 kHz (which is very close to the frequency of mode 3) with the participation of the fundamental frequency f 1 = 13.6 kHz (which is very close to the frequency of mode 6). For the reader comprehension, combinations 2f 0 , f 0 + f 1 and 2f 0 + f 1 are also present inducing the smaller contributions of stable modes that are visible in Figures 11(b).

For the adhesion cases, two attractors are identified. The first is given by case 9 allA while all the other initial conditions tested belong to the second attractor. To facilitate the reader's understanding, it is noted that the convergence to the same final state was verified for modes 3 and 6 in cases 9 3 × 33A and 9 3 × 9A but it requires a much longer time interval in order to be achieved. In both adhesion and separation cases (see Figures 10(a), 10(c) , 10(e) and 10(g) for adhesion cases, and Figures 10(b), 10(d), 10(f) and 10(h) for separation cases), the four initial conditions ending on the same attractor have the same transient behavior as that observed for the cases with 5 unstable modes.

Conclusion on system behavior

In the previous sections, the effect of functioning and initial conditions on the response of a simple friction destabilized system is studied. Dependence on operating conditions is observed once again. More precisely, the system is strongly influenced by the kind of first non linear event (adhesion or separation) that occurs when equilibrium is lost. The effect of the initial conditions is less systematic. Depending on the parameters, different behaviours are identified with two main classes. For the first, obtained for lower friction coefficients when only five modes are unstable, the non-linear system has a final state that depends only on the operating conditions and not on the initial conditions tested here. For the second, In the nine unstable modes cases, different final states for the same operating conditions can be observed.

The use of the complex modal projection allows identifying the contributions of the linear complex modes in the non-linear solution. This technique proves to be useful for monitoring the evolution of the response, especially in the modal competition phase when multiple unstable modes are present. Indeed, in most cases, despite the saturation due to adhesion or separation, the non linear modes do not change greatly in comparison to their linearized contribution obtained for a given operating point via the stability analysis. Moreover, in many cases, peak to mode M-MAC shows that the response at the fundamental frequencies of the final limit cycle are very close to some of the unstable linear modes. However, large oscillations of the modal projections are sometimes found, revealing strong alterations of the nonlinear modes, especially in the transient phase and/or when many modes are present at a similar energy level. CWT analysis of the responses show rather classical results with harmonics and combinations of fundamental frequencies in both transient and stationary phases. Only monoharmonic periodic limit cycles are found in the five unstable mode cases. In other cases, some limit cycles are clearly bi-harmonic. In the latter cases, peak to mode M-MAC proves that non linear modes may significantly change. As only one or two unstable modes remain in the final stationary solution at a high energy level, an interesting point is the selection mechanism of this dominant mode(s). In these numerical simulations, no simple rule could be found. At this stage, the modal competition therefore appears to be a purely transient non linear phenomenon. However, the sparse nature of the modal projections is encouraging regarding the objective of reducing the system. This will be investigated in the following section of the paper.

Implementation of efficient reduction

In this section, model reduction is considered. As the non-linear solving process implies several resolutions of a linear system at each time step, reducing the size of the system is a major concern in order to obtain reasonable computation time let alone perform a computation. One of the most natural methods for reducing the size of the system is computation on a reduced basis B.

It should be noted that reduction of non linear dynamic systems has received much attention in recent years. The concept of non linear modes [START_REF] Kerschen | Nonlinear normal modes, part i: A useful framework for the structural dynamicist[END_REF], which may be defined as a bi-dimensional invariant in the phase space, has been used to perform efficient reduction in certain applications [START_REF] Pesheck | Accurate reduced-order models for a simple rotor blade model using nonlinear normal modes[END_REF][START_REF] Touzé | Hardening/softening behavior in non-linear oscillations of structural systems using non-linear normal modes[END_REF][START_REF] Villa | The invariant manifold approach applied to non-linear dynamics of a rotor-bearing system[END_REF][START_REF] Sinou | Extension of the center manifold approach, using rational fractional approximants, applied to non-linear stability analysis[END_REF][START_REF] Sinou | Methods to reduce non-linear mechanical systems for instability computation[END_REF][START_REF] Sinou | Stability analysis and non-linear behaviour of structural systems using the complex non-linear modal analysis[END_REF]. Other reduction tools have been proposed, like numerical techniques based on invariant inertial manifolds and Proper Orthogonal Decomposition [START_REF] Steindl | Methods for dimension reduction and their application in nonlinear dynamics[END_REF][START_REF] Lenci | Dimension reduction of homoclinic orbits of buckled beams via the non-linear normal modes technique[END_REF]. However, at present, all these methods have required considerable computational effort. The linear reduction methods proposed in the following section need no significant supplementary computations. They are based on the a priori selection of specific linear modes, from the previous analysis of the exact non-linear response of the system. The aim is to evaluate the efficiency of different reduction bases and compare them.

Computation on a reduced basis B consist in searching an approximated solution uapp of the problem in the subspace of R n spanned by basis B. By using the following relations

(ü app , uapp , u app ) = B (ü red , ured , u red ) and ü⋆ app , u⋆ app , u ⋆ app = B (ü ⋆ red , u⋆ red , u ⋆ red ) (26)
The mechanical system

u⋆ T app Mü app + u⋆ T app C uapp + u⋆ T app Ku app = u⋆ T app (f + r c ) (27)
can be rewritten in the following form

M RED üred + C RED ured + K RED u red = f red + r red (28) with
(f red , r red ) = B T (f , r c ) and (M RED , C RED , K RED ) = B T (M, C, K) B (29)
where u red is the coordinate of the approximation in the reduced basis. This reduced equation of the problem may be used instead of Equation ( 14) with the same algorithm for the solving processes. In the following sections, different a priori reduction bases are defined and evaluated on the previous test cases. The quality of the reduced solutions are compared by using the reference computations.

Definition of two kinds of reduction basis

In the definition of the reduction bases, two assumptions are pursued. The first is that since there are strong non-linearities on the contact zone, all the physical contact degrees of freedom have to be kept. The second assumption is that modal reduction is based on the complex modes from the stability analysis due to the fact that these modes have proved capable of giving an accurate representation of the reference solution, and that real modes of the coupled system without friction are sometimes unable to give a correct representation of the complex modes. It is tempting to perform a reduction of the contact interface as well. Without prior knowledge, such reduction is difficult on the solution or on an approximation of it, as would be done with a POD approach [START_REF] Niroomandi | Accounting for large deformations in real-time simulations of soft tissues based on reduced-order models[END_REF] that requires a simulation on the full system for representative small time intervals. To facilitate comprehension, conserving all the physical degrees of freedom on the contact interface also ensures correct representation of contact reaction in the reduction basis. Correct representation of the non-linear force in the reduction basis is a requirement for reduction [START_REF] Boyere | Réduction de modèle en dynamique linéaire et nonlinéaire : Méthode de Ritz, Code Aster Manual[END_REF].

Therefore all the reduced bases tested can be written as:

B = C DOF G DOF = CDOF 0 S GDOF (30)
where C DOF is a basis of the contact physical degrees of freedom (DOF) or of the static constrained modes at the contact interface (i.e. the static response of the structure to a unitary relative displacement on a contact node) and G DOF is a basis of generalized degrees of freedom that are the restriction on the internal degrees of freedom of the complex modeshapes with friction of the system. If C DOF is a basis of static interface modes, part S is a non zero matrix. The inclusion of complex modes in the reduced basis is done by including both the real and the imaginary parts of the complex mode. Co-linear components are then eliminated. The number of basis vectors for an n complex modes basis is therefore at most 2n. Practically, this number is between 1.2n and 1.5n. Two main types of bases are now defined: the first is a classical modal truncation basis, where G DOF is the trace of the complex modes of the complete system below a cut-off frequency -usually taken as 1.5 or 2 times the highest frequency of interest -, but due to strong non-linearities, it is not certain whether conventional truncation rules can be applied. The second is based on the fact that due to non-linearity, the non-linear solution of the self-excited system will be composed by some fundamental frequencies, their harmonics and combinations. The fundamental frequencies are most likely to belong to the frequencies of the unstable modes. That is why G DOF in this type of base is taken to be the trace of the unstable modes plus the trace of the complex modes around the harmonic frequencies of the unstable modes and their combinations.

It is noteworthy that this type of base gives an extension of the space spanned with potentially much higher frequency than usual truncation. However, depending on the number of unstable modes and their distribution, the number of generalized degrees of freedom may not be high enough due to gaps in the frequency band.

Bases of the first type will be referred to as F C n , were n is the number of modes included. Bases of the second type will be denoted IM o , where o is the harmonic of the unstable modes and their order of interaction. Superscript s means that static modes are included in the basis.

Figure 12 shows the modal content retained for different bases used in the following parts of the paper for cases with 1, 5 and 9 unstable modes. The modal sparsity of IM o bases is very clear. One can also see that IM o bases may contain higher frequency content than F C 100 basis (see for examples IM 3 in Figure 12(b), IM 2 and IM 3 in Figure 12(c)).

Performance of reduction bases 7.2.1. Presentation of case studies

In the following sections, it has been chosen to focus only on the six cases 1 9A, 1 9S, 9 allA, 9 33S, 5 9A and 5 allS to validate the robustness of the different bases. All the studied cases and the associated conditions are given in Table 2. For the reader comprehension, these tests give an overview of the effectiveness or otherwise of all bases in configurations of varying complexity with one, 5 and 9 unstable modes (for µ = 0.075, µ = 0.19 µ = 0.35), separation or adhesion and different initial conditions. Indeed, the following case studies will show that both types of bases are able to perform efficient and accurate reduction but that some bases are limited by the robustness of their results, meaning that the accuracy of the reduced solution may depend on the initial conditions chosen and not only on the reduced system itself.

For the reader comprehension, the following cases are not processed in order of increasing the number of unstable modes, but in order of increasing difficulty and complexity to have satisfactory results in terms of reduction bases.

Finally, quantitative indicators were defined for testing and comparing the various bases. They were defined on the FFT of the stabilized non-linear stationary solution. If not stated otherwise, all the FFT presented are kinetic power spectra computed on the last 3 ms of the solution given by:

E c (ω) = 1 2 FFT( u, ω) ⋆T M FFT( u, ω) (31)
For each reduced computation, the main peaks of the FFT are paired with the reference ones by using two criteria: frequency proximity and peak-to-peak M-MAC. Once pairing is done, frequency shift, δf induced by reduction, amplitude error δa and peak-to-peak M-MAC are used to quantify error. The peak-to-peak M-MAC are given with a precision of 10 -2 . The first main peaks of the FFT are retained if the secondary ones are no less than 20 dB below the first. When peak-to-peak M-MAC between the reference and reduced solutions is below 0.6, pairing is considered to have failed and other indicators (frequency and amplitude error) on the unpaired peaks are not given.

Reduction for cases 1 9A and 1 9S

First of all, these two cases represent the case with one unstable mode and an initial condition only on mode 9 with adhesion or separation as a first occurring non-linear event. The reduction results on the simplest cases are very encouraging and show that it is quite possible to obtain accurate results with a good level of reduction. Table 3 gives the three main peak amplitudes and the associated frequencies of the FFT of the reference solution up to the frequency resolution of the FFT. The error indicator given in Table 4 shows that agreement between the reference solution and all the reduced bases is very good on the limit cycle for both cases 1 9A and 1 9S. For example, Figure 13(a) shows that there is very good agreement between the reduced and the reference phase plots on the contact zone for the first case 1 9A. Evolutions of the modal projection of mode 9 (see Figure 13(b)) also illustrate that all the reduced bases are efficient. The peak-to-peak M-MAC between the reference peak and the reduced peaks always has a modulus of 1. These quantitative indicators increase confidence in the ability of this method to perform model reduction for friction destabilized systems. Then, considering more specifically the case with separation (i.e. case 1 9S), the error increases slightly for all bases without static modes, and especially for the F C 100 basis as shown in Table 4. It may be noted that the subspaces spanned by all the IM bases of the system with only one unstable mode are strictly included in the subspace spanned by F C 100 and that the error of the latter is greater. This fact highlights the difficulty of performing accurate non-linear model reduction since adding modal vectors to the reduction basis may decrease the quality of the solution computed.

Even though very efficient reduction can be achieved for models with only one unstable mode (i.e. for the reader comprehension, it can be seen that the basis IM 1 only has two generalized degrees-of-freedom corresponding to the real and imaginary parts of the unstable mode and thus leads to an impressive reduction), cases with separation exhibit a considerable difference of behavior between the bases that have static modes and the others. Error indicators for the F C s 100 basis remain at the same level for adhesion and separation cases, whereas the performance of the F C 100 basis can be seen to decrease. It can be seen that including static residue in the basis seems to be a very efficient way of ensuring quality for both types of non-linear event. It is also interesting to see that the IM bases perform better than F C 100 basis for separation, though slightly worse for adhesion. Figure 13(c) and 13(d) show that due to separation, reduction leads to visible error on the phase plots. This fact is more particularly visible for an internal node, as indicated in Figure 13(d).

Reduction for case 9 allA

Studying the case with 9 unstable modes with operating conditions leading to adhesion gives interesting results. The reference main peak amplitudes and frequencies are given in Table 3. This system appears easy to reduce since almost all the bases give very good results for the indicators shown in Table 5. Basis IM 1 leads to a wrong cycle that corresponds to the other attractor observed in Section 6.

Figure 14(a) shows very good agreement for all the bases except IM 1 for the power density spectrum. The phase plot in Figure 14(b) shows some defects in the reduction performed with IM 2 and F C 100 . However IM s 2 and F C s 100 agree well with the reference. It should be noted that comparison is not easy for quasi-periodic cycles. The choice is made here to focus on cycle size, the number of loops and global shape. Finally, Figures 14(c), 14(d), 14(e) and 14(f) illustrate the modal projection for modes 3, 6, 9 and 33 respectively. It confirms that most bases agree very well with the reference solution. The behavior of the IM 2 basis on the transient phase between 5 ms and 30 ms shows excessive contributions of modes 9 and 33 but the system returns to the correct solution. Moreover, as a consequence of the reduction, the transient evolution of non-linear behavior which is composed of all the modal contributions can be time-shifted in comparison to the reference solution. In this present case, the different phases are simply in advance, but their respective amplitudes are correct. In conclusion, although the proposed reduced bases give very good results for the final stationary state, there may be a speed-up phenomenon in the transient behavior of the non-linear system.

Reduction for case 9 33S

Here we present a case with 9 unstable modes and an initial condition only on mode 9 with seperation as a first occurring non-linear event. Table 3 gives the main peak amplitudes and frequencies for the reference computation. This case is interesting due to the significant participation of combinations of the two fundamental Case Peak number Frequency (kHz) Amplitude (J) (dB) 5) clearly show the efficiency of the IM bases in increasing order: a significant decrease of the frequency error IM 1 to IM 3 is observed. Amplitude and MAC errors also show the same trend. Error indicators also outline the efficiency of static modes for increasing reduction accuracy as bases IM s 2 and F C s 100 are both very close to the reference computation.

The power spectrum density of the final state for each basis is shown in Figure 15(a). All the bases are capable of reaching the main fundamental frequency peak at f 0 = 15.67 kHz and the harmonics of order 2 and 3. The only basis able to reproduce the peak at f 1 is F C s 100 . Therefore, it is also the only basis showing peaks with combinations of f 0 and f 1 . Basis IM s 2 also exhibits very small secondary peaks corresponding to more visible ones in the reference.

The phase plot comparison in Figure 15(b) clearly shows that the harmonic enrichment of the IM bases makes it possible to obtain a phase plot very similar to the reference even though the amplitude is 35% too low. Basis IM s 2 is the closest to the reference phase plot, with an accurate amplitude and shape. F C s 100 has a very correct shape but a higher amplitude than in the reference.

Finally, Figures 15(c),15(d), 15(e) and 15(f) show the temporal evolution of the modal projection on modes 3, 6, 9 and 33 for the different bases. For modes 3 and 6 (see Figures 15(c) and 15(d)), the only bases that obtain the correct final amplitude for the mode are the two bases with static modes F C s 100 and IM s 2 . However, their behavior around t = 5 ms during the transient phases is not in good agreement with the reference with a significant overestimation of the growth rate of the mode. The amplitude of this mode with bases F C 100 and IM 3 at t = 5.5 ms is closer to the reference, but the mode is stabilized at this amplitude for the rest of the time interval. IM 2 exhibits the same stabilization as IM 3 and F C 100 , but with greater amplitude. With basis IM 1 , the amplitude decreases after t = 5.5 ms until t = 10 ms, where the amplitude is stabilized. For mode 9 (see Figure 15(e)), the mode with the greatest contribution, all the bases give a good amplitude with a perfect match for bases IM 3 , F C s 100 and IM s 2 and with the same small shift for the three other bases. The shift for the other bases occurs very early when the amplitude of mode 9 is still negligible and may result from numerical noise. One of the most noticeable differences between the two sets of bases is that the F C 100 , IM 1 and IM 2 bases have a smoother transition toward the stabilization of the amplitude around 5 ms. For mode 33 (see Figure 15(f)), all the bases give the same result until t = 4 ms. Then IM 1 has a strong oscillation and the amplitude of mode 33 is stabilized at 0.14. Bases IM 2 and F C 100 have an abrupt increase of amplitude, but not as much as in the reference, and then they stabilize at t = 0.28s. Basis IM 3 follows the reference solution exactly until t = 11 ms and fails to reproduce the slight drop in amplitude visible in the reference. All the amplitudes of mode 33 with bases F C 100 , IM 1 , IM 2 and IM 3 are stabilized without any oscillation. Bases F C s 100 and IM s 2 are very close to the reference, even though they give a lower amplitude than the reference after t = 5 ms. They are both able to reproduce the drop that occurs in the reference at t = 11 ms, but they anticipate it at t = 8.5 ms.

Reduction or case 5 9A

This case represents the case with 5 unstable modes and an initial condition only on mode 9 with adhesion as a first occurring non-linear event. It shows that choosing a "correct" basis may be difficult for accurate reduction. To facilitate the reader's comprehension, Figures 16(a) and 16(b) illustrate the FFT of the final non-linear state and the corresponding phase plot. Table 3 gives the values (i.e. frequency and amplitude) of the three main peaks for the reference non-linear stationary solution. Table 5 shows the efficiency of the different bases tested. It can first be seen that the basis IM 3 fails to obtain better results than IM 2 . It is even less accurate than IM 1 . Both bases IM 1 and IM 3 tend to stabilize the driven phase of mode 9. It is also clear that inclusion of static modes helps obtaining the correct frequency for the two main peaks: IM s 2 basis is perfectly correlated with the reference on the second peak frequency when IM 2 exhibits a frequency shift of about 1%. Figure 16(a) shows that static modes significantly change the amplitude between the peaks of the FFT. However, they do not necessarily improve its accuracy since IM s 2 has a much greater error than IM 2 . The phase plots on Figure 16(b) show that bases F C 100 , IM 2 and F C s 100 agree reasonably with the reference on the contact zone. IM s 2 has greater amplitude on the contact zone than the reference and the two last bases give a wrong shape corresponding to mode 9 instead of mode 6.

Comparison of the modal amplitudes for modes 3, 6 and 9 on Figures 16(c), 16(d) and 16(e) shows that bases IM 1 and IM 3 have stabilized the cycle composed of mode 9 alone, with a steady amplitude for modes 3 and 6 at very low energy levels. Another striking phenomenon is the fact that reduction may lead to considerably accelerating the evolution of the solution: saturations of modes 3 and 6 occur before 0.06 s for bases F C 100 , F C s 100 , and IM s 2 and around 0.075 s for IM 2 and only after 0.1 s for the reference solution. This speed up relies on two effects: first, the amplitudes of these modes increase strongly when mode 9 reaches saturation. This increase is overestimated for all the reduced bases showing this speed up. Also, the growth rates of the modes are also very often wrongly estimated: they are a little too high except for IM 2 for which the growth rate is a little too low. These wrong estimations can also explain the failure of bases IM 1 and IM the increase of modes 3 and 6 when mode 9 saturates and the growth rates are also underestimated at almost zero.

] [ 1 1 1 ] IM s 2 [ 0 0 - 0 ] [ 0.73 2.17 - 3.62 ] [ 1 1 0.04 ] 5 9A F C 100 [ 0 0 ] [ 0.13 1.19 ] [ 1 1 ] IM 1 [ - - ] [ - - ] [ 0 0.2 ] IM 2 [ 0 1.16 ] [ -0.69 0.45 ] [ 1 0.73 ] IM 3 [ - - ] [ - - ] [ 0.01 0.2 ] F C s 100 [ 0 0 ] [ -0.01 0.26 ] [ 1 1 ] IM s 2 [ 0 0 ] [ 2.1 2.81 ] [ 1 0.66 ] 5 allS F C 100 [ - - ] [ - - ] [ 0 0 ] IM 1 [ 0 - ] [ -19.33 - ] [ 0.98 0.05 ] IM 2 [ - - ] [ - - ] [ 0 0 ] IM 3 [ - - ] [ - - ] [ 0.02 0.02 ] F C s 100 [ 0 0 ] [ -0.09 0.51 ] [ 1 0.97 ] IM s 2 [ - - ] [ - - ] [ 0 0 ]
Finally, Figures 17(a) and 17(b) show the temporal evolutions of tangential velocity and displacement on the contact zone for the reference computation and bases IM 2 and F C s 100 . When comparing the displacements (Figure 17(b)), it can be seen that the transient behavior is well represented by both bases. Especially the first phase, with the growth of mode 9 and the surge at t = 0.5 ms. The second phase is anticipated on both reduced bases, but besides the anticipation it seems very close to the reference computation. The comparison of the velocities (Figure 17(a)) again shows very good agreement in the first phase for both bases and also a very similar transition between phases, even though transition is anticipated in the reduced computations. However, in the second phase, basis IM 2 clearly exhibits an amplitude of greater velocity than the reference.

Reduction for case 5 allS

This case represents the case with 5 unstable modes and an initial condition on all unstable modes with separation as a first occurring non-linear event. This case appears interesting due to the fact that reduction appears quite difficult: the only basis giving really satisfactory results is F C s 100 . The characterization of the three main peaks is given in Table 3. Table 5 also shows that IM 1 is the second basis in terms of representation of the final state with good accuracy in the frequency and shape of the first peak but with a significantly wrong amplitude. All the other bases fail to have a peak-to-peak M-MAC criterion greater than 0.6 and most of them simply have an value near zero for the M-MAC criterion.

Comparison of the FFT on Figure 18(a) shows that all the bases except F C a predominance of mode 9. IM 1 has a bi-harmonic cycle with modes 3 and 9. Phase plots on Figure 18(b) shows that the reduction has induced in most bases a very wrong estimation of behavior on the contact zone with exaggerated amplitudes.

Modal projection on modes 3, 6 and 9 (see Figures 18(c), 18(d) and 18(e)) confirms that most of the bases produce a cycle composed of mode 9. Modes 3 and 6 tend to have fluctuating amplitudes, meaning that the modes of the system represented by the reduced model are much more affected by non-linearity than the reference system.

Conclusions of the performance of reduction bases

This review of the performance of reduction bases on the different cases brings to light the following results: both the modal truncation and the harmonic bases are able to give an accurate result on the four first cases. It is clear that static modes are a very efficient way to improve reduction accuracy in the presence of separation. However, even if both types of basis are able to give good results, the latter two cases also show limitations of the harmonic bases and illustrate the fact that when reducing a non-linear system, increasing the reduction basis does not necessarily improve the accuracy of the basis, e.g. in case 5 allS, IM 1 is the best of the IM bases, despite the fact that it is always present in each of them.

Finally, the efficiency of the proposed reductions in term of computational time is now discussed. For the different cases presented in the previous section, Table 6 gives the ratio between the computation time of the reduced system and the reference solution (by performing the same initial conditions in all cases). Unsurprisingly, the smallest basis as the smallest ratio, and IM bases tend to have increasing speedup ratio with increasing of the number of unstable modes. Moreover, It can be observed that separation tends to slightly increase the ratio for the bases without static modes. For θ = 0.5, this integration scheme induces no numerical damping for linear problems and is unconditionally stable. However, when dealing with shocks between elastic media, any other value than θ = 1 is unable to guarantee the inelastic shock predicted in the continuous framework. A workaround to this issue is to evaluate the relative displacement for the computation of the contact reactions with another integration scheme given by Equation (A.6). This methods is called the modified θ-method.

ũi+1 = u i+1 + τ (1 -θ) ui+1 (A.6)
Using this modified estimation of displacement introduces a contact reaction before the shock occurs. This contact reaction is such that the shock occurs exactly on a time step. Imposition of the inelastic shock with the modified θ-method introduces numerical damping in the integration scheme which depends on the spatial discretization of the contact interface. The algorithm describing the solver used in Section 3.2 is given by Algorithm 1.

Algorithm 1 Fixed point algorithm for nonlinear contact dynamics Require: J, f 1..N , crit, iter max 1: for i = 1 to N do f mtd ← compute sndm(f i,i+1 , u i , ui ) ( un , ut , u n , u t ) ← local projection k ui+1 , k u i+1 10:

(r n , r t ) ← local projection k r i+1 11:

r n ← proj R-r n -ρ n ( k un + 1 τ b) 12:
r t ← proj Kµ r t + nr n -ρ t k utv g 13:

k+1 r i+1 ← global projection(r n , r t ) end while 18: end for

Figure 1 :

 1 Figure 1: System under study (a) Contact problem on an elastic body (b) Simple test case: 2D elastic layer

  Sc u⋆ • rdS ( u, r) verifying Equation (1 and (u, r) verifying Equation (2) on S c U = {u, u = u d on ∂Ω u } and U 0 = {u, u = 0 on ∂Ω u } ε and ε ⋆ are the symmetric gradients of u and u ⋆

  0 on S c r e n ≤ 0 and r e t = -µr e n v g v g = -µr e n t on S c u e n r e n = 0 on S c ε e and ε ⋆ are the symmetric gradients of u e and u ⋆ where (n, t, b) is a local basis on the contact zone, with n being the outward normal, t is the sliding direction: t = vg vg and b = n ∧ t.

  b -rn t on S e c ε and ε ⋆ are the symmetric gradients of û and u ⋆

Figure 2 :

 2 Figure 2: Divergence rate with respect to frequency (× : unstable modes, o: stable modes) -(a) µ = 0, (b) µ = 0.19, (c) µ = 0.35, (d) µ = 0.5

Figure 3 :

 3 Figure 3: Mode shapes of the unstable modes with µ = 0.5 (mode number, frequency, divergence rate and LT (or RT) for left (or right) travelling waves) -(a) mode 3, 10.52 kHz, 6, 5% LT, (b) mode 6, 14.25 kHz, 7.7% LT, (c) mode 9, 16.36 kHz; 11.6% LT, (d) mode 11, 18.5 kHz, 8.1% RT, (e) mode 13, 20.83 kHz, 1.0% RT, (f) mode 15, 21.45 kHz, 3.4% LT, (g) mode 20, 24.15 kHz; 5.4% LT, (h) mode 24, 27.25 kHz, 3.3% LT, (i) mode 31, 31.20 kHz, 2.9% RT, (j) mode 33, 32.16 kHz, 1.8% LT, (k) mode 42, 34.01 kHz, 0.8% LT, (l) mode 50, 38.29 kHz, 0.4% LT

6. 3 .

 3 Cases with 9 unstable modes 6.3.1. Case 9 9A

Figure 4 :Figure 5 : 2 0Figure 6 :Figure 7 :Figure 8 :Figure 9 :Figure 10 :Figure 11 :

 45267891011 Figure 4: Case 5 9A -(a) Tangential velocity on a contact node (b) Tangential displacement on a contact node (gray line: mean value) (c) Phase plot on contact zone (gray: complete history, black: final cycle, dashed black: first cycle) (d) CWT of tangential velocity on contact zone, (e) Complex modal projection where "i" defines the unstable mode number for "Mi", (f) Peak to mode M-MAC at f 0

32

 32

Figure 12 :Figure 13 :

 1213 Figure 12: Modal content of different bases and stability analysis. (red: unstable modes, blue: stable modes, black: selection of modes) -(a) 1 unstable mode, (b) 5 unstable modes, (c) 9 unstable modes

Figure 14 :Figure 15 : 2 -Figure 16 :Figure 17 :Figure 18 :

 14152161718 Figure 14: Bases comparison for case 9 allA, +: reference, o: F C 100 , * : IM 1 , •: IM 2 , ×: F C s 100 , : IM s 2 -(a) FFT of the final state, (b) Phase plot, (c) Modal projection of mode 3, (d) Modal projection of mode 6,(e) Modal projection of mode 9, (f) Modal projection of mode 33

4 : k r i+1 ← r i 5 : 7 :ε ← crit + 1 8 :

 4578 k ui+1 ← solve J u = f mtd + k r i+1 6: b ← compute modified gap(u i ) while ε ≥ crit and iter max ≥ k do 9:

14 :

 14 k+1 ui+1 ← solve J u = f mtd + k+1 r i+1

  Thus comparison between |β i | and |β j | does indeed provide information on which mode has greater participation, but no conclusion can be obtained by comparing |β i | and |β j | + |β k |. The solution around the equilibrium point is given as follows:

Table 2 :

 2

		Number of	µ	Selected modes for	Initial	Non-linear Reduction	Initial
		unstable modes		initial condition	nonlinear event	behavior		condition
	1 9A	1	0.075	9	adhesion		X
	1 9S	1	0.075	9	separation		X
	5 9A	5	0.19	9	adhesion	X	X	X
	5 9S	5	0.19	9	separation	X		X
	5 20A	5	0.19	20	adhesion			X
	5 20S	5	0.19	20	separation			X
	5 3x9A	5	0.19	3 and 9	adhesion			X
	5 3x9S	5	0.19	3 and 9	separation			X
	5 3x20A	5	0.19	3 and 20	adhesion			X
	5 3x20S	5	0.19	3 and 20	separation			X
	5 allA	5	0.19	all	adhesion			X
	5 allS	5	0.19	all	separation		X	X
	9 9A	9	0.35	9	adhesion	X		X
	9 9S	9	0.35	9	separation			X
	9 33A	9	0.35	33	adhesion			X
	9 33S	9	0.35	33	separation	X	X	X
	9 3x9A	9	0.35	3 and 9	adhesion			X
	9 3x9S	9	0.35	3 and 9	separation			X
	9 3x33A	9	0.35	3 and 33	adhesion			X
	9 3x33S	9	0.35	3 and 33	separation			X
	9 allA	9	0.35	all	adhesion	X	X	X
	9 allS	9	0.35	all	separation			X

Description of the different cases. The symbol X denotes the part of the paper where the case is studied: " Nonlinear behavior" for the study of the transient analysis; "Reduction" for the performance of reduction bases; "Initial condition" or the study of the influence of initial conditions

6.2. Cases with 5 unstable modes 6.2.1. Case 5 9A

This case represents the case with 5 unstable modes (i.e. µ = 0.19) and an initial condition only on mode 9 with adhesion as a first occurring non-linear event.

Table 3 :

 3 Kinetic Power Spectrum Peaks -reference solutions for each case frequencies f 0 = 15.67 kHz and f 1 = 27.67 kHz in the final state. Error indicators on the two main peaks (see Table

Table 4 :

 4 Error indicators on first peak for the different bases -Cases 1 9A and 1 9S

	Case base	δf δa( dB)
	1 9A F C 100 0	0.38
	IM 1	0	0.28
	IM 2	0	0.28
	IM 3 F C s 100	0 0	0.28 0.004
	1 9S F C 100 0	0.032
	IM 1	0	0.053
	IM 2	0	0.053
	IM 3 F C s 100	0 0	0.053 0.005

  3 . They underestimate

	Case	basis	δf (%)	δa(dB)	P2P M-MAC
	9 allA F C 100 [ 0 IM 1 [ -IM 2 [ 0 F C s 100 [ 0 IM s 2 [ 0 9 33S F C 100 [ 2.13 2.13 -2.13 ] [ 7.14 0 0 0 ] [ -0.05 ---] [ -0 0 0 ] [ -0.06 0 0 0 ] [ 0.01 0 0 0 ] [ -0.02 IM 1 [ 4.26 4.26 -4.26 ] [ 7.4 IM 2 [ 2.13 2.13 -2.13 ] [ 7.01 IM 3 [ 0 0 -0.71 ] [ 0.28 F C s 100 [ 0 0 0 0 ] [ 0.08	0.09 -0.15 -0.08 -0.15 -2.46 -7.95 -2.92 -0.02 -0.68 ] [ 1 -] [ 0 -0.05 -0.29 ] [ 1 0.34 0.32 ] [ 1 0.25 0.33 ] [ 1 2.58 -4.57 ] [ 0.98 0.98 0.04 0.97 ] 1 1 ] 0.22 0 ] 1 1 ] 1 1 ] 1 1 ] ] [ 0.98 0.99 0.04 0.14 ] -27 -3.4 1.63 -6.34 ] [ 0.98 0.98 0.04 0.93 ] 0.25 5.28 1.71 ] [ 1 1 0.04 ] 0.24 0.43 0.89

Table 5 :

 5 Error indicators (Case 9 allA at f 0 = 14.33 kHz, 2f 0 = 28.67 kHz, f 1 = 10 kHz, and f 0 -f 1 = 4.33 kHz; Case 9 33S at f 0 = 15.67 kHz, 2f 0 = 31.33 kHz, f 1 = 27.67 kHz and 3f 0 = 47 kHz; Case 5 9A at f 0 = 14.33 kHz and 2f 0 = 28.66 kHz; Case 5 allS at f 0 = 10.33 kHz and 3f 0 = 31.33 kHz)

Table 6 :

 6 Computational time reduction

-(a) FFT of the final state, (b) Phase plot, (c) Modal projection of mode

3, (d) Modal projection of mode 6,(e) Modal projection of mode 9

Conclusion

A friction destabilized system was studied in this paper in which the destabilization mechanism corresponds to the mode coupling phenomenon due to friction force with a Signorini-Coulomb contact law. The non-regularized problem with finite element discretization was studied by including a stability analysis and the transient and stationary dynamics of the non-linear self-excited system. Parametric studies were conducted with operating conditions and initial conditions as parameters. They showed that for systems with low numbers of unstable modes, the limit cycle depends only on the operating parameters. For systems with more unstable modes, for each functioning condition, two limit cycles were found depending on the initial conditions.

Different reduction bases were then proposed and tested. They were based on the complex modes resulting from the stability analysis. The comparison of the reduced model computation with the reference model showed that accurate model reduction for friction destabilized system is possible but requires great care. The inclusion of the unstable modes and modes near the harmonics and combinations of the unstable mode frequencies in the basis is not sufficient to ensure an accurate prediction independent of the initial condition considered. All the stable modes in the frequency range must be included to ensure an accurate prediction of the non-linear response. If separation at the interface is present, the reduction basis must contain a correct representation of the static residue.

Finally, even with the most complete bases, it was observed that model reduction may also lead to an alteration of the transient evolution of modal participations in the temporal solution, often leading to an advance phenomenon.

Appendix A. Integration scheme and solving algorithm

The integration scheme used for the dynamics solver is based on the classical θ-method integration scheme with θ = 0.5. In this section integration rules used in this scheme are recalled and description of the modified θ-method for dealing with shocks is given.

Integrating Equation ( 14) over the whole time interval [0; T ] leads to Equation (A.1).

with (u, u, r c ) verifying Equations ( 4) and ( 5) at each mesh node. Using a regularly spaced time discretization with a time step τ , and using Equation (A.2) for approximating the integration of contact reaction and Equation (A.3) integration of other quantities between time steps t i and t i+1 , time discretized equations governing the system is given by Equation (A.4).

with (u, u, r c ) verifying Equations ( 4) and ( 5) at each mesh node. Thus the discrete time integration of the system is obtained using Equation (A.5)


A Nguyen

L Bergé

P González De Alaiza Martínez

I Thiele

S Skupin

THz Emissions from Air-Plasmas Created by Midand Far-Infrared Two-Color Femtosecond Pulses

We study THz emission by two-color femtosecond filaments in air using mid-to far-infrared pump wavelengths. 3D numerical simulations show that 10.6-µm laser pulses can produce THz fields with unequaled mJ energies and GV/m amplitudes.

Introduction

Terahertz (THz) science is attracting a broad interest because many molecular fingerprints belong to the THz range [START_REF] Tonouchi | Cutting-edge terahertz technology[END_REF]. A well-known technique to create efficient THz emitters is to focus a two-color femtosecond light pulse into air, creating a plasma that acts as frequency converter [START_REF] Kim | Coherent control of terahertz supercontinuum generation in ultrafast laser-gas interactions[END_REF]. Recent studies [START_REF] Clerici | Wavelength scaling of terathertz generation by gas ionization[END_REF][START_REF] Nguyen | Spectral dynamics of THz pulses generated by two-color laser filaments in air: the role of Kerr nonlinearities and pump wavelength[END_REF] showed that increasing the pump wavelength can enhance the laser-to-THz conversion efficiency which is usually limited to about 10 -4 for a 800 nm pump pulse. Moreover, it is known that the pump power contained in a filament is limited to a few times the critical power for selffocusing P cr before undesired spatial breakup and multiple filamentation kicks in [START_REF] Campillo | Periodic breakup of optical beams due to selffocusing[END_REF]. Because P cr scales quadratically with the pump wavelength, increasing the latter makes filaments with much higher energies possible. Longer laser wavelengths, e.g., 3.9 µm, are nowadays accessible in the femtosecond regime and energetic, 100-fs CO 2 laser pulses operating at 10.6 µm should be available in the coming years [START_REF] Pogorelsky | BESTIA -The next generation ultra-fast CO 2 laser for advanced accelerator research[END_REF]. We here focus on intense, broadband THz fields produced by such long-wavelength pulses. By means of a three-dimensional (3D) unidirectional solver [START_REF] Nguyen | Spectral dynamics of THz pulses generated by two-color laser filaments in air: the role of Kerr nonlinearities and pump wavelength[END_REF], we study THz generation by two-color femtosecond pulses ionizing air. Our numerical results evidence unprecedented THz energy yields and conversion efficiencies achieved by 10.6-µm pump pulses along extended propagation distances.

THz spectra and supercontinuum generation of mid-and far-infrared pulses

In our simulations, the input light field has a Gaussian profile in space and time with 2.5-mm 1/e 2 transverse radius and 100-fs FWHM duration. It is composed of the fundamental (FH) pulse operating at the carrier wavelength λ 0 and 10% of second harmonic (SH) in energy. The overall beam is linearly focused by a converging lens of 30-cm focal length. All pulses have an input power equal to 1.73P cr of the FH pulse and the dispersion map includes HITRAN data for CO 2 and water absorption [START_REF]Hitran on the web[END_REF]. Peak intensities of ∼ 200 TW/cm 2 produce electron densities above 10 17 cm -3 . When increasing the pump wavelength and pump energy (P cr ∝ λ 2 0 ), longer filaments are promoted and they start to self-focus at shorter propagation distances. Figure 1(a) illustrates the resulting THz energy yields below 10 THz that reach the 0.1 mJ level for 3.9-µm pumps and several mJ for 10.6-µm lasers. Figures 1(b,c,d) detail the THz spectra computed at the distance of maximum THz field production, i.e., where the on-axis THz fields attain their maxima. Because an important question has been the role of Kerr-induced four-wave mixing compared to that of photocurrents [START_REF] Andreeva | Ultrabroad terahertz spectrum generation from an air-based filament plasma[END_REF], Figs. 1(b,c,d) also compare the spectra computed with or without the plasma terms before the linear focus. The presence of plasma clearly increases the THz spectral intensity by at least three orders of magnitude and shifts the frequency centroid down to 0.2 THz for the three pump wavelengths. Hence, plasma generation -even when it involves low electron densities ∼ 10 15-17 cm -3 -is here the key player in THz pulse generation.

Figures 1(e,g,i) display the entire optical spectrum developed for the three FH wavelengths along the propagation axis. When increasing the pump wavelength, high-order harmonics are produced more efficiently and form a broader supercontinuum through the plasma nonlinearity and by self-and cross-phase modulation. They are able to contribute to photoionization. For a near-IR pump (0.8 µm), the distance of maximum THz generation coincides with that of maximum plasma generation close to focus. For larger wavelengths, maximum THz generation shifts towards shorter propagation distances. For a pump wavelength of 10.6 µm, THz production is maximum at z ≈ 20 cm only, well before maximum plasma generation sets in, which is mainly due to high-order harmonics that do no longer contribute to THz generation. Figures 1(f,h,j) reveal that on-axis THz fields induced by a 10.6-µm pump are two orders of magnitude higher than with a near-IR pump pulse, and twice as high as with a 3.9-µm pump pulse. Moreover, the transverse THz field profile broadens with the pump wavelength (not shown). The THz conversion efficiency obtained at 10.6 µm is thus three times the one at 3.9 µm, that is 3% versus 1%, respectively. In addition, we report much smaller angles of the conical THz emission and a better coupling between the two colors along extended propagation ranges, which renders filamentation setups using CO 2 lasers particularly promising for further applications in THz spectroscopy.

Conclusion

In summary, we numerically showed that two-color ultrashort CO 2 laser filaments can provide THz conversion efficiencies up to 3%. Because far-infrared femtosecond filaments contain hundreds of mJ and are able to self-channel over long distances along which THz production is stable, they can supply THz energy yields of several mJ. Such energy levels are record values for 100-mJ pump pulses operating at moderately high (non-relativistic) laser intensities.

Fig. 1 .

 1 Fig. 1. (a) THz energies numerically computed for two-color Gaussian pulses with FH wavelengths λ 0 = 0.8, 3.9 and 10.6 µm. Vertical dotted lines indicate the distances of maximum THz production. (b,c,d) THz spectra at these distances, namely, (b) z = 30 cm, (c) z = 27 cm and (d) z = 19.5 cm. Gray curves correspond to the same plotted quantities without plasma. Note the absorption line ∼ 6 THz due to water vapor. (e,g,i) Intensity spectra along the propagation axis. (f,h,j) On-axis THz fields obtained by filtering (20th-order Butterworth) the total field in the frequency range ν < 10 THz.