L Charroyer

O Chiello

J-J Sinou

Self-excited vibrations of a non-smooth contact dynamical system with planar friction based on the shooting method

Self-excited vibrations of a non-smooth contact dynamical system with planar friction based on the shooting method

Lucien Charroyer, Olivier Chiello, Jean-

Introduction

Despite the large amount of past and present research on the subject, squeal noise is still a current issue for engineering and scientific communities. Over recent decades, many investigations considering various specific features and developments have been proposed to study the problem of friction-induced vibration. Overviews can be found for instance in [START_REF] Kindkaid | Automotive disc brake squeal[END_REF][START_REF] Ouyang | Numerical analysis of automotive disc brake squeal: a review[END_REF][START_REF] Spurr | A theory of brake squeal[END_REF].

Regarding the squeal noise occurring in several engineering applications, different models and mechanisms have been proposed and developed [START_REF] Ibrahim | Friction-induced vibration, chatter, squeal and chaos part ii: Dynamics and modelling[END_REF][START_REF] Ibrahim | Friction-induced vibration, chatter, squeal and chaos part 1: mechanics of contact and friction[END_REF]. Different kind of instabilities have been distinguished, especially stick-slip, variable dynamic friction coefficient, sprag-slip and mode coupling. The mode-coupling mechanism is commonly considered as one of the main cause of squeal noise and many works are based on this assumption. According to this hypothesis, it has been demonstrated that squeal could occur with a constant friction coefficient in relation to the sliding velocity. In this case, the determination of squeal events is generally performed in two steps. The first stage, the Complex Eigenvalues Analysis (CEA), is useful to determine the occurrence of squeal events. It is based on the stability analysis which consists in computing the complex eigenvalues and modes of the linearized problem around the quasi-static steady sliding state. If one of the real parts of the eigenvalues is positive, the associated mode may induce the origin of squeal phenomena. Even if this first step is essential in a design process aimed at clarifying the propensity of squeal events, it is not sufficient to accurately describe the self-sustained vibrations leading to the squeal noise.

Indeed, these vibrations are the results of purely non-linear effects, especially those which may occur in the contact zone, for instance contact/separation (loss of contact) or stick/slip phenomena. Thus, the second stage corresponds to a full non-linear analysis of the problem and aims at computing the vibratory level and spectrum of the self-excited vibrations. More realistic results can be achieved via this non-linear analysis [START_REF] Abubakar | Complex eigenvalue analysis and dynamic transient analysis in predicting disc brake squeal[END_REF] due to the fact that CEA may lead to an underestimation or an over-estimation of the unstable modes observed in the nonlinear regime [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]. In particular, the real parts of the complex modes provide by CEA which traduce their divergence rate in the linear regime do not inform about their amplitude contributions in the nonlinear regime. Calculations of the non-linear vibrations are mostly performed by numerical integrations in the time domain [START_REF] Lorang | Instabilité vibratoire des structures en contact frottant : Application au crissement des freins de TGV[END_REF][START_REF] Lorang | Stability and transient analysis in the modelling of railway disc brake squeal[END_REF][START_REF] Massi | Linear and non-linear numerical approaches[END_REF][START_REF] Loyer | Study of nonlinear behaviors and modal reductions for friction destabilized systems. application to an elastic layer[END_REF][START_REF] Oberst | Nonlinear transient and chaotic interactions in disc brake squeal[END_REF][START_REF] Meziane | Experimental and numerical investigation of friction-induced vibration of a beam-on-beam in contact with friction[END_REF] from given initial conditions closed to the equilibrium or using progressive loadings. In most cases, both transient and steady state vibrations can be identified in the time histories but more complex behaviors may also be observed [START_REF] Loyer | Study of nonlinear behaviors and modal reductions for friction destabilized systems. application to an elastic layer[END_REF][START_REF] Oberst | Nonlinear transient and chaotic interactions in disc brake squeal[END_REF][START_REF] Chevillot | Nonlinear transient vibrations and coexistences of multi instabilities caused by friction in an aircraft braking system[END_REF]].

However, these numerical integrations are often computationally too expensive so that they may be considered impractical in several engineering applications. Assuming the existence of steady state limit cycles, alternative methods have been developed for directly estimating the steady state regime of the non-linear dynamic solutions, like the center manifold, multiple scales, normal forms or harmonic balance methods. The Harmonic Balance Method (HBM) is the most widespread approach. In the specific case of unstable autonomous systems, Coudeyras et al. [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][START_REF] Coudeyras | Periodic and quasi-periodic solutions for multi-instabilities involved in brake squeal[END_REF] notably proposed a new treatment called the Constrained Harmonic Balance Method (CHBM) for computing the periodic or quasi-periodic solutions of non-linear systems prone to mono or multi-instabilities involved in brake squeal. As explained by Nacivet and Sinou [START_REF] Nacivet | Modal amplitude stability analysis and its application to brake squeal[END_REF] the use of this non-linear method in the context of finite element models for engineering applications remains rare due to the requirement of many computational developments. In the mono-instability case, the determination of a steady state solution amounts to finding a periodic solution of the problem which may be viewed as a normal non-linear mode of the problem. Kerschen et al. provide a very detailed description of the concepts of non-linear normal modes and the state-of-the-art in terms of computational implementation [START_REF] Kerschen | Nonlinear normal modes, part i: A useful framework for the structural dynamicist[END_REF][START_REF] Peeters | Nonlinear normal modes, part ii: Toward a practical computation using numerical continuation techniques[END_REF]. A more recent review can be found in [START_REF] Renson | Numerical computation of nonlinear normal modes in mechanical engineering[END_REF]. In order to compute isolated periodic solutions, the most used approaches are the shooting method, the previously mentioned HBM and the orthogonal collocation method. The shooting method is a popular numerical technique [START_REF] Slater | A numerical method for determining nonlinear normal modes[END_REF][START_REF] Ribeiro | Non-linear forced vibrations of thin/thick beams and plates by the finite element and shooting methods[END_REF][START_REF] Seydel | Practical Bifurcation and Stability Analysis[END_REF][START_REF] Stoykov | Numerical computation of periodic responses of nonlinear large-scale systems by shooting method[END_REF][START_REF] Kuether | A numerical approach to directly compute nonlinear normal modes of geometrically nonlinear finite element models[END_REF][START_REF] Wang | Bifurcations of nonlinear normal modes via the configuration domain and the time domain shooting methods[END_REF] which iteratively finds the initial conditions and the period that realize both the periodic motion and a phase condition. It defines as a two-point boundary-value problem using the periodicity condition as explained in the reviews of Kerschen [19] or Nayfeh [START_REF] Nayfeh | Applied Nonlinear Dynamics: Analytical, Computational, and Experimental Methods[END_REF]. It is easy to implement since it only requires the results of some numerical integrations of the problem. However, two main limitations have been identified. The first one is the sensitivity of the local iterative shooting algorithm to its initialization: the method requires a good initial guess in order to converge. The second limitation is linked to the Jacobian computation, which makes the algorithm prohibitively expensive for models with many degrees-of-freedom like large scale finite element models.

In this paper, it is proposed to test the performance of the shooting method in the case of non-linear autonomous systems subjected to flutter instability and especially for non-smooth contact dynamical systems with planar friction prone to friction-induced mono-instability. In order to overcome the first disadvantage of the shooting method (i.e. the sensitivity of the local iterative shooting algorithm), an original approximation is proposed and discussed for the initialization of the iterative algorithm. This approximated initial guess is computed using an approach combining the results of the linear stability analysis for the shapes and the period with a non-linear power balance for the amplitude. In this study, the proposed shooting method is applied to a three degree-of-freedom system sliding on a plane and submitted to planar friction with a constant frictional coefficient (with non-regularized Signorini unilateral contact and Coulomb friction laws). The second disadvantage of the shooting method (i.e. the Jacobian computation in the case of large scale models) is thus not treated and lies out of the scope of the paper. The paper is divided in three parts. In the first one, the mechanical system, the model of frictional contact and the numerical integration scheme are presented. In the second part, the shooting method with appropriate initial guess is detailed. Finally, numerical simulations and results provided by both the proposed shooting approach and a full time integration are compared and discussed.

Modeling and formulation of the problem

Figure 1 shows the frictional mechanical system under consideration. This system has been previously presented in detail in [START_REF] Charroyer | Parametric study of the mode coupling instability for a simple system with planar or rectilinear friction[END_REF]. It is composed of a mass m in frictional contact with a rigid plane moving with a constant rectilinear velocity V .

γ denotes the angle between the direction t of the imposed velocity of the plane and the coordinate x of the mechanical system. The three springs (i.e. k 1 , k 2 and k n ) hold the mass against the moving plane and an external force F is applied to the mass. The set of equations describing the dynamics of the system can be written as:

M Ü + C U + KU = R + F (1)
where

U = U x U y U z T
, U and Ü are respectively the displacement, velocity and acceleration vectors along the three directions x, y and z. Expressions of the mass matrix M and damping matrix C are given by M = mI and C = cI, respectively, where I defines the identity matrix. The stiffness matrix K is given by:

K =      k 1 0 k 1 tan α 0 k 2 k 2 tan β k 1 tan α k 2 tan β k n      (2)
with

k 1 = k 1 cos 2 α, k 2 = k 2 cos 2 β and k n = k n + k 1 tan 2 α + k 2 tan 2 β.
The vector of the external force F is written as:

F = 0 0 -F n T with F n > 0 (3)
Finally, R includes the contact normal reaction and the friction force applied on the center of mass due to frictional contact with the plane. The contact is assumed to be unilateral and the friction is governed by Coulomb's law. The vector of the frictional contact forces can be expressed as:

R = P t T R t + P b T R b + P n T R n (4)
where R t , R b and R n are, respectively, the two components of the friction force and the contact reaction which are applied on the mass, according to the directions t, b and z. The matrices P t , P b and P n are used to transpose forces from the contact to the global coordinate system. They are given by

P b = {-sin γ cos γ 0} (5)
P t = {cos γ sin γ 0} (6)
P n = {0 0 1} (7)
To deal with the unilateral contact, a non-regularized Signorini law is chosen:

         U n ≥ 0 R n ≥ 0 U n R n = 0 (8)
where U n corresponds to the gap between two solids. It is defined as positive when contact is open.

In dynamics, due to this unilateral contact law, non-regular impacts phenomena occur when the normal gap reaches zero (contact) from a positive value (loss of contact or separation). In this case, an impact law on the velocity must be added since the velocity after the impact is undetermined. In this study an inelastic impact law is chosen so that the velocity after the impact is null. This law is implicitly imposed by the numerical integration scheme exposed at the end of the section.

In the context of friction-induced vibrations due to mode coupling, the importance of the normal oscillations may often lead to loss of contacts and impacts at the friction interface. This has been identified as a crucial phenomenon for the saturation of self-sustained vibrations characterized by slip-separation waves, especially in brake squeal problems [START_REF] Lorang | Instabilité vibratoire des structures en contact frottant : Application au crissement des freins de TGV[END_REF][START_REF] Loyer | Study of nonlinear behaviors and modal reductions for friction destabilized systems. application to an elastic layer[END_REF]. Moreover, inelastic impacts may dissipate some energy since the local kinetic energy moves instantaneously from a positive to a null value. This is particularly true in the discrete problem with only one contact node studied in the paper

The planar friction force is given by the Coulomb's law with a contact friction coefficient µ as stated below:

   Ut -V Ub    = 0 ⇒    R t R b    ≤ µR n    Ut -V Ub    = 0 ⇒    R t R b    = -µR n        Ut -V Ub               Ut -V Ub        (9)
where V defines the value of the imposed velocity V . Ut = P t U and Ub = P b U are the mass velocities according to the directions t and b.

Figure 2 shows the unilateral contact (i.e. a non-regularized Signorini law) and Coulomb's friction laws for the system under study.

As explained for instance in [START_REF] Jean | The non-smooth contact dynamics method[END_REF][START_REF] Alart | A mixed formulation for frictional contact problems prone to newton like solution methods[END_REF], a reformulation of the contact and friction laws can be performed in terms of projections on the positive real set proj R + and on the Coulomb cone proj D(µRn) in order to facilitate the numerical implementation of the contact equations. The following relations are then obtained:

R n = proj R + (R n -ρ n U n ) ∀ρ n > 0    R t R b    = proj D(µRn)      R t R b    -ρ t    Ut -V Ub      ∀ρ t > 0 (10)
For given initial conditions, equations of motion (1) to [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], together with frictional contact laws (10) can be solved numerically using an appropriate integration scheme and a non-linear algorithm. In this study, a modified θ-method (with θ = 0.5) that has been previously proposed by Jean [START_REF] Jean | The non-smooth contact dynamics method[END_REF] and Vola et al. [START_REF] Vola | Consistent time discretization for a dynamical frictional contact problem and complementarity techniques[END_REF] is employed. This is an implicit first-order integration scheme which is adapted for dynamic contact problems, especially when loss of contacts and inelastic shocks occur. In addition, at each time iteration, a fixed point algorithm is used in order to verify the frictional contact laws [START_REF] Massi | Linear and non-linear numerical approaches[END_REF]. This modified theta-method has been previously tested for the prediction of self-excited vibration of mechanical systems with loss of contacts and a non-regularized Signorini law (see for example [START_REF] Lorang | Instabilité vibratoire des structures en contact frottant : Application au crissement des freins de TGV[END_REF][START_REF] Lorang | Stability and transient analysis in the modelling of railway disc brake squeal[END_REF]. The detailed and complete formulation of this specific method can be found for instance in [START_REF] Loyer | Study of nonlinear behaviors and modal reductions for friction destabilized systems. application to an elastic layer[END_REF] and the discussion about the dissipation properties of the integration scheme in [START_REF] Acary | Energy conservation and dissipation properties of time-integration methods for nonsmooth elastodynamics with contact[END_REF]. The full time integration method consists in starting the numerical integration with initial conditions close to the quasi-static steady-sliding state and stopping the process when the steady state limit cycle is reached. The aim of the shooting method presented in the next section is to directly compute the limit cycle in the case of mono-instability.

Shooting method with appropriate initial guess

In this section we propose to introduce the shooting method for determining the steady state non-linear responses of self-excited systems subject to friction-induced mono-instability. The shooting method for autonomous system with an unknown period is first presented. Secondly, the approach used to estimate an appropriate initial guess for the shooting algorithm is described and discussed.

Shooting method applied to an autonomous system with mono-instability

The shooting method is an iterative method based on a numerical time integration procedure which determines a periodic solution by iterating on the initial conditions and the period of the vibrational response. In the context of self-excited systems, the fundamental frequency of the non-linear steady-state solution is unknown. Indeed it may differ from the frequency of the unstable mode that may be calculated for instance via a stability analysis. Thus, the implementation of the shooting method requires not only the determination of the periodic solution (i.e. the displacement and velocity vectors of the system) but also the fundamental frequency.

The problem defined in Equation ( 1) is rewritten in state space form (i.e.

Z = U U T , Z ∈ R n ): Ż = G (Z) (11)
where

G (Z) = E -AZ (12)
A =   0 -I M -1 K M -1 C   (13)
E = 0 M -1 (R + F) T (14)
Assuming that a periodic solution of the system (11) exists, this solution satisfies the following conditions:

Z (t + T ) = Z (t) (15)
Z (t + τ ) = Z (t) if 0 < τ < T (16)
with T the period of the autonomous system.

According to the shooting method, the initial value problem given by:

Ż = G (Z, t) Z (t = 0) = Z 0 (17)
can be converted into a boundary value problem defined by:

Ż = G (Z, t) Z (T ) = Z (t = 0) = Z 0 (18)
where Z 0 corresponds to the chosen initial conditions and T defines the period of the system. Thus, the shooting method consists in seeking an initial condition Z 0 and a solution Z (T, Z 0 ) with a minimal unknown period T (i.e. the period value is unknown due to the fact that we consider an autonomous system) such that:

H (Z 0 , T ) = Z (Z 0 , T ) -Z 0 = 0 (19)
This last relation describes the condition of periodicity of the system.

In order to find the initial condition Z 0 and the period T that verify condition of periodicity ( 19), an iterative Newton-Raphson algorithm is used. It consists in applying incremental corrections ∆Z 0 i and ∆T i at the i th iteration verifying:

H i + ∂H ∂Z 0 Z0 i ,T i ∆Z 0 i + ∂H ∂t Z0 i ,T i ∆T i + HOT = 0 (20)
where HOT corresponds to neglected higher order terms. In this equation, H i , Z 0 i and T i denote the values of H, Z 0 and T at the i th iteration. The corrected initial condition and period are then given by:

Z 0 i+1 = Z 0 i + ∆Z 0 i (21)
T i+1 = T i + ∆T i (22)
As shown by Equation [START_REF] Renson | Numerical computation of nonlinear normal modes in mechanical engineering[END_REF], this procedure requires the calculation of the two following partial derivatives:

∂H ∂t Z0 i ,T i = ∂Z(Z 0 i , t) ∂t t=T i = G Z(Z 0 i , T i ) (23)
∂H ∂Z 0 Z0 i ,T i = ∂Z ∂Z 0 Z0 i ,T i -I (24)
The Jacobian matrix ∂Z ∂Z 0 (Z 0 , t) of the solution with respect to the initial conditions is the monodromy matrix. It is a T -periodic n × n matrix. It must be noted that the stability of periodic solutions of the system of non-linear differential equations ( 18) can be easily determined using the Floquet theorem and the monodromy matrix. This requires the calculation of the complex eigenvalues ρ i of this matrix, called Floquet multipliers. Floquet multipliers represent the dimensionless numbers that indicate the period-to-period increase or decrease of the perturbation. If a Floquet multiplier has a modulus greater than one, the periodic orbit is unstable, otherwise if all Floquet multipliers have a modulus less than 1 the periodic orbit is stable [START_REF] Nayfeh | Applied Nonlinear Dynamics: Analytical, Computational, and Experimental Methods[END_REF][START_REF] Hale | Dynamics and Bifurcations[END_REF].

In this study, the monodromy matrix is computed using a finite difference approximation but other techniques could be used [START_REF] Nayfeh | Applied Nonlinear Dynamics: Analytical, Computational, and Experimental Methods[END_REF]. Considering the following initial conditions:

Z(0) + k = Z 0 + δe k (25) Z(0) - k = Z 0 -δe k ( 26
)
where e k represents the k th column of identity matrix, the k th column of the Jacobian matrix is computed by using the relation:

∂Z ∂Z 0 k Z0 i ,T i = Z Z 0 i + δe k , T i + k -Z Z 0 i -δe k , T i - k 2δ (27)
At this point it is necessary to address the problem of the influence of non-smoothness of the system on behavior of infinitesimal perturbation. Indeed, function H is non smooth due to the possible velocity discontinuities in case of impacts.

However, from a theoretical point a view, the function is piecewise smooth (and even piecewise linear if we exclude the non linearity of the friction force direction) and the problem only occurs at singular points when initial conditions lie very close to the discontinuities. In practice, the computation of the Jacobian matrix of the algorithm by finite differences may be sensitive to this non smoothness. In order to avoid crossing the discontinuities, parameter δ in the previous equations must be sufficiently small. An interesting point in that the algorithm seems to converge even when larger values of delta (crossing the discontinuities) are considered. This secant approximation may even accelerate the convergence in some cases. Consequently, the finite difference method may appear as an advantageous method for non smooth systems because the convergence of the algorithm, even if it is not proved for such a a secant approximation of the Jacobian, may be controlled by an appropriate setting of δ.

System (20) has 2 × N equations for 2 × N + 1 unknowns with N the number of degree of freedom. This system is completed by one equation determining the system phase. Several phase conditions can be found in [START_REF] Nayfeh | Applied Nonlinear Dynamics: Analytical, Computational, and Experimental Methods[END_REF]. In this study the condition of orthogonality between the increment of the initial conditions ∆Z 0 and G is chosen:

G (Z(t = T )) ∆Z 0 = 0 (28)
Finally, the following matrix system has to be solved at each iteration i th :

  ∂Z ∂Z0 Z0 i ,T i -I G Z(Z 0 i , T i ) G T Z(Z 0 i , T i ) 0      ∆Z 0 i ∆T i    =    Z 0 i -Z(Z 0 i , T i ) 0    ( 29
)
The procedure is repeated until the two following convergence criteria are satisfied:

H i / Z 0 i+1 < 1 ( 30
)
|∆T i |/|T i+1 | < 2 (31)
As previously explained the example treated in the proposed study concerns a non-smooth contact dynamical system (non-regularized Signorini unilateral contact and Coulomb friction laws) with planar friction. Indeed the implementation of a robust nonlinear technique to obtain the limit cycles is challenging, the main difficulty resulting in the fact that the limit cycles are directly linked to a non-smooth contact configuration (no other nonlinearity is considered to saturate the growing oscillations). So a first-order θ-method time integration scheme is developed and used during the shooting process for the computation of the vibrational nonlinear responses over a period T i . For more details, the interested reader is referred to the following papers [START_REF] Lorang | Instabilité vibratoire des structures en contact frottant : Application au crissement des freins de TGV[END_REF][START_REF] Lorang | Stability and transient analysis in the modelling of railway disc brake squeal[END_REF][START_REF] Loyer | Study of nonlinear behaviors and modal reductions for friction destabilized systems. application to an elastic layer[END_REF]. A fixed point algorithm is also used in order to verify the frictional contact laws at each time integration during the shooting process.

Initialization of the shooting method

One of the most important issues when using the shooting method is to introduce a relevant initial condition in order to enhance the computational efficiency of the proposed shooting algorithm. Considering an energetic point of view, the proposed initialization process can be derived by analyzing the physical self-excited behavior occurring during the development of the steady state vibration. Indeed, for a system subject to mode coupling instabilities, the vibrational amplitudes growth until the occurrence of severe non-linearities (i.e. the generation of contact/separation or slip/stick cycles in our case). These non-linearities induce a saturation of the frictional contact reactions and a decrease of the average injected power ratio due to the contact reactions, compared with the linear case. In the same time, the average power ratio dissipated by viscous damping remains almost constant. This physical phenomenon is responsible for the stabilization of the solution and a convergence towards the steady state non-linear response of the frictional mechanical system.

The proposed approach for the initialization of the shooting method is decomposed into two steps. Firstly a simple linear choice for the approximation of the steady state self-excited solution of the system is proposed. Secondly the calculation of the power balance in the system is developed considering the non-linear saturation of the frictional contact reactions. This leads to the estimation of the initial conditions for the numerical computation of the shooting method.

This proposed process is now developed and discussed in details.

First of all, the solution (that is composed by the displacement U(t) and the velocity U(t)) is approximated by assuming that it is periodic and governed by the unstable mode (that is previously calculated via the stability analysis, see [START_REF] Charroyer | Parametric study of the mode coupling instability for a simple system with planar or rectilinear friction[END_REF] for more details):

U(t) = q Φe i (λ)t (32)
U(t) = q i (λ)Φe i (λ)t ( 33
)
where λ is the eigenvalue of the unstable mode, Φ the associated eigenvector, and q the modal amplitude. (.) and (.)

denote the real and imaginary parts of the variables.

The most relevant value for q (that is denoted q 0 in the following) is searched in order to define the following initial guess:

Z 0 0 =    q 0 (Φ) + U 0 q 0 (λΦ)    ( 34
)
T 0 = 2π/ (λ) (35)
as the initial conditions of the shooting method. U 0 corresponds to the steady sliding equilibrium of the non-linear system (see [START_REF] Charroyer | Parametric study of the mode coupling instability for a simple system with planar or rectilinear friction[END_REF] for more details). It may be noted that Equations [START_REF] Acary | Energy conservation and dissipation properties of time-integration methods for nonsmooth elastodynamics with contact[END_REF] and [START_REF] Hale | Dynamics and Bifurcations[END_REF] correspond to the purely dynamic part of the approximated self-excited vibration of the non-linear system, considering only the participation of one unstable mode.

The formulation given in Equation (34) also assumes that the equilibrium is not drastically affected during the initial transient vibrations of the non-linear system. This assumption, is certainly restrictive and may sometimes be unverified (see [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 example). However, in this first step we are only looking for a simple approximation of the system that allows us to enhance the computational efficiency of the shooting method via the introduction of a relevant initial condition. The final non-linear sought solution is determined by the shooting method.

The second step of the proposed process concerns the calculation of the value q 0 . To achieve this, we propose to balance the average power on one period T of the steady state. The instantaneous power balance of the system can be written in the form:

∂E ∂t + P C = P R + P F (36)
where ∂E ∂t defines the variation of the mechanical energy, P C the power dissipated by damping, P R the power injected or dissipated by the contact reactions and P F the power injected by the external force. Considering a steady state regime, the variation of the mechanical energy and the power due to static forces like the external force are null over one period T . Considering only dynamic forces, averaging over one period and normalizing the powers by the average mechanical energy Ē, the power balance may be expressed in function of the average power ratio τC dissipated by damping and the average power ratio τR injected by the frictional contact reactions:

τC = τR (37) with τC = PC Ē (38) τR = PR Ē (39) PC = 1 T T 0 UT C Udt = 1 2 q 2 (λ) 2 Φ * T CΦ (40) PR = 1 T T 0 UT (R -R 0 )dt (41) Ē = 1 T T 0 ( 1 2 UT M U + 1 2 U T KU)dt = 1 4 q 2 (λ) 2 Φ * T MΦ + 1 4 q 2 Φ * T KΦ (42)
where R 0 denotes the frictional contact reaction of the quasi-static equilibrium U 0 (see [START_REF] Charroyer | Parametric study of the mode coupling instability for a simple system with planar or rectilinear friction[END_REF] for more details) and R defines the vector of the frictional contact forces (expression given in Equation ( 4)). The subscript -denotes an average quantity over the period. As expected, using a linear approximation for the displacements and velocities, the average power ratio dissipated by damping is independent of the modal amplitude q:

τC = 2 (λ) 2 Φ * T CΦ (λ) 2 Φ * T MΦ + Φ * T KΦ (43)
The average power ratio due to frictional contact efforts τR depends on the modal amplitude q and is calculated numerically by saturating the reactions corresponding to the linearized sliding contact problem. Assuming a sliding contact, the dynamic part of the contact reaction is linearly linked to the displacements via P n KU. An approximation of the saturated frictional contact reactions is then given by :

R n = proj R+ (R n0 + (P n KU))    R t R b    = -µR n        Ut -V Ub               Ut -V Ub        if    V -Ut -Ub    = 0, 0 otherwise (44)
We recall that the vector of the frictional contact forces is given by R

= P t T R t + P b T R b + P n T R n .
Finally, the value of q 0 verifying relation (37) is calculated numerically and the shooting algorithm is initialized using Equations ( 34) and (35). The global approach is summarized on Figure 3.

Figure 3: Overview of the shooting approach with appropriate initial guess

Numerical results

The numerical results are decomposed into two main parts: the stability analysis and the non-linear analysis for the computation of the steady state non-linear vibration. The numerical simulations and the results regarding the stability analysis are briefly presented in this paper. A comprehensive study has been previously performed by the authors on this subject in [START_REF] Charroyer | Parametric study of the mode coupling instability for a simple system with planar or rectilinear friction[END_REF]. In the second part of the section, the calculation of the transient non-linear self-excited vibrations is presented. More specifically, the efficiency of the shooting method is illustrated. All the proposed results are analyzed and compared to a reference calculation based on the full time integration scheme as presented in the first section and details for instance in [START_REF] Loyer | Study of nonlinear behaviors and modal reductions for friction destabilized systems. application to an elastic layer[END_REF].

Regarding the shooting method, both values 1 and 2 of the convergence criteria (see Equations ( 31) and ( 31)) have been chosen equal to 10 -3 . The numerical studies are performed with the following parameters: m = 0.01 kg, f 1 = 100

Hz, f 2 = 75 Hz, k i = m × f i N/m, k n = 20 N/m, η = 0.02, c = m × ηf 1 , α = 30 • , β = 60 • and F n = 10 N.

Stability analysis

In order to predict the squeal propensity of the mechanical system, its stability is investigated. It consists in analyzing the behavior of small perturbations around the quasi-static sliding equilibrium. A bilateral contact between the mass and the plane is assumed. The frictional reactions are linearized around the static equilibrium and the complex eigenvalues of the system are calculated. If the real part of all eigenvalues are negative, the system is stable. On the other hand, if at least one of the eigenvalues has a positive real part, then the system is unstable and prone to self-excited vibrations. All the details of this analysis can be found in [START_REF] Charroyer | Parametric study of the mode coupling instability for a simple system with planar or rectilinear friction[END_REF].

Figure 4 shows the stability areas for two parameters (i.e. the friction coefficient µ and the sliding direction γ) in the case of two specific configurations. In the first case, a sliding velocity V = 3 m.s -1 is considered whereas in the second case, V = 0.75 m.s -1 . This choice of configurations are explained in the next section. For each configuration the main unstable area is in the bottom right quarter of Figures 4 and corresponds to directions of the sliding velocity between -90 • and 0 • .

For details on the influence of various parameters of the mechanical system under study and more explanations regarding the stability analysis, the reader may refer to [START_REF] Charroyer | Parametric study of the mode coupling instability for a simple system with planar or rectilinear friction[END_REF].

-

Non-linear vibration

Based on the previous results concerning the stability analysis of the mechanical system under study, the three following cases are chosen for the non-linear simulation:

• case 1: (µ, γ, V ) equals to (0.2244, -30 • ,3 m.s -1 ) -marked with a square on Figure 4(a),

• case 2: (µ, γ, V ) equals to (0.75, -30 • ,0.75 m.s -1 ) -marked with a circle on Figure 4(b),

• case 3: (µ, γ, V ) equals to (0.40, -60 • ,3 m.s -1 ) -marked with a diamond on Figure 4(a),

These three reference cases are used to demonstrate the efficiency and robustness of the proposed shooting method to predict the estimate self-excited vibrations but also to discuss on the limitation and the precautions of use for the proposed shooting method (more specifically in case 3). The cases 1 and 3 are chosen because they lead to self-excited vibrations with contact/separation phenomena while the second case illustrates the occurrence of stick-slip vibrations. When the system is in an unstable area, a perturbation of the static equilibrium generate self-excited vibrations. The vibrational amplitudes increase linearly until occurrence of severe non-linearities (ie. contact/separation phenomenon for the first and third cases and stick-slip phenomenon for the second case) and the establishment of steady state non-linear periodic motions.

For each case, results provided by the proposed shooting method are compared with those obtained from a full time integration. All the results are presented by plotting the phase diagrams of the solutions. The initial guess provided by the approach explained in Section (3.2) is indicated by a blue star on these phase diagrams. The red curves correspond to the results of the direct integration scheme whereas the blue one corresponds to the results of the proposed shooting method.

All the integrations have been performed using 1024 time steps per period which leads cases 1 and 2 to time steps of about 10 µs. This value is the result of a convergence study. The differences with a smaller time step (2048 per period) are about 0.99991 in terms of mean MAC (Modal Assurance Criterion) and less than 0.25 % in terms of mean mechanical energy on the period.

Figure 5 and 6 illustrate the phase diagrams according to the three directions (i.e. x, y and z) for the first and second cases, respectively. It is clearly shown that the results based on the proposed shooting method are in perfect agreement with those obtained from the direct time integration. Moreover, the starting point used by the shooting method is very close to the final estimated limit cycle. This demonstrates the relevance of the use of the power balance criteria proposed in Section 3.2 to be used as a first indicator of the limit cycle amplitudes and more particularly to initiate the calculation of the shooting method. It is worth noting that the non-use of this starting point can lead to difficulties in estimating the limit cycles via the shooting method. This illustrates the efficiency of the proposed non-linear method. It is also observed that the limit cycles obtained may be very complex, with pronounced non-linear behavior such as, for instance, the shock discontinuities of the limit cycles on Figure 5(a) and 5(b). Stick-slip phenomena are also reproduced with a good precision for the second case (see Figures 6(a) and (b)). For cases 1 and 2, the modulus of the eigenvalues of the monodromy matrices are lower than 1 which means that the steady state self-excited responses of the mechanical system are stable.

Considering case 3, Figures 7 show the limit cycles provided by the shooting method compared with the phase diagrams obtained with the direct time integration. Time evolutions of the vertical displacement and velocity given by the direct integration are also shown on Figures 8. Two stages can be distinguished. In a first stage (before about t = 0.15 s), the solution diverges from the equilibrium and seems to stabilize close to the periodic limit cycle given by the shooting method. In a second stage (from about t = 0.15 s), increasing oscillations around this limit cycle are observed with a definitive stabilization on a new limit cycle with about twice the period of the first cycle. This can be explained by the fact that the calculated limit cycle via the shooting method is unstable, which is confirmed by the study of Floquet multipliers of the monodromy matrix provided by the shooting method. In addition, a bifurcation phenomenon of period doubling

-3 -2 -1 0 1 2 3 4 X-Displacement (m) ×10 -3 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 X-Velocity (m.s -1 ) (a)
Figure 5: Phase diagrams provided by the direct time integration procedure (red) and limit cycles given by the shooting method (blue) in the direction x (a), y (b) and z (c) for case 1 can be detected. Indeed, by varying the bifurcation parameter γ (i.e. the direction of the sliding velocity) from -90 • to 0 • , one characteristic Floquet multiplier crosses the unit circle through -1 for γ equal to -70 • and -48 • , which traduces a period-doubling bifurcation of the limit cycle for γ = [-70 • ; -48 • ] as stated for instance in [START_REF] Nayfeh | Applied Nonlinear Dynamics: Analytical, Computational, and Experimental Methods[END_REF]. A period doubling bifurcation in a dynamical system is a bifurcation in which the system switches to a new behavior with twice the period of the original system. In order to find this new limit cycle, the shooting method is used but considering an initial period

T 0 = 4π (λ)
with λ the eigenvalue of the unstable mode. This leads to a stable limit cycle (according to the monodromy matrix) with phase diagrams presented in Figures 9. This new limit cycle is in perfect agreement with the reference limit cycle performed via the direct temporal integration. All the results can lead to several comments. Firstly the shooting method makes it possible to find an unstable limit cycle which is not possible using the full time integration. The shooting

Figure 6: Phase diagrams provided by the direct time integration procedure (red) and limit cycles given by the shooting method (blue) in the direction x (a) and y (b) for case 2 method will not necessarily give the limit cycle obtained by the temporal integration if the input period that is initially chosen for the computation is not close to the period of the stable limit cycle. It generally finds the limit cycle (that may be stable or unstable) that is close to the provided initial conditions whereas the direct integration scheme leads to the stable limit cycle only. However, the common use of Floquet theory which test the stability of a limit cycle solution can

give information concerning the new behavior of the original system.

Another interesting point is the study of the limit cycle periods. For each case, it can be noted that the fundamental frequency of the steady state self-excited oscillations (calculated for both the direct integration time and the shooting method) is slightly different from the frequency of the unstable mode provided by the stability analysis as indicated in Table 1. These results highlight the need to consider non-linear analysis not only to estimate the self-excited vibrations but also to accurately predict the value of the fundamental frequency of the steady state self-vibration during squeal events.

This also demonstrates the need to consider the frequency of the self-excited vibration of the system as an unknown. We recall that the initial input condition for the frequency in the shooting process corresponds to the frequency of the unstable mode.

Stability analysis Direct time integration Shooting method

Case 1: (µ, γ, V )=(0.2244, -30

Figure 7: Phase diagrams provided by the direct time integration (red) and unstable limit cycles given by the shooting method in the direction x (a), y (b) and z (c) for case 3

Finally, the convergence of the proposed shooting method is discussed. In our present study, the calculation of the steady state non-linear response is performed by considering the following values for the convergence criterion: 1 = 10 -3

and 2 = 10 -3 . Only 3 iterations are needed for the convergence of the shooting method in the three cases, excepted when the limit cycle with the double period is searched. In this specific case, the initial guess is less adapted and the algorithm needs 7 iterations in order to converge.

Finally Table 2 gives a comparison of the periods integrated by using the shooting method and the direct integration to obtained the self-excited steady-state vibration. To provide a valuable comparison between the shooting method and the direct integration, the same initial conditions (i.e. the initial guess derived in Section 3.2) are used. It appears that the number of needed periods for the calculation of the results is lower with the proposed shooting method than the direct Periods integrated with the shooting method 39 39

Periods integrated with the direct time integration ≈ 60 2

Table 2: Number of iterations for the convergence of the proposed shooting method and comparison of the number of integrated periods for each method time integration for Case 1. However, for Case 2, the results are calculated faster by using the time integration. This may be explained by the fact that the time interval of transient oscillations before the appearance of the steady-state vibrations is very short for Case 2. For Case 1, the transient phase before the appearance of self-excited vibrations is still important because of the contact/separation phenomena. Moreover, it may be recalled that, in contrast to the shooting method, one of the major disadvantages of the direct time integration method is that it is not known a priori when the limit cycles are reached and therefore when the process of the modified θ-method integration scheme can be stopped. Another advantage of the shooting method is that at each shooting step, the integrations which are necessary to determine the Jacobian can be computed in parallel. This is not the case for the direct method.

Parametric study with the direction of the sliding velocity γ

In the previous section, we have noted that some configurations of the system require the search for a stable limit cycle with a period doubled. A parametric study is now presented to determine when the doubling period bifurcation occurs. In [START_REF] Charroyer | Parametric study of the mode coupling instability for a simple system with planar or rectilinear friction[END_REF], it has been highlighted that one of the most important parameters determining the stability of a system subjected to

Figure 9: Phase diagrams provided by the direct time integration (red) and stable limit cycles given by the shooting method in the direction x (a), y (b) and z (c) for case 3 planar friction is the direction of the sliding velocity. This parameter is then chosen in our study. The value of the sliding velocity is equal to 0.75 m.s -1 and the friction coefficient µ = 0.4. The direction of the sliding velocity γ varies from -10 • to -80 • , as illustrated by the red line on the Figure 10.

For the starting point (γ = -10 • ), two calculations are performed by applying the shooting method initialized with the presented previously energy criterion. The first one allows one to determine the T-periodic limit cycles and the second one the potential 2T-periodic limit cycles. Then a sequential continuation strategy is determined. The direction of the sliding velocity is decreased by For a direction of the sliding velocity between -68 • and -50 • , it is noticed that there is a real eigenvalue moving out of the unit circle in the negative direction (cross by -1, as indicated in Figure 12) with the remaining eigenvalues staying inside the unit circle. So the T-periodic solution is unstable on this angle range (i.e. one characteristic multiplier satisfy |ρ| > 1 for γ = [-68 • ; -50 • ]) and the period-doubling bifurcation occurs (i.e. the system switches to a new behavior with twice the period of the original system). By calculating eigenvalues of the monodromy matrix for the 2T-periodic solution for γ = [-68 • ; -50 • ], it is then found that all eigenvalues are staying inside the unit circle, and this leads to stable 2T-periodic solution. This phenomenon explains the difference between the curves previously shown in Figure 11 for the T-periodic and 2T-periodic solutions. For the reader comprehension, the maximal amplitudes which correspond to the stable solutions are marked by a black star on Figure 11.

Conclusion

A shooting method with a specific initial guess based on an energetic criterion is proposed in order to calculate the non-linear self-excited vibration of mechanical systems subject to friction-induced vibration. The effectiveness of the proposed non-linear method is tested on a non-smooth contact dynamical system with planar friction. For such systems, a specific implicit non dissipative integration scheme has to be used due to the fact that the nonlinear self-excited vibrations are directly linked to the non-smooth contact configuration. So the use of a first-order θ-method time integration scheme in connection with the shooting method is essential in order to accuracy estimate the saturation of the growing oscillations are generated when loss of contacts and inelastic shocks occur. It is shown that the results provided by the shooting method correlate well with those obtained with a full time integration. It demonstrates the potential of the shooting method with a specific initial guess to predict the self-excited vibrations of non-smooth contact dynamical systems. This study illustrates also the fact that unstable steady state non-linear solutions of the nonsmooth contact dynamical system with planar friction and specific bifurcations like the period-doubling phenomenon can be successfully detected by analyzing the complex eigenvalues of the monodromy matrix which is calculated during the shooting process.

Finally it is illustrated that adhesion or separation events occurs (i.e. the generation of contact-separation or slip-stick limit cycles) for a non-smooth contact dynamical system with planar friction by considering a Coulomb's law with a contact friction coefficient.

Figure 1 :

 1 Figure 1: Description of the mechanical system

Figure 2 :

 2 Figure 2: Unilateral contact and Coulomb's friction laws

Figure 4 :

 4 Figure 4: Stability area for (a) V = 3 m.s -1 (b) V = 0.75 m.s -1

Figure 8 :

 8 Figure 8: Time evolutions of (a) displacement and (b) velocity provided by the direct time integration in the z direction for case 3

Figure 10 :

 10 Figure 10: Illustration of the range of the parametric study

Figure 11 illustrates

 11 Figure 11 illustrates the evolution of the maximal amplitude of the displacement with the direction of the sliding velocity. The blue and red curves represent the T-periodic and 2T-periodic limit cycles which differ for an angle between -68 • and -50 • . This means that we have the coexistence of the two limit cycles in this study area. For an angle between [-80 • ; -70 • ] and [-48 • ; -10 • ], it can be noted that the shooting method for the calculation of the potential 2T-periodic solution finds the T-periodic solution because of the non-existence of the 2T-periodic solution. To determine the stability of the T-periodic and 2T-periodic solutions, the calculation of the eigenvalues of the monodromy matrix is carried out.

Figure 12 illustrates

 12 Figure 12 illustrates evolutions of eigenvalues for the T-periodic solutions. It is noted that the T-periodic limit cycles are stable for an angle between [-80 • ;-70 • ] and [-48 • and -10 • ] (i.e. according to the Floquet's theory the solution is stable due to the fact that all the characteristic multipliers satisfy |ρ j | < 1 or all the eigenvalues are staying inside the unit circle).

3 Figure 11 :λFigure 12 :

 31112 Figure 11: Evolution of the maximal amplitude of the displacement with the direction of the sliding velocity γ -blue line = T-periodic solution ; red line = 2T-periodic solution ; black star = stable solution

Table 1 :

 1 Frequency of the self-excited vibration provided by the stability analysis (from the eigenvalue of the unstable mode), the direct time integration and the shooting method

	• ,3 m.s -1 )	91.88 Hz	91.96 Hz	91.95 Hz
	Case 2: (µ, γ, V )=(0.75, -30 • ,0.75 m.s -1 )	93.40 Hz	90.84 Hz	90.79 Hz
	Case 3: (µ, γ, V )=(0.40, -60 • ,3 m.s -1 )	94.94 Hz	46.38 Hz	46.93 Hz

Acknowledgement

The authors would like to thank the financial support provided by Celya (Centre Lyonnais d'Acoustique, Université de Lyon).

J.-J. Sinou acknowledges the support of the Institut Universitaire de France.