Van-Vuong Lai

Olivier Chiello
email: olivier.chiello@ifsttar.fr

Jean-Franc ¸ois Brunel

Philippe Dufrénoy

Full finite element models and reduction strategies for the simulation of friction-induced vibrations of rolling contact systems

Keywords: Curve squeal, Friction-induced vibrations, Rolling contact, Finite element method, Structural instability, Mode coupling, Railway noise

The aim of this paper is to develop a full Finite Element (FE) computational method for the dynamics of unstable frictional rolling contact systems in order to calculate reference solutions in many applications, especially curve squeal for railway transportation but also roller bearing or metal cutting. The proposed method is characterized by the use of a fine FE discretization of the contact surface in an Eulerian frame, nonlinear frictional contact laws and model reduction techniques. An application to the frictional rolling contact between two annular cylinders is presented in both quasi-static and dynamic cases with mode coupling instabilities. The validation of the approach in quasi-static conditions is carried out by comparison with CONTACT software. Stability and transient results

show that the technique is able to simulate friction-induced vibrations at high frequencies. Reduced models are tested and show a good agreement with the full model.

Introduction

A large variety of mechanical system (road vehicle tyres, railway wheels and roller ball bearings, gear box, belt/pulley, lead-screw drives, metal cutting, etc.) generates rolling contact forces which are transferred within a finite area of contact between the rolling element and the substrate. The main reason to make use of rolling contacts in different applications is the low resistance to motion. Even if friction mechanisms are rather weak, they may be an important source of vibration and noise.

As an example in the railway field, the squeal noise emitting by rail-bound vehicles in tight curves (radius lower than 200m) is characterized by high sound pressure levels (130 dB at 0.9m from the wheel) at pure medium and high frequencies. In urban areas where tight curves are numerous, squeal may affect many passengers and local residents and it is necessary to reduce this noise. The modeling of curve squeal can contribute towards an understanding of the generation mechanism and the effects of the different kinematic and mechanical parameters.

State of the art brings many models able to simulate curve squeal, distinguished according to mechanisms leading to squeal, wheel/rail contact models and solution domains (time or frequency).

Although longitudinal wheel slip and wheel flange contact have originally been cited as a cause of curve squeal, they have been discredited as a main energy input in several works [START_REF] Rudd | Wheel/rail noise-part ii: Wheel squeal[END_REF][START_REF] Remington | Wheel/rail squeal and impact noise: What do we know? what don't we know? where do we go from here?[END_REF][START_REF] Vincent | Curve squeal of urban rolling stockpart 1: State of the art and field measurements[END_REF]. Most of mechanisms proposed in literature put forward the high lateral slip of the wheel on the rail-head as the main cause of curve squeal. Indeed, in tight curves, a steady lateral sliding motion is imposed to the wheel due to the misalignment between wheel and rail (angle of attack or yaw angle). A widespread assumption is that the friction forces generated by the sliding motion may lead to structural instability and self-sustained vibration of the wheel/rail system. Two instability mechanisms have been proposed: falling friction and mode coupling. The falling friction was first mentioned by Rudd [START_REF] Rudd | Wheel/rail noise-part ii: Wheel squeal[END_REF]. A falling slope in the velocity-dependent friction law can be mathematically expressed as negative damping leading to an unstable behavior. A large number of squeal models are based on this mechanism, as in the former paper of Rudd [START_REF] Rudd | Wheel/rail noise-part ii: Wheel squeal[END_REF] or the more recent works of Heckl et al. [START_REF] Heckl | Curve squeal of train wheels, part 1: mathematical model for its generation[END_REF][START_REF] Heckl | Curve squeal of train wheels, part 2: Which wheel modes are prone to squeal?[END_REF] which explain in detail the mechanism by using frequency-domain and time-domain models. The consideration of mode coupling instability in curve squeal models is more recent. Although such a mechanism may occur with a constant friction coefficient, it needs to take into account the vertical dynamics of the wheel/rail system in order to work. Mathematically, this is the nonsymmetric stiffness matrix due to friction which generates the instability. Models proposed in [START_REF] Glocker | Curve squealing of trains: Measurement, modelling and simulation[END_REF][START_REF] Pieringer | A numerical investigation of curve squeal in the case of constant wheel/rail friction[END_REF][START_REF] Zenzerovic | An engineering time-domain model for curve squeal: Tangential point-contact model and green's functions approach[END_REF] consider the mode-coupling mechanism whereas models proposed in [START_REF] Ding | Effects of rail dynamics and friction characteristics on curve squeal[END_REF][START_REF] Ding | An assessment of mode-coupling and falling-friction mechanisms in railway curve squeal through a simplified approach[END_REF] considers both falling friction and mode coupling.

For the modeling of rolling contact, there are mainly three kinds of model: point contact models, Kalker's type models with discretized contact surface and finite element models. Extensive reviews can be found for instance in [START_REF] Wriggers | Computational contact mechanics[END_REF][START_REF] Jacobson | Rolling contact phenomena[END_REF]. The first two types have been developed in the context of wheel/rail contact and have been investigated in detail by Kalker using special algorithms and formulations based on linear elasticity. Equivalent point-contact models are based on analytic formulas or heuristic laws of frictional rolling contact such as Hertz theory [START_REF] Hertz | Über die berührung fester elastischer körper[END_REF] for the normal problem and Kalker's linear theory [START_REF] Kalker | Wheel-rail rolling contact theory[END_REF], Vermeulen and Johnson's model [START_REF] Vermeulen | Contact of nonspherical elastic bodies transmitting tangential forces[END_REF], Shen et al. model [START_REF] Shen | A comparison of alternative creep force models for rail vehicle dynamic analysis[END_REF] or CHOPAYA law proposed by Ayasse et al. [START_REF] Ayasse | Handbook of railway vehicle dynamic -Chapter 4: Wheel-Rail Contact[END_REF] for the tangential problem. They are used in quasi-static computations with smooth surfaces, even in the case of large slips. In slightly nonlinear dynamic cases, close to pure rolling, equivalent point stiffness and damping can also be used and give good results, especially for wheel/rail rolling noise modeling [START_REF] Thompson | Railway noise and vibration: mechanisms, modelling and means of control[END_REF]. In other cases, including large slips and high-frequency dynamics, surface contact models are more adapted. They need a discretization of the contact zone. Kalker's simplified theory [START_REF] Kalker | A fast algorithm for the simplified theory of rolling contact[END_REF] implemented in the algorithm FASTSIM and Kalker's variational theory [START_REF] Kalker | Wheel-rail rolling contact theory[END_REF] implemented in CONTACT software are the most used models but other methods exist (cf. for instance Pieringer's model [START_REF] Pieringer | A numerical investigation of curve squeal in the case of constant wheel/rail friction[END_REF]) with slightly different assumptions or resolution techniques. However, some simplifications are generally performed in these models especially in the semi-analytical computation of the local contact flexibilities or influence functions (Boussinesq and Cerruti elastic half-space assumptions, contact/friction decoupling). The impact of these simplifications is rather unknown in the case of friction-induced vibrations characterized by large slips, high frequency dynamics and very fast evolutions at the contact zone. The alternative is to use finite element models. It makes it possible to consider also large deformations and plastic behaviors. The modeling of contact and friction by the finite element method is an active search field. Several variational formulations, discretization methods and specific non linear algorithms exist in the literature [START_REF] Wriggers | Computational contact mechanics[END_REF]. Lagrange multipliers or penalty methods are the most used techniques in order to take into account the constraints corresponding to contact and friction. Dynamics problems require special care since unilateral contact can lead to non-smooth behaviors like localized impacts between the structures [START_REF] Sinou | Mode coupling instability in friction-induced vibrations and its dependency on system parameters including damping[END_REF][START_REF] Charroyer | Parametric study of the mode coupling instability for a simple system with planar or rectilinear friction[END_REF][START_REF] Loyer | Study of nonlinear behaviors and modal reductions for friction destabilized systems. application to an elastic layer[END_REF]. In the case of rolling contact, it can be useful not to apply a Lagrangian description but to use a rotating reference frame, notably in order to optimize the meshes of the structures. Finite element treatment based on pure Eulerian or Alternate Lagrangian Eulerian (ALE) formulations for rolling contact problems can be found in [START_REF] Wriggers | Computational contact mechanics[END_REF][START_REF] Jacobson | Rolling contact phenomena[END_REF]. Most applications are quasi-static but such formulations could be advantageously used in high-frequency dynamics cases. Works dealing with wheel/rail contact analysis by FE methods are quite recent. Bogdansky et al. [START_REF] Bogdanski | Numerical stress analysis of rail rolling contact fatigue cracks[END_REF] solved the normal problem of static elastic wheel/rail contact in three dimensions with the FE method. The normal problem of contact of elastic-plastic material was investigated by Wiest et al. with a 3D static FE model [START_REF] Wiest | Assessment of methods for calculating contact pressure in wheel-rail/switch contact[END_REF]. Other solutions of the normal problem with FE can be found but the treatment of the tangential problem of wheel/rail rolling contact with FE is not frequent. Interesting results are given in the study of Toumi [START_REF] Toumi | Finite element analysis of the frictional wheel-rail rolling contact using explicit and implicit methods[END_REF] as well as the work of Zhao and Li [START_REF] Zhao | The solution of frictional wheel-rail rolling contact with a 3d transient finite element model: Validation and error analysis[END_REF] who treated wheel/rail rolling contact with an explicit FE method. By approaching a quasi-static state, the FE model is notably validated against Hertz theory and CONTACT in both normal and tangential solutions, for the case of wheel treadrail head contact. In these studies, dynamics and contact equations are solved in a Lagrangian frame. All the potential contact area in the rolling direction needs to be meshed with elements of sufficient small size to insure the desired precision. In addition, the meshes of the two potential contact surfaces are not compatible. As a consequence, the treatment of the contact is not straightforward. For instance, Wiest [START_REF] Wiest | Assessment of methods for calculating contact pressure in wheel-rail/switch contact[END_REF] and Toumi [START_REF] Toumi | Finite element analysis of the frictional wheel-rail rolling contact using explicit and implicit methods[END_REF] used a commercial software where contact between the two bodies is defined using a "masterslave" algorithm. Two types of contact discretization (node-to-surface and surface-to-surface) are used and the results differ depending on the user choice. Efficient algorithms exist for both formulations but, in any case, nonlinear transient FE simulations with frictional contact still remain very expensive in terms of CPU time and memory size when Lagrangian reference frame are considered.

Existing studies can also be distinguished according to the type of researched solutions. Two possible approaches can be used: complex eigenvalue analysis (CEA) and time-integration analysis. On the one hand, stability analyses look for solutions of the linearized dynamic equations governing the system [START_REF] Heckl | Curve squeal of train wheels, part 2: Which wheel modes are prone to squeal?[END_REF][START_REF] Van Ruiten | Mechanism of squeal noise generated by trams[END_REF][START_REF] De Beer | Squeal noise of rail-bound vehicles influenced by lateral contact position[END_REF][START_REF] Xie | Introduction of falling friction coefficients into curving calculations for studying curve squeal noise[END_REF][START_REF] Squicciarini | Curve squeal in the presence of two wheel/rail contact points[END_REF]. By solving a complex eigenvalue problem or by using the Nyquist criterion, they allow to determine the potential unstable modes and corresponding frequencies of the wheel/rail system and are useful to test the influence of the parameters on the onset of squeal. On the other hand, the dynamic nonlinearities related to the contact evolution at the frictional interfaces can be taken into account using time integration methods and the steady-state amplitudes as well as the full spectrum of squeal can be determined [START_REF] Heckl | Curve squeal of train wheels, part 1: mathematical model for its generation[END_REF][START_REF] Glocker | Curve squealing of trains: Measurement, modelling and simulation[END_REF][START_REF] Pieringer | A numerical investigation of curve squeal in the case of constant wheel/rail friction[END_REF][START_REF] Fingberg | A model of wheel-rail squealing noise[END_REF][START_REF] Périard | Wheel-rail noise generation: curve squealing by trams[END_REF][START_REF] Chiello | Curve squeal of urban rolling stock -part 3: Theoretical model[END_REF][START_REF] Huang | Squeal prediction for a bogied vehicle in a curve[END_REF][START_REF] Brunel | Transient models for curve squeal noise[END_REF]. As mentioned by many authors (cf. [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 instance), both approaches are needed. Indeed, the complex eigenvalue problem may overestimate or underestimate the number of unstable modes. The time integration method can also lead to the appearance of new frequency peaks related to harmonics of the fundamental frequencies of unstable modes. Finally, frequencies of the unstable modes could be different from those resulting from the time integration. The aim of this paper is to develop a full Finite Element (FE) computational method for the dynamics of frictional rolling contact systems in order to calculate reference solutions of complex systems, for curve squeal noise for example. Continuous equations of the problem are derived around the stationary position of rolling in an Eulerian reference frame. Local unilateral contact and Coulomb friction laws apply on the rolling interface.

Weak formulation and FE discretization lead to a nonlinear discrete system which is solved in the time domain by a numerical integration and a point-fixed algorithm at each time step. In addition to the transient approach, a stability analysis performed around the sliding equilibrium position allows to determine unstable modes and frequencies.

In order to obtain reasonable computation times, two reduction strategies are proposed. An application of the method to the frictional rolling contact between two annular cylinders is presented in both quasi-static and dynamic cases with mode coupling instabilities.The validation of the approach in quasi-static conditions is carried out by comparison with CONTACT software. Reduction strategies are validated by comparison with the full method.

The main contribution of the paper lies in the whole strategy proposed in order to efficiently compute dynamic transient solutions of rolling/sliding contact systems at high frequencies using finite elements. This strategy is based on the combination of existing techniques: an Eulerian reference frame, numerical tools adapted to the treatment of non-smooth frictional contact laws and reduction methods. The general features of the rolling contact formulations, the non-smooth contact formulations and the reduction methods are not new as part of the finite element method. However, their application to 3D high frequency dynamics of rolling systems with friction is less frequent, essentially due to the computational practical challenge of combining rolling phenomena, nonlinear frictional contact and time integration schemes. The paper provides some tools, which can help the community to compute numerical solutions for this kind of problem, especially in the case of friction-induced vibrations.

Finite Element (FE) formulation in the time-domain

Contact problem in the Eulerian frame

Two bodies in rolling contact are considered, as shown in Fig. 1. The conventions frequently used of rolling contact problems are adopted (cf. for instance [START_REF] Johnson | Contact mechanics[END_REF]). The z-axis is chosen to coincide with the common normal to the two surfaces in contact, the longitudinal x-axis corresponds to the rolling direction and the y-axis refers to the lateral direction. The whole structure is decomposed into two subdomains Ω = Ω 1 ∪ Ω 2 where subscripts 1 and 2 stand for the upper and the lower bodies respectively. In the absence of deformation and sliding, material particles of each surface move through the contact region in a direction parallel to x-axis with rolling speed V. In case of sliding, relative lateral velocities ∆V y = V y1 -V y2 and longitudinal velocities ∆V x = V x1 -V x2 have to be considered. The bodies may also have a relative angular velocity ∆ω z = ω z1 -ω z2 around their common normal (or spin).

∂Ω u f Ω 1 Ω 2 ∂Ω f1 S C Z (commun normal) Y Δω z ¡¢ u£ 1 ¡¢ f u 0 ΔVy ΔVx X Figure 1: Two bodies in rolling contact
According to [START_REF] Johnson | Contact mechanics[END_REF], in the Eulerian frame which moves with the point of contact, the relative sliding instantaneous velocities (or creep velocities) between the two deformable bodies at a fixed point of the potential contact interface S c are given by:

ṡx = v x1 -v x2 = ∆V x -∆ω z y + V ∂u x1 ∂x - ∂u x2 ∂x + (u x1 -ux2 ) ṡy = v y1 -v y2 = ∆V y + ∆ω z x + V ∂u y1 ∂x - ∂u y2 ∂x + (u y1 -uy2 ) (1)
where u(x, y, z, t) and v(x, y, z, t) denote respectively the displacement and velocity fields of the structure in the Eulerian frame, subscripts x and y stand for the spatial components of the fields and the notation u = ∂u ∂t refers to the time partial derivative. The terms involving V ∂u ∂x represent the deformation contributions due to rolling in the Eulerian frame whereas the terms involving u simply represent the dynamic contributions.

In addition, a displacement u d is prescribed on the part of the boundary of the domain δΩ u = δΩ u1 ∪ δΩ u2 and a surface load f s is applied on the part

δΩ f = δΩ f 1 ∪ δΩ f 2 .
In Ω, a volume load f d can be also considered.

Contact laws

To deal with unilateral contact on the interface, a non-regularized Signorini law is chosen:

∆u n -g ≤ 0 r n ≤ 0 (∆u n -g)r n = 0 (2)
where ∆u n = u n1 -u n2 is the normal relative displacement on the potential contact interface, r n is the normal contact stress and g is the initial gap. This law simply conveys that (i) there is no interpenetration between the two bodies, (ii) the interface only undergoes compression and (iii) it respects the condition of complementarity. A equivalent semi-regularized form of the law uses the projection on the negative real set (cf. for instance [START_REF] Jean | The non-smooth contact dynamics method[END_REF]):

r n = Proj R -(r n -ρ n ∆u n ) ∀ρ n > 0 (3)
where Proj R -(x) = min(x, 0) and ρ n is a positive real number called normal augmentation parameter.

To deal with frictional contact, a non-regularized Coulomb law with a constant friction coefficient µ is used:

r t ≤ -µr n r t = -µr n ⇒ ∃λ > 0, ṡt = -λr t r t < -µr n ⇒ ṡt = 0 ( 4
)
where ṡt is the vectorial creep velocity with components ṡx and ṡy defined in equation ( 1) and r t is the tangential stress due to friction. This law states that the frictional reaction cannot be greater than a limit -µr n . If it reaches the limit, the material particle slides and has a direction opposed to the direction of the reaction (slip state). In other cases, the relative velocity is null (stick state).

A equivalent semi-regularized form of the law uses the projection on the Coulomb cone (cf. for instance [START_REF] Jean | The non-smooth contact dynamics method[END_REF]):

r t = Proj C (r t -ρ t ṡt ) ∀ρ t > 0 ( 5
)
where Proj C (x) = min( µ|r n | x , 1)x and ρ t is a positive real number called tangential augmentation parameter.

Dynamics equations and FE discretization

A major consequence of the Eulerian frame is to induce convective terms in the equations of motion. These terms are only significant in the case where the rotational speed is in the order of magnitude of the eigenfrequencies ( [START_REF] Chambrette | Stability of a beam rubbed against a rotating disc[END_REF]). It is assumed that it is not the case in the present application and the effects of the convective terms are neglected in the following developments.

The displacement field u must verify the equations of continuum mechanics [START_REF] Raous | Numerical characterization and computation of dynamic instabilities for frictional contact problems[END_REF] together with the local formulations of contact equations Eqs. ( 3) and [START_REF] Heckl | Curve squeal of train wheels, part 2: Which wheel modes are prone to squeal?[END_REF]. In order to introduce FE approximations, the principle of virtual power is used for the system dynamics and the contact laws are written in weak forms [START_REF] Kudawoo | Problèmes industriels de grande dimension en mécanique numérique du contact: performance, fiabilité et robustesse[END_REF]:

find u ∈ U and r such as ∀u * ∈ U o and ∀r * Ω ρu * . ü dΩ + Ω (u * ) : σ(u) dΩ = Ω u * .f s dΩ + ∂Ω F u * .f d ds + S c u * .r ds S c r * n r n ds = S c r * n Proj R -(r n -ρ n ∆u n ) ds S c r * t r t ds = S c r * t Proj C (r t -ρ t ṡt ) ds (6)
where (u) and σ(u) stand respectively for the symetric gradient and the stress tensor corresponding to any displacement field u, the notation ü = ∂ 2 u ∂t 2 refers to the double time partial derivative,

U = {u|u = u d on ∂Ω u } and U 0 = {u|u = 0 on ∂Ω u }.
By using appropriate shape interpolation functions for unknown and test displacement and reaction fields (u, r, u * , r * ), finite element discretization of Eq. ( 6) directly gives:

M Ü + C U + KU = F + P n T R n + P t T R t R n = S c N n T Proj R -N n H -1 n R n -ρ n (P n U -G) ds R t = S c N t T Proj C N t H -1 t R t -ρ t Ṡt (P t U, P t U) ds (7)
where U, R n and R t denote respectively the vectors of nodal displacements, equivalent normal reactions and equivalent friction forces, M, C, K are respectively the mass, damping and stiffness matrices of the structure without contact, G is the vector of nodal initial gaps and N n , N t are the shape function vectors on the contact interface.

In addition, P n , P n are matrices allowing to pass the contact reactions from the local relative frame to the global frame whereas

H n = S c N n T N n ds and H t = S c N t
T N t ds are transformation matrices from nodal to equivalent forces. Finally, Ṡt (P t U, P t U) denotes the vector of nodal creep velocities which can be determined linearly from local displacement and velocity vectors, taking into account quasi-static creep velocities (cf. Eq. ( 1)).

Reduction strategies

As the size of the system is often large and the nonlinear solving process implies several resolutions of a linear system at each time step, reducing the size of the system is necessary to obtain reasonable computation times. The principle is to search an approximated solution U = Bq r of the problem spanned by a reduced basis B which leads to a reduced dynamics equation:

M r qr + C r qr + K r q r = B T (F + P n T R n + P t T R t ) (8)
where M r , C r , K r = B T (M r , C r , K r )B and the size of the system is reduced to the number of modes in basis B.

In the field of friction-induced vibrations, the use of a priori reduction bases remains a commonly used technique for industrial systems but the performance of the bases depends on the application. In the case of a simplified brake squeal application, Brizard et al. proposed several optimized reduction bases for the stability analysis [START_REF] Brizard | Performances of some reduced bases for the stability analysis of a disc/pads system in sliding contact[END_REF].

For the same application, Lorang and Chiello [START_REF] Lorang | Stability and transient analysis in the modelling of railway disc brake squeal[END_REF] showed a strong participation of the unstable complex mode shapes determined by the stability analysis in the nonlinear solution, suggesting that reduction bases based on unstable complex modes could be used in the nonlinear problem. However, Sinou [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] observed that the contributions of the harmonic components and the combination of frequency components can not be neglected. Loyer et al. [START_REF] Loyer | Study of nonlinear behaviors and modal reductions for friction destabilized systems. application to an elastic layer[END_REF] performed intensive simulations in order to test two types of bases. The degrees of freedom at contact interface are kept in these bases based on the complex modes resulting from the stability analysis. The study notably shows that an accurate model reduction for friction-destabilized system is possible but requires great care. In particular, all the stable modes in the frequency range must be included to ensure an accurate prediction of the nonlinear response. In addition, if separation at the interface is present, the reduction basis must also contain corresponding static boundary modes.

It should also be noted that some other performant techniques exist for the reduction of non linear dynamic systems. To determine relevant reduction bases, it is thus possible to use singular value decomposition (SVD)

of Proper Orthogonal Decomposition [START_REF] Steindl | Methods for dimension reduction and their application in nonlinear dynamics[END_REF] of some representative sets of the full temporal solution to extract the major shapes. The concept of nonlinear modes [START_REF] Kerschen | Nonlinear normal modes, part i: A useful framework for the structural dynamicist[END_REF] was also developed in recent years and proves to be efficient in some applications. For autonomous systems with mono-instabilities (only one unstable mode), the benefit is small since the searched periodic limit cycles are precisely the nonlinear modes.

In this paper, the proposed basis is rather classical and is based on an a priori selection of free-interface normal modes enriched by static local modes.

Component Mode Synthesis (CMS) with free-interface modes

A first reduction basis includes free-interface normal modes of the structure and static attachment modes B = [Φ Φ s ]. Matrix Φ contains the real solutions of the free and undamped system:

(K -ω 2 M)U = 0 (9)
whereas Φ s make the base statically complete by adding static solutions to unitary forces (normal and tangential) on the contact interface:

KΦ s = P c T ( 10
)
where

P c T = [P n T P t T ]
. In order to improve the numerical efficiency, the basis is then orthogonalized in relation to stiffness matrix K before use. It is well-known that static attachment modes are useful to provide a better approximation of the dynamics by comparison with normal modes alone. It is also worth noticing that the size of the basis is n m + 3n c where n m is the number of normal modes taken into account in the frequency range of interest and n c is the number of facing nodes on the interface.

Contact static approximation

A stronger approximation consists in neglecting the dynamic terms in the reduced equations corresponding to attachment modes. However, bases Φ and Φ s are not orthogonal which is a necessary condition in order to separate the contributions of normal modes and attachment modes. This can be fixed by using residual attachment modes defined by:

Φs = Φ s -Φ(Φ T KΦ) -1 Φ T P c T (11)
These modes are the static displacement responses to unit contact reactions after the elimination of the contribution of normal modes. It can be easily verified that ΦT s KΦ = 0. The use of basis B = [Φ Φs ] in Eq. ( 8) together with the elimination of the dynamic terms relating to attachment modes gives:

Φ T MΦ q + Φ T CΦ q + Φ T KΦq = Φ T (F + P n T R n + P t T R t ) ΦT s K Φs q s = ΦT s (F + P n T R n + P t T R t ) (12)
where q and q s are the generalized coordinate vectors corresponding respectively to normal modes and residual attachment modes such that U = Φq + Φs q s . A total decoupling between the generalized equations corresponding to free-interface normal modes and residual static attachment modes is obtained.

In addition, noticing that:

ΦT s K Φs q s = Φ T s -P c Φ(Φ T KΦ) -1 Φ T K Φ s -Φ(Φ T KΦ) -1 Φ T P c T q s = P c Φ s -P c Φ(Φ T KΦ) -1 Φ T P c T q s = P c Φs q s = P c (U -Φq) (13)
and including the above expression in Eq. ( 12) gives:

Φ T MΦ q + Φ T CΦ q + Φ T KΦq = Φ T (F + P n T R n + P t T R t ) P c U = ΦT s (F + P n T R n + P t T R t ) + P c Φq (14)
It must be emphasized that, with this formulation, contact displacements P c U can be directly calculated from expanded forces and normal modes coordinates.

This reduction strategy consists in solving the global dynamics using free-interface normal modes and adding a local static residual flexibility controlled by matrix ΦT s in the expression of the contact displacements. This allows to determine the contact displacements with a high precision level which is needed to solve the frictional contact laws and calculate the reaction forces. This method gets very close to models such Pieringer's one [START_REF] Pieringer | A numerical investigation of curve squeal in the case of constant wheel/rail friction[END_REF] combining modal dynamics and static local resolutions of the kind of CONTACT software. However, in the proposed method, the local flexibility is calculated numerically by finite elements instead of using semi-analytical influence functions derived from Boussinesq's approximations. In particular, it takes into account the static normal/tangential coupling at the contact zone. It also depends on the number of free-interface modes in the basis since it corresponds to a residual flexibility.

Stability analysis

The aim of the stability analysis is to address the mechanism of self-excited vibration due to frictional contact through the determination of the evolution of small perturbations around the steady sliding equilibrium. In stable cases, the perturbations vanish and no vibration occur. In unstable cases, some perturbations tend to diverge which can lead to self-sustained vibrations. Such an analysis is performed by a linearization of the non-linear equations around the equilibrium.

Quasi-static equilibrium

The quasi-static or steady sliding equilibrium is first obtained by neglecting the dynamic terms in Eqs. [START_REF] Rudd | Wheel/rail noise-part ii: Wheel squeal[END_REF] and [START_REF] Pieringer | A numerical investigation of curve squeal in the case of constant wheel/rail friction[END_REF]:

KU e = F + P T n R e n + P T t R e t R e n = S c N T n Proj R -N n H -1 n R e n -ρ n (P n U e -G) ds R e t = S c N T t Proj C N t (H -1 t R t -ρ t Ṡe t ) ds (15)
where Ṡe t = Ṡt (P t U e , 0) is the vector of nodal quasi-static creepage velocities. Assuming that some solutions of Eq. ( 15) exist and can be calculated, it notably provides the status of the nodes on the contact interface as a function of the equivalent normal reactions r e n and friction forces r e t :

• if r e n = 0, the facing nodes are not in contact,

• if r e n 0 and r e t < -µr e n the facing nodes are sticking,

• if r e n 0 and r e t = -µr e n the facing nodes are sliding ; the sliding direction of the friction force is then given by t = r e t / r e t .

Complex Eigenvalue Analysis (CEA) in case of full steady sliding

In this paper, stability analysis is only carried out in the case of full steady sliding (no sticking region) and maintained contact configuration: it is thus assumed that for each node in contact at equilibrium, bilateral contact and sliding Coulomb friction laws apply. In order to perform the stability analysis, these laws have to be linearized.

On effective contact region, the linearized forms of Eqs. ( 3) and ( 5) with the above assumptions can be written (cf. for instance [START_REF] Brizard | Performances of some reduced bases for the stability analysis of a disc/pads system in sliding contact[END_REF]):

∆u n = 0 r t = -µr n t -c b ub b (16)
where c b = -µr e n / ṡe t is a damping term due to the linearisation of the sliding direction of the friction force and b is the tangential direction orthogonal to t.

Searching a discrete solution of the form U e + Ũ exp(λt) where Ũ stands for the complex displacement vector corresponding to small harmonic perturbations around the equilibrium, the linearized form of equation ( 7) leads to a constrained non symmetric eigenvalues problem:

λ 2 M + λ(C + C b ) + K Ũ = ( PT n + µ PT t ) Rn Pn Ũ = 0 (17)
where C b is the damping matrix provided by the linearisation of the sliding direction of friction force [START_REF] Loyer | Study of nonlinear behaviors and modal reductions for friction destabilized systems. application to an elastic layer[END_REF][START_REF] Brizard | Performances of some reduced bases for the stability analysis of a disc/pads system in sliding contact[END_REF][START_REF] Lorang | Tgv disc brake squeal[END_REF] and Pn , Pt are new projection matrices such that Pn is the restriction of P n on nodes in the effective contact region at equilibrium whereas Pt is the restriction of P t on components in direction t on nodes in the effective contact region.

Complex modes and eigenvalues of the problem are then calculated. Modes corresponding to eigenvalues with positive real part are unstable. The divergence rate of a mode is notably defined as Re(λ)/Im(λ) where (Re(λ), Im(λ)) are respectively the real and imaginary parts (pulsation) of the mode. This rate corresponds to a negative damping rate.

Reduced CEA formulations

Solving such a large non-symmetric eigenvalue problem needs model reduction. Several reduction bases have been notably tested by following the methodology proposed by Brizard [START_REF] Brizard | Performances of some reduced bases for the stability analysis of a disc/pads system in sliding contact[END_REF]. The two reduction strategies proposed in the nonlinear case (section 3) are adapted here for the stability analysis.

Component Mode Synthesis with free-interface modes

The first reduced basis includes free-interface normal modes of the structure and static attachment modes B = [Φ Φ s ] as proposed in section 3 but defined on effective contact region by using Pn and Pt instead of P n and P t . In order to solve the eigenvalue problem with this reduced basis, the constraints are directly added to the generalized coordinates associated to the modes (cf. [START_REF] Craig | On the use of attachment modes in substructure coupling for dynamic analysis[END_REF][START_REF] Craig | A review of time-domain and frequency-domain component mode synthesis method[END_REF] for symmetric normal problems and [START_REF] Loyer | Study of nonlinear behaviors and modal reductions for friction destabilized systems. application to an elastic layer[END_REF][START_REF] Brizard | Performances of some reduced bases for the stability analysis of a disc/pads system in sliding contact[END_REF] for non-symmetric friction problems).

Contact static approximation

The stronger approximation consisting in neglecting the dynamic part of the attachment modes is also tested for the stability analysis. The same technique used in section 3 is applied. The reduction is first performed with basis B = [Φ Φs ] composed of free-interface component modes Φ and static attachment residual modes Φs defined on effective contact region. The dynamic effects corresponding the residual attachment modes are then neglected leading to the following constrained eigenvalue problem:

λ 2 Φ T MΦ + λΦ T (C + C b )Φ + Φ T KΦ q = Φ T ( PT n + µ PT t ) Rn Pc Ũ = ΦT s ( PT n + µ PT t ) Rn + Pc Φ q Pn Ũ = 0 ( 18
)
where q is the generalized vector corresponding to free-interface normal modes.

From the two last lines of Eq. ( 18), normal reactions Rn can be expressed as a function of the generalized vector q:

Rn = -I n ΦT s ( PT n + µ PT t ) -1 Pc Φ q (19)
where I n is the Boolean localization matrix such that Pn = I n Pc . Finally, using the above expression in the reduced eigenvalue problem gives:

λ 2 Φ T MΦ + λΦ T (C + C b )Φ + Φ T (K + K c )Φ q = 0 (20)
where K c is a non-symmetric stiffness matrix taking account the effect of the local residual flexibility of the structure due to normal and friction forces:

K c = ( PT n + µ PT t ) I n ΦT s ( PT n + µ PT t ) -1 Pn (21)
This reduction strategy consists in solving the global dynamics using free-interface normal modes and adding a local static residual stiffness controlled by matrix K c . This matrix plays the same role than 2-DOF Hertzian stiffness/Coulomb friction relations in the case on equivalent point contact models but is calculated by finite elements.

Numerical methods for nonlinear simulations

For the determination of quasi-static and transient solutions from Eqs. ( 7) and ( 15), with or without reduction, the following numerical methods are used.

Integration scheme for the dynamics

For the computation of the transient solution, the chosen time integration method is a modified θ-method. This is a first-order scheme developed by Jean [START_REF] Jean | The non-smooth contact dynamics method[END_REF] and appropriate to unilateral contact dynamics. It notably allows to compute quasi-inelastic shocks. A value θ = 0.5 is chosen to avoid numerical damping. The detail of the scheme with an application to friction-induced vibrations can be found for instance in Loyer's work [START_REF] Loyer | Study of nonlinear behaviors and modal reductions for friction destabilized systems. application to an elastic layer[END_REF].

Non-linear resolution

For the computation of the quasi-static solution and dynamics solutions at each time step, an iterative fixed point algorithm on equivalent contact reactions and friction forces is used with a stop criterion based on forces convergence [START_REF] Loyer | Study of nonlinear behaviors and modal reductions for friction destabilized systems. application to an elastic layer[END_REF]. This algorithm is appropriate to the formulation of the frictional contact laws as nonlinear projections. The main advantage of the fixed point algorithm is that the integrator matrix remains constant at each iteration. In the simulations presented in this paper, the augmented parameters ρ n , ρ t are chosen as the smallest eigenvalue of the integrator matrix condensed on the contact degrees of freedom [START_REF] Khenous | Hybrid discretization of the signorini problem with coulomb friction. theoretical aspects and comparison of some numerical solvers[END_REF]. More optimized parameters can be found in [START_REF] Loyer | Study of nonlinear behaviors and modal reductions for friction destabilized systems. application to an elastic layer[END_REF] but have not been used.

Results

In this section, an application of the methodology presented in the previous sections is proposed for the case of rolling contact of two annular cylinders with lateral creepages. The model is first described. Second, the model is validated in quasi-static conditions with CONTACT software. Considering that there is no reference approach in the dynamic case, this quasi-static validation is essential in order to show that the proposed method is able to accurately reproduce the creep phenomenon characterized by an effective contact area separated in two (stick and slip) zones. A stability analysis is then carried out. Transient calculations are finally performed with and without reduction and results are discussed and compared.

Description of the model

A rolling contact between two same annular cylinders is considered, as shown in Fig. 2. The material behavior is assumed to be linearly elastic, isotropic and undergoing small deformations. A Rayleigh damping with coefficients α and β is introduced in the model. The values of α and β have been chosen in order to obtain modal damping factors comprised between 0.5 % and 5 % in the frequency range of interest of [100-5000] Hz. These quite high values are useful to show that the found instabilities are strong and do not disappear with any small damping in the structures (cf. section 6.4). The material and geometrical parameters of the cylinders are listed in Tab. 1.

A vertical displacement u z0 is applied at the hub of the upper cylinder whereas a rigid constraint is applied at the lower cylinder. For the discretization by finite elements, compatible meshes on the interface are considered: facing nodes of the two cylinders in potential contact have identical tangential coordinates x and y. The global meshes of the cylinders have been chosen in order to respect the dynamic criterion of 10 elements per wavelength for the free normal modes in the frequency range of interest. The meshing of the contact zone have to allow an accurate description of the slip/stick transition in the area. Consequently, the size of the elements in the contact zone have to be much smaller than the dimensions of the effective area (at least an order of magnitude). On the other hand, the theory gives the dimensions of (4.5, 4.5) mm for the contact area. The convergence of the quasi-static and dynamic solutions gives an optimal size of 1 mm for the elements at the contact zone as shown in Fig. 3. Actually, in the proposed application, the two cylinders are fully symmetrical and it would have been possible

to model only one cylinder using symmetry relations.

Free-interface normal modes

The 80 first natural frequencies and corresponding free-interface modes have been calculated up to 7510 Hz.

The eigenmodes are classified according to their predominant deformations (axial, radial or circumferential) and their numbers (n, m) of nodal diameters and nodal circles. Two modes are represented in Fig. 4. This classification is necessary to determine the optimal mesh according to a targeted number of element per wavelength for the free normal modes. However, no convergence study has been performed in order to examine the effect of these meshes on the full model with frictional contact (nor for stability analysis either for self-sustained vibrations).

Quasi-static results

In this section the quasi-static rolling contact of the two cylinders with lateral creepage is considered. The contact zone is supposed to be laterally centered. As shown in Fig. 2b, the geometry of the lateral profiles of the two cylinders has been curved in order to obtain an effective 3-dimensional contact. For each cylinder the applied radius of curvature is 0.5 m. The rolling is performed in the -x direction with V = 10 m/s. An vertical displacement u z0 = 0.1 mm is applied leading to a resultant vertical contact force of about 75 kN. A friction coefficient µ = 0.3 is considered. Longitudinal and spin creepages are set to zero (∆V x = ∆ω z = 0).

For an imposed lateral creepage ∆V y /V = 0.3%, the longitudinal distributions of the normal contact stresses P n and frictional stresses P t on the center line of the contact zone (y = 0) are presented in Fig. 5. As expected, a stick zone occurs at the leading edge of the contact and a slip zone occurs at the trailing edge of the contact.

The comparison of the results obtained with the proposed full FE method and the results provided by CONTACT software shows a good agreement with some small differences for the tangential stresses in the transition between the stick zone and the slip zone.

-0.01 -0.005 0 0.005 0.01

x(m) As in the finite element method proposed, Kalker's variational theory [START_REF] Kalker | Wheel-rail rolling contact theory[END_REF] implemented in CONTACT software is based on a discretization of the contact zone combined with unilateral contact and Coulomb friction laws.

The differences mainly lie in the two simplifications proposed by Kalker in order to calculate the distribution of stresses and deformations in the contact area: local flexibilities computed by using Boussinesq and Cerruti (static) elastic half-space assumptions and decoupling between the normal and tangential problems. Another difference concerns the interpolation functions in the contact zone (linear in the proposed method versus constant in CON-TACT software). This can explain the small differences between the results obtained with both methods. As fewer approximations are involved in the proposed method than in Kalker's theory, the results obtained with the finite element method are supposed to be more accurate. However, in such a quasi-static configuration, the differences can simply come from the interpolation functions.

Stability results

Stability analysis is carried out in case of full sliding using a higher creepage ∆V y /V = 1% by using 160 free-interface normal modes (80 for each cylinder). In order to solve the non symmetric eigenvalue problem, the two reduction strategies presented in section 4.3 (classical CMS and contact static approximation) are tested and compared with reference results obtained with an iterative method minimizing the force residue, like in [START_REF] Brizard | Performances of some reduced bases for the stability analysis of a disc/pads system in sliding contact[END_REF].

In the case where the contact point is centered, no instability is found whatever the friction coefficient. In order to obtain instabilities the contact point is laterally shifted, as shown in Fig. 6. This offset increases the lateral/vertical coupling which is a necessary condition to obtain mode-coupling instabilities. The lateral position of the contact point seems to be one of the most important parameter for the occurrence of instabilities. Indeed, it bring out a stronger coupling between the normal and tangential dynamics due to the asymmetrical geometry of the system compared with the centered position.

Divergence rate (%) In order to confirm the mode-coupling instability and to test the behavior of the reduced basis presented in section 4.3, the evolution of the eigenvalues corresponding to the two complex modes involved in the instability (i.e.

the bifurcation curves) are represented in Fig. 9 without damping (by imposing C = C b = 0) and in Fig. 10 with damping. As expected, the increase of the friction coefficient brings the two modes closer until their frequencies coalesce and one mode becomes unstable.

In order to analyze the case with damping, a brief discussion about the known effects of damping on coalescence patterns is needed. Indeed, these effects are not straightforward as shown for instance [START_REF] Sinou | Mode coupling instability in friction-induced vibrations and its dependency on system parameters including damping[END_REF][START_REF] Hoffmann | Effects of damping on mode-coupling instability in friction induced oscillations[END_REF][START_REF] Fritz | Investigation of the relationship between damping and mode-coupling patterns in case of brake squeal[END_REF]. In particular, a non-uniform repartition of damping between the modes involved in the coalescence can destabilize the system in some cases, especially trough a smoothing effect. In the present system, the internal viscous damping in the structures is proportional and the two modes involved in the instability are equally damped at µ = 0. Thus, in this special case, only a lowering effect is expected, which tends to stabilize the mechanical system by shifting vertically the bifurcation curves in comparison with the undamped case. This can be observed in Fig. 10b. To be rigorous, it should be noted that the linearization of the planar friction force in Eq. ( 16) leads to an added damping effect (matrix C b ), which is nor proportional either equally distributed between the modes or could be a smoothing factor as show in [START_REF] Charroyer | Parametric study of the mode coupling instability for a simple system with planar or rectilinear friction[END_REF]. However, in the present application, this effect seems to be negligible compared to the high internal damping. Concerning the evolution of the frequencies in Fig. 10a, the results also comply with previous observations, showing a dispersion of the curves after the merging point.

With regard to reduction strategies, the first validation has been performed by comparing the real and imaginary parts of complex pulsations provided by the different bases in the CEA for the reference friction coefficient µ = 0.3.

As recommended in Brizard's work, the highest frequency of the modes taken into account has to be equal to 1.5

times the upper frequency of the frequency range of interest, i.e. 7546 Hz. The number of free normal modes in the CMS and CSA bases are chosen with this criterion. Thus, 80 free normal modes per disc have been retained, i.e. 160 modes in all. To go further, the validation has been done for µ = 0 to 0.6 but only for the complex pulsations of the two modes involved in the instability (bifurcation curves). The relative errors on the frequencies and the absolute errors on the divergence rates of complex modes obtained with the two proposed reduction bases in comparison with the results obtained with the iterative method are presented in Figs. Frequency, b: Divergence rate)

Transient results

Transient results corresponding to the unstable case founded in the previous section are determined using a numerical time integration from given initial conditions. In all the following results, the integration starts from the equilibrium i.e. the initial displacements are the displacements obtained from the quasi-static solution and the initial velocities are null. At this point, it is important to recall that the initial conditions can influence the final limit cycle. However, it is especially true in the case of several unstable modes as shown for instance in [START_REF] Loyer | Study of nonlinear behaviors and modal reductions for friction destabilized systems. application to an elastic layer[END_REF].

When only one mode is unstable as in the present case, it is also true but particularly when initial conditions are The time series of the lateral velocity of a point outside the contact zone with coordinates (-0.01, 0, -0.01),

i.e. at 1 cm from the contact zone in the vertical and longitudinal directions, is presented in the Figure 17a. The 360 spectrogram in Fig. 17b allows to observe which frequencies are present in the solution and when they appear. A fundamental frequency is founded at f 0 = 1612.9 Hz which gets very close to the frequency of the unstable complex mode provided by stability analysis (1609 Hz). In the steady-state response, harmonic frequencies f k = k f 0 appear. Although in this case, the fundamental frequency is not so different from the frequency of the unstable Tangential contact resultant force (kN) Figure 16: Status of nodes in the potential contact zone: no contact zone (red), slip zone (blue) and stick zone (green). The Arabic notations "1", "2", "3", "4", "5" and "6" refer to Fig. 15 where the points are equitably distributed over time in the limit cycle complex mode, these results show the advantages of the nonlinear time integration method compared with the stability analysis. Indeed, the amplitude of the vibrations and the harmonic frequencies cannot be determined by the complex eigenvalue analysis.

Reduced models

The two reduction strategies presented in section 3 are tested using 160 free-interface normal modes in the frequency range up to 7510 Hz and 1323 residual attachment modes corresponding to 1323 degrees of freedom on the potential contact interface. For the classical CMS method, a reduced integrator matrix of size 1473 × 1473 is then obtained while the size of the initial matrix is 239358 × 239358. This reduced matrix is full but may be factorized before the iterative computation. The problem is therefore reduced significantly. For the contact static approximation, the size of the integrator matrix is reduced to 160 × 160 which induces a high gain in computation times.

Concerning the precision of the reduced models, the lateral velocity of the point outside the contact zone with coordinates (-0.01, 0, -0.01) obtained with the reduced basis are presented in Figs. acceleration or a deceleration of the transient response [START_REF] Loyer | Study of nonlinear behaviors and modal reductions for friction destabilized systems. application to an elastic layer[END_REF]. In the present case, the reduced solutions are clearly accelerated, which explains the differences in amplitudes. In the steady-state response (Fig. 19b), the fundamental frequencies obtained with both reduction bases are 1615.5 Hz. The corresponding error (0.16 %) is small. In addition, the amplitude of the three solutions are very similar, with an error of 2.3% on the maximum amplitude.

The corresponding spectra are also given in Fig. 20 showing a good agreement between the reduced solutions and the full solution: an error of 0.15 dB at the fundamental frequency and errors of 0.02, 0.73 and 1.12 dB at the harmonic frequencies.

In addition, it should be noted that the computation times resulting from the use of the reduction bases CMS and CSA are respectively ten and fifteen times shorter than for the computation of the full solution. Although the performance of the reduction bases are rather impressive, it is difficult to generalize this result to other cases or applications since it strongly depends on the frequencies involved in the nonlinear solution.

In the present case, the nonlinear effects seem to be localized to the contact area, which is coherent with the approximations induced by the two reduction strategies. The amplitude of the harmonics in the velocity spectrum of the point outside the contact zone are quite small compared with the amplitude of the peak at the fundamental frequency. For cases with a greater contribution of harmonics, it might be necessary to complete the basis with normal modes at higher frequencies.

To be complete, it should also be noted here that several non linear methods exist for the direct computation of the steady-state solution, even if they have not been used much in the case of rolling contact systems. In particular, when only one mode is unstable, nonlinear frequency-domain methods like for instance the constrained harmonic balance method (CHBM) proposed by [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] could undoubtedly be an efficient alternative to the time integration method. Other techniques like the shooting method or the orthogonal collocation method could also be used [START_REF] Steindl | Methods for dimension reduction and their application in nonlinear dynamics[END_REF][START_REF] Kerschen | Nonlinear normal modes, part i: A useful framework for the structural dynamicist[END_REF][START_REF] Charroyer | Self-excited vibrations of a non-smooth contact dynamical system with planar friction based on the shooting method[END_REF]. However, the main objective of the paper does not concern the performance of the time integration scheme but the use of a full finite element approach to study the dynamics of rolling/sliding contact systems at high frequencies, notably in order to calculate reference solutions. Solutions provided by time integrations are often used as reference solutions even when effective nonlinear methods like the CHBM can be used. Moreover, the time integration method has the advantage of providing the transient part of the solution and allows the steadystate solution to be linked with initial conditions. Finally, a third reason for using the time integration method is that it could be used without adaptation for systems exhibiting several unstable modes in the complex eigenvalue analysis, which is not a rare case, for instance in curve squeal problems.

Conclusion

In this paper, a method is proposed for the modeling and analysis of the high-frequency friction-induced instabilities of two structures in rolling contact in time and frequency domains. A full finite element formulation around the stationary position in an Eulerian reference frame is derived with a fine discretization of the contact surface combined with unilateral and Coulomb friction laws with constant friction coefficient. Appropriate numerical techniques are used in order to solve the nonlinear discrete equations in quasi-static or dynamic conditions. In addition to the transient approach, a stability analysis performed around the full sliding equilibrium position allows to determine unstable modes and frequencies.

In order to reduce the computational effort, reduction strategies are proposed for both domains. The first technique uses a classical CMS reduced basis including free-interface normal modes and static attachment modes.

A second technique consists in simply adding a residual static contact flexibility to the free-interface normal modes when solving the frictional contact equations (contact static approximation). It significantly differs from classical reduction bases used in previous papers and is well suited to this kind of problem where the contact area is small compared to the structures and proves to be efficient in computation cost.

The method is tested in the case of frictional rolling contact between two annular cylinders. The quasi-static results show a good agreement with the ones obtained with Kalker's CONTACT software. In case of full sliding, the stability analysis brings out a mode coupling instability when the contact zone is laterally shifted from the center of the cylinders. In the unstable configuration, the numerical integration provides solutions in the time domain which are coherent with the stability results. In particular, the unstable frequency determined by this analysis is very close to the fundamental frequency observed in the transient solution. The status of the nodes at the contact zone in the self-sustained vibrations clearly show slip-stick limit cycles.

Concerning the performance of the reduction strategies, the approximate results obtained for the stability analysis show a good agreement with the reference ones. In the time domain, the amplitudes of the reduced solutions are slightly greater than the full solution in the linear phase. However in the steady-state phase, solutions are very similar. This allows to considerably reduce the computational times, especially when the contact static approximation is used.

This methodology may be used later in a large variety of complex friction-mechanical system, especially curve squeal.

Figure 2 :

 2 Figure 2: Contact between two annular cylinders with lateral curvatures and contact point in the center (a:Front view with zoom on the contact zone, b: Side view)

Figure 3 :

 3 Figure 3: FE model of a cylinder with fine mesh on the contact zone

Figure 4 :

 4 Figure 4: Examples of free eigenmodes of a cylinder (a: Axial mode (0,2) (316 Hz), b: Radial mode (0,3)(3695 Hz))

Figure 5 :

 5 Figure 5: Quasi-static frictional stresses along the line y = 0 for ∆V y /V = 0.3% and µ = 0.3

Figure 6 :

 6 Figure 6: Lateral offset of the contact point

Figure 7 :

 7 Figure 7: Divergence rate of complex eigenvalues for ∆V y /V = 1% and µ = 0.3

Figure 8 :Figure 9 :Figure 10 :

 8910 Figure 8: Unstable mode for ∆V y /V = 1% and µ = 0.3

11 and 12 .Figure 11 :

 1211 Figure 11: Relative errors of the frequency and absolute errors on the divergence rate of complex modes without damping for ∆V y /V = 1% (a:

Figure 12 :

 12 Figure 12: Relative errors of the frequency and absolute errors on the divergence rate of complex modes with damping for ∆V y /V = 1% (a: Frequency, b: Divergence rate)

Figure 13 :

 13 Figure 13: Evolution of the normal contact resultant force

Figure 14 :Figure 15 :

 1415 Figure 14: Evolution of the tangential contact resultant force

Figure 17 :

 17 Figure 17: Tangential velocity of a point outside the contact zone with coordinates (-0.01, 0, -0.01) (a: Temporal evolution, b: Spectrogram (time-frequency))

  18 and 19 in comparison with the full solution. The fundamental frequencies of the three solutions are very close to each other. With the same initial conditions, the amplitudes of the reduced solutions in the transient response are greater than the full solution. It is well-known that the alteration of the transient evolution of modal participations in the transient state can lead to an

Figure 18 :Figure 19 :

 1819 Figure 18: Lateral velocity in a point outside the contact zone with coordinates (-0.01, 0, -0.01) obtained with the reduced basis (a: Component Mode Synthesis, b: Contact static approximation)

Figure 20 :

 20 Figure 20: Power spectrum (dB re 1 m/s) of the steady-state lateral velocity in a point outside the contact zone with coordinates (-0.01, 0, -0.01) obtained with the reduced basis

Table 1 :

 1 Material and geometrical parameters of the cylinders size of the problem can become very large when the number of contact dofs increases. In this application, Hertz

	Young's modulus	205 GPa
	Poisson's ratio	0.3
	Density	7800 kg/m 3
	Rayleigh's damping coefficients (α, β)	(1, 10 -6 )
	Internal radius	0.1 m
	External radius	0.5 m
	Thickness	0.05 m

Acknowledgements

This work was carried out within the CERVIFER project backed by the ADEME organization and the Haut-de-France region. It was also supported by the LabEx CeLyA of Université of Lyon, operated by the French National Research Agency (ANR-10-LABX-0060/ANR-11-IDEX-0007). The authors gratefully acknowledge the support of these institutions.