L Charroyer

O Chiello

J - J Sinou

Estimation of self-sustained vibration for a finite element brake model based on the shooting method with a reduced basis approximation of initial conditions

The objective of this paper is to discuss and propose an original nonlinear method for the estimation of nonlinear vibrations for mechanical systems subject to friction-induced noise and vibration. To fulfill such an objective the computation of nonlinear dynamic steady-state solutions of autonomous non-smooth contact systems prone to mono-instability is performed by developing an extension of the shooting method. This adaptation consists in using a reduced basis in the iterative process of this nonlinear method in order to seek the unknown initial conditions that verify condition of periodicity of the nonlinear solution of the problem. Efficiency of the proposed approach is illustrated through numerical examples for the prediction of self-sustained vibrations of a railway braking system.

Introduction

This work is part of research intended to mitigate or eliminate the squeal noise produced by railway vehicles during braking [START_REF] Berndt | Experimental and theoretical investigation of brake squeal with disc brakes installed in rail vehicles[END_REF][START_REF] Lorang | Tgv disc brake squeal[END_REF][START_REF] Lorang | Stability and transient analysis in the modeling of railway disc brake squeal[END_REF][START_REF] Brizard | Performances of some reduced bases for the stability analysis of a disc/pads system in sliding contact[END_REF][START_REF] Sinou | A global strategy based on experiments and simulations for squeal prediction on industrial railway brakes[END_REF][START_REF] Sinou | Non smooth contact dynamics approach for mechanical systems subjected to friction-induced vibration[END_REF]. Indeed the squealing of disc brakes installed on rail vehicles is a very unpleasant noise for passengers and residents in the station due to the fact that brake noise can reach very high sound levels up to 105 dB(A).

The proposed study is focused on the prediction of the occurrence of squeal noise and the estimation of self-sustained vibration for mechanical systems subjected to friction-induced vibration. In recent years extensive reviews have been carried out to summarize the state of knowledge of experimental and theoretical investigations for the problem of friction-induced vibration in brake systems [START_REF] Crolla | Brake noise and vibration: state of art[END_REF][START_REF] Ouyang | Numerical analysis of automotive disc brake squeal : a review[END_REF][START_REF] Ibrahim | Friction-induced vibration, chatter, squeal and chaos part ii : Dynamics and modelling[END_REF][START_REF] Kinkaid | Automotive disc brake squeal[END_REF]. It is now recognized and admitted that brake squeal is a complex phenomenon. Since the beginning of the 21th century various types of model and different mechanism of squeal generation have been built from a theoretical point of view [START_REF] Tonazzi | Instability scenarios between elastic media under frictional contact[END_REF][START_REF] Tison | Improvement in the predictivity of squeal simulations: Uncertainty and robustness[END_REF][START_REF] Nobari | Uncertainty quantification of squeal instability via surrogate modelling[END_REF][START_REF] Lü | A unified approach for squeal instability analysis of disc brakes with two types of random-fuzzy uncertainties[END_REF][START_REF] Treimer | Uncertainty quantification applied to the mode coupling phenomenon[END_REF][START_REF] Stender | Impact of an irregular friction formulation on dynamics of a minimal model for brake squeal[END_REF]. The purpose of such numerical studies is focused on a better understanding of the generation of friction-induced vibration as well as the prediction of its frequencies for different mastered or not operating conditions. However, the prediction of squeal frequencies and nonlinear friction-induced vibrations in real industrial applications remains difficult.

The goal of the present study is to provide a numerical strategy for the prediction of brake squeal and more specifically for the determination of the self-sustained nonlinear vibration. For such purpose the study focuses on the modeling of self-sustained structural vibrations in presence of non-smooth frictional contact and the understanding of the squeal generation mechanism for a railway disc brake. It is assumed that 1 L. Charroyer, O. Chiello and J-J. Sinou, Estimation of self-sustained vibration for a finite element brake model based on the shooting method with a reduced basis approximation of initial conditions, Journal of Sound and Vibration, 468, 115050, 2020. doi.org/10.1016/j.jsv. 2019.115050 vibration is due to structural mode coupling [START_REF] Hoffmann | A minimal models for studying properties of the mode-coupling type instability in friction induced oscillations[END_REF][START_REF] Charroyer | Parametric study of the mode coupling instability for a simple system with planar or rectilinear friction[END_REF]. This coupling leads to instability of the sliding quasistatic equilibrium of the system and the development of self-sustained vibrations. Two steps are necessary to calculate these vibrations. Firstly, the occurrence of the vibrations is determined through the stability analysis of the quasi-static equilibrium. Secondly, a nonlinear analysis is developed to calculate the vibration amplitudes and the spectral content of the self-sustained vibrations. It is well known that this second step is very challenging and can be very complex to implement in the case of large finite element models. This nonlinear analysis is generally performed by using a numerical time integration scheme from given initial conditions that are close to the equilibrium. Taking a look at the vibration time histories obtained by this technique, two different states may be generally distinguished. In the transient state, nonlinear forces appear and gradually stabilize the unstable solution whereas in the steady state, vibrations are stabilized with potential strongly nonlinear events at the friction interface such as loss of contact, stick or impacts [START_REF] Charroyer | Self-excited vibrations of a non-smooth contact dynamical system with planar friction based on the shooting method[END_REF].

However, direct time integration is not very appropriate in the case of large numerical models with many degrees of freedom. Indeed, the computational cost due to calculation of the transient state is very high although the determination of the steady state is generally sufficient from a practical point of view. Consequently, an interesting numerical strategy lies in the direct estimation of the steady state nonlinear solution. Assuming the existence of steady state limit cycles, alternative methods have been developed like widespread harmonic balance methods (HBM). For unstable autonomous systems, the Constrained Harmonic Balance Method (CHBM), proposed by 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], is dedicated to the computation of periodic or quasi-periodic solutions of nonlinear systems prone to mono or multi-instabilities involved in brake squeal. However, its application to industrial engineering applications remains rare due to the requirement of many computational developments [START_REF] Nacivet | Modal amplitude stability analysis and its application to brake squeal[END_REF]. 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 nonlinear mode of the problem [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][START_REF] Renson | Numerical computation of nonlinear normal modes in mechanical engineering[END_REF]. In order to compute these 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] Stoykov | Numerical computation of periodic responses of nonlinear large-scale systems by shooting method[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] Slater | A numerical method for determining nonlinear normal modes[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 is easy to implement since it only requires the results of some numerical integrations of the problem.

The present study proposed a new method for estimating the steady state regime of nonlinear dynamic solutions in case of mono-instability. This numerical method is based on the classical shooting technique combined with a specific initial guess determined using an energetic criterion as proposed in the work of Charroyer et al. [START_REF] Charroyer | Self-excited vibrations of a non-smooth contact dynamical system with planar friction based on the shooting method[END_REF]. In addition, an enhancement of the method is proposed in the case of large finite element models. In order to make this nonlinear technique more efficient in terms of computational cost, an original approximation of its iteration process is undertaken by the use of a reduced basis for the initial conditions.

The paper is organized as follows: firstly, a description of the finite element model of railway disc brake with non-smooth contact formulation is briefly presented. Secondly, the proposed alternative version of the shooting method is discussed. Finally numerical results are investigated and the estimations of the steady states via the proposed strategy are compared with those obtained by a direct time integration.

Formulation of the problem

Fig. 1 illustrates the brake system model under study. The simple finite element model is derived from a more complex one, which was developed to find instabilities in a TGV (French high-speed train) braking system [START_REF] Sinou | A global strategy based on experiments and simulations for squeal prediction on industrial railway brakes[END_REF]. This brake system is composed of the two main components, which contribute to the brake squeal phenomenon: one thick annular disc and a set of 6 small cylindrical pins applied on either sides of the disc. The flexibility of the pads' supporting structures (i.e. the backplates and its supports) is modeled with spring elements. The upper edge of the hub/disc assembly is clamped. The full finite element model has 31974 degrees of freedom (27090 degrees-of-freedom for the hub/disc assembly and 4884 degrees-of-freedom for all the lining). The nonlinear interface that corresponds to the frictional contact zone between the disc and the 12 small cylindrical pins includes 228 pairs of nodes, which represents 1368 degrees of freedom. All the material and geometrical properties of the brake system are given in Tab. 1. A Rayleigh viscous damping model is used with different parameters for each substructure (i.e. the hub-disc assembly, the pads and the supporting structures) in order to take into account the dissipation in this simple brake system. The general form is given by C = αM + βK where α and β are the Rayleigh damping coefficients given in Tab. 1. The choice of material parameters (elasticity, Rayleigh damping) and operational conditions (boundary conditions, rotation speed, braking force) is realistic and has been carried out based on various measurements and FE model updating made on a TGV brake system [START_REF] Sinou | A global strategy based on experiments and simulations for squeal prediction on industrial railway brakes[END_REF]. In particular, the characterization of the dynamics of the separate subsystems (disc, pads and supporting backplates) has been performed by experimental modal analysis.

Finally the equations of motion of the brake system can be written in a discrete form as follows:

M Ü + C U + KU = F + R = F + P n T R n + P t T R t + P b T R b (1)
where M, C and K are mass, damping and stiffness matrices. F corresponds to the external force. R n , R t and R b define respectively the contact normal reaction and the two tangential friction forces (along the directions t and b at a particular contact point) applied on each cylindrical pin by the disc, whereas P n , P t and P b are the projection matrices on the relative displacements between the disc and the pads along the normal and tangential directions defined at each contact node, respectively. t defines the principal tangential direction given by the sliding velocity and b denotes the second tangential/radial direction such that b = n ∧ t, as illustrated in Fig. 2. A Coulomb law with a constant friction coefficient µ is used and a non-regularized Signorini law is chosen in order to deal with unilateral contact on the interface (see [START_REF] Sinou | A global strategy based on experiments and simulations for squeal prediction on industrial railway brakes[END_REF][START_REF] Sinou | Non smooth contact dynamics approach for mechanical systems subjected to friction-induced vibration[END_REF][START_REF] Loyer | Study of nonlinear behaviors and modal reductions for friction destabilized systems. application to an elastic layer[END_REF] for more details). An usual reformulation of the contact and friction laws in terms of projections on the positive real set proj R + and on the Coulomb cone proj D(µRn) at each mesh node is used in order to facilitate the numerical implementation of the contact equations. Expressions of the contact normal reaction R n and the two tangential friction forces R t and R b are given by [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]:

R n = proj R + (R n -ρ n G n ) (2)
R t R b = proj D(µRn) R t -ρ t ( Ġt -V) R b -ρ t Ġb (3)
where V defines the vector of the amplitude of the imposed velocity at each interface node. (G n , G t , G b ) correspond to the relative displacements in the local coordinate system associated at each contact node (i.e. (G n , G t , G b ) = (P n , P t , P b )U). ρ n (ρ t , respectively) is a diagonal matrix that is composed of positive real numbers called normal augmentation parameters (tangential augmentation parameters, respectively). The over-dot denotes the time derivative. It can be noted that contact/separation and sticking/sliding events at each node of the frictional interface of cylindrical pins can be generated during transient simulations.

3 Nonlinear simulation

Preamble

The global strategy for the prediction of friction-induced vibrations is generally decomposed into two main parts: the stability analysis for the propensity of brake squeal and the nonlinear analysis for the computation The stability analysis is classically conducted by considering small perturbations about the sliding equilibrium point (see [START_REF] Sinou | A global strategy based on experiments and simulations for squeal prediction on industrial railway brakes[END_REF][START_REF] Sinou | Non smooth contact dynamics approach for mechanical systems subjected to friction-induced vibration[END_REF] for more details). This consists in linearizing the equations governing the nonlinear dynamical behavior around this sliding equilibrium point by considering that the small imposed disturbance does not change the contact state at each frictional interface. Complex eigenvalues and complex modes associated with this linearized problem are calculated to determine if the equilibrium point is stable or unstable. If at least one eigenvalue has a positive real part, the equilibrium sliding point is unstable, which implies that the small imposed perturbation of equilibrium leads to increase in vibrations. This growth rate of linear vibration is governed by one or more unstable modes and its validity is restricted by the appearances of new nonlinear events due to contact and no-contact events at each friction interface.

So, in order to predict the self-sustained vibrations, a numerical resolution of the complete nonlinear system has to be performed in addition to the stability analysis by implementing an appropriate integration scheme and a nonlinear algorithm. For the prediction of nonlinear self-sustained vibration of mechanical systems with loss of contacts and a non-regularized Signorini law, the modified θ-method is classically used with a value θ = 0.5 to avoid numerical damping.

More details and the complete formulation of the modified θ-method can be found for instance in [START_REF] Sinou | Non smooth contact dynamics approach for mechanical systems subjected to friction-induced vibration[END_REF][START_REF] Loyer | Study of nonlinear behaviors and modal reductions for friction destabilized systems. application to an elastic layer[END_REF]. Even if the time integration is effective in determining nonlinear self-excited vibratory levels, it is generally time consuming. Therefore it can be considered to use nonlinear method for the prediction of the steady state nonlinear responses of self-excited systems. In the following section, we propose to develop the shooting method for determining limit cycle in the case of mono-instability. This nonlinear method has the advantage to directly compute the limit cycle without calculating transient squeal events. It has been previously proposed by Charroyer et al. [START_REF] Charroyer | Self-excited vibrations of a non-smooth contact dynamical system with planar friction based on the shooting method[END_REF] in the case of an academic mechanical system with 3 degrees-of-freedom.

In the following, the use of the shooting method with a specific initial guess will be described. More specifically, the main originality of the proposed shooting method lies in the use of a reduced basis for the iteration process of this nonlinear method while the resolution of the nonlinear problem is still performed on the complete basis of the mechanical problem. This new formulation and adaptation of the shooting method will be discussed in Section 3.3.

Shooting method

In this section, the shooting technique coupled with the modified θ-method is briefly presented in order to predict the nonlinear periodic solution for an autonomous system with mono-instability. For more details, the interested reader can refer to the work of Charroyer et al. [START_REF] Charroyer | Self-excited vibrations of a non-smooth contact dynamical system with planar friction based on the shooting method[END_REF].

Eq. ( 1) can be written in the following form:

Ż = g (Z) with Z = U U T and g (Z) = E -AZ (4)
where

A = 0 -I M -1 K M -1 C and E = 0 M -1 (R + F) T (5)
The solution of the nonlinear system is denoted Z (Z 0 , t) where Z 0 corresponds to the initial conditions such that Z (Z 0 , t = 0) = Z 0 . The objective of the shooting method is to find out an initial condition Z 0 and the minimal period T in such a way that Z (Z 0 , t) is periodic. This can result in an iterative process for minimizing the following residual function:

H (Z 0 , T ) = Z (Z 0 , T ) -Z 0 (6)
where Z (Z 0 , T ) is calculated by the modified θ-method on a discrete time basis of the interval [0; T ]. In order to find the initial condition Z 0 and the period T such that H (Z 0 , T ) = 0, incremental corrections ∆Z 0 i and ∆T i are applied at the i th iteration verifying:

H i + ∂H ∂Z 0 Z 0 i ,T i ∆Z 0 i + ∂H ∂t Z 0 i ,T i ∆T i + HOT = 0 (7)
where HOT corresponds to neglected higher order terms. H i , Z 0 i and T i correspond to the values of H, Z 0 and T at the i th iteration. The calculation of the two partial derivatives defined in Eq. ( 7) are given by:

∂H ∂t Z 0 i ,T i = ∂Z(Z 0 i , t) ∂t t=T i = g Z(Z 0 i , T i ) (8) ∂H ∂Z 0 Z 0 i ,T i = ∂Z ∂Z 0 Z 0 i ,T i -I 2N (9)
In this last expression, I 2N is the 2N × 2N identity matrix where N is the total number of degrees of freedom of the system and ∂Z ∂Z 0 is a 2N × 2N matrix representing the variation of the solution with the initial conditions and is called monodromy matrix. It must be computed at each iteration.

This iterative process is repeated until two convergence criteria are satisfied:

H i / Z 0 i+1 < 1 ( 10
)
|∆T i |/|T i+1 | < 2 (11)
where 1 and 2 are chosen residuals.

Adaptation of the shooting method based on a reduction basis for the initial conditions

One of the main difficulties of applying the shooting method to a finite element model is the number of degrees of freedom of the model, which can be very important. The main problem lies in the size of the vector of initial conditions and its use in the iterative process (i.e. the size of the monodromy matrix) to find initial conditions and a period that ensure a periodic movement of the autonomous system. In the following a new approach is developed, thus making it possible to reduce the number of degrees of freedom on which the process of the shooting method iterates.

We propose to seek the initial conditions in a reduced phase space of the form:

Z 0 = U 0 U0 = U stat 0 + T 0 0 0 T 0 Q 0 ( 12
)
where T 0 is a reduced basis of size m of the displacement field. Q 0 the vector of amplitudes for the reduced displacements and velocities. U stat defines a static field vector. The idea is to iterate on the reduced vector Q 0 of size m instead of the full vector of size N (with N >> m). The choice of the reduction basis is crucial for the success of the method. In this paper, a very simple base including only the static deformation and the unstable mode provided by the stability analysis is proposed (see section 5.2) but other more complete bases can be tested including stability modes or component modes in the whole frequency range. Some bases adapted to the friction-induced vibration problems are notably proposed in [START_REF] Brizard | Performances of some reduced bases for the stability analysis of a disc/pads system in sliding contact[END_REF].

It is important to note that only the initial conditions are projected on the basis T 0 and therefore only the iterative process of the shooting method uses the reduction. The process for the determination of the nonlinear solution Z (Z 0 , t) is kept unchanged: it is calculated by the modified θ-method on a discrete time basis of the interval [0; T ]. This implies that there is no approximation of the solution Z (Z 0 , t) through this reduction.

From a practical point of view, T 0 is orthonormalized with respect to the mass matrix M and orthogonalized with respect to the stiffness matrix K, such that:

T 0 0 0 T 0 Q 0 = T 0 0 T Q (13)
where T T MT corresponds to the identity matrix and T T KT is a diagonal matrix. Q corresponds to the vector of sought modal amplitudes in displacements and velocities. It is writing in the form

Q = q 1 . . . q n p 1 . . . p n T (14)
The associated basis T is defined by

T = Φ 1 . . . Φ m .
As previously explained, the use of this reduced basis of initial conditions implies only modifications of the iterative process of the shooting method. At each iteration, the corrections made to the vector Q and to the period T are calculated by projecting the system (7) on the basis T such that:

H q i + ∂H q ∂Q Q i ,T i ∆Q i + ∂H q ∂t Q i ,T i ∆T i = 0 ( 15
)
where H q corresponds to the projection of the condition of periodicity on the basis T. This vector of size 2m is given by:

H q j = Φ T j M(U(Q, t = T ) -U stat ) -q j H q j+m = Φ T j M U(Q, t = T ) -p j ∀j : 1 ≤ j ≤ m (16)
This new procedure implies also the calculation of the two following partial derivatives H q (Q, T ) and ∂H q /∂Q.

The first partial derivative H q (Q, T ) of size 2m is now given by:

∂H q ∂t Q i ,T i = ∂Z q (Q i , t) ∂t t=T i = g q Z(Q i , T i ) ( 17
)
where g q is the projection of the equation of motion into state variable on the basis T and Z q is the projection of the state vector on the basis T. The complete expression of g q is given by:

g q j = Φ T j M U(Q, T ) g q j+m = Φ T j F + R(Q, T )-C U(Q, T ) -K(U(Q, T ) -U stat ∀j : 1 ≤ j ≤ m (18)
The second partial derivative ∂H q /∂Q of size 2m × 2m is given by

∂H q ∂Q Q i ,T i = ∂Z q ∂Q Q i ,T i -I 2m (19)
where I 2m is the 2m × 2m identity matrix and ∂Z q ∂Q corresponds to the reduced monodromy matrix (i.e.

projected on the basis of T): it represents the variation of the reduced amplitudes with reduced initial conditions Q. It is a T -periodic matrix of size 2m × 2m. For more clarity, it is noted J q in the rest of this paper. Each column of this matrix J q can be calculated using a numerical time integration of the solution from perturbed reduced initial conditions of the form:

Q + δQ k e k ∀k : 1 ≤ k ≤ 2m (20)
where e k corresponds to the k th column of the identity matrix, Q k defines the k th component of Q and δ is the perturbation. Then, each column of the matrix J q can be calculated by using the following relations:

J q | k j = Φ T j M(U(Q + δQ k e k , T ) -U(Q, T )) δQ k J q | k j+m = Φ T j M( U(Q + δQ k e k , T ) -U(Q, T )) δQ k        ∀j : 1 ≤ j ≤ m (21)
This iterative process is repeated by considering the corrected initial condition and period given by:

Q i+1 = Q i + ∆Q i (22)
T i+1 = T i + ∆T i (23)
until the two convergence criteria are satisfied:

m j=1 Φ T j KΦ j H q j 2 + H q j+m 2 m j=1 Φ T j KΦ j q j 2 + p j 2 < 1 ( 24
)
|∆T i |/|T i+1 | < 2 (25)
where 1 and 2 are chosen residuals. The first criterion is an energetic adaptation of the criterion previously proposed in Eq. [START_REF] Kinkaid | Automotive disc brake squeal[END_REF]. It ensures a similar convergence for both displacements and velocities.

It can be noted that expressions [START_REF] Hoffmann | A minimal models for studying properties of the mode-coupling type instability in friction induced oscillations[END_REF] and ( 19) correspond to approximations of the two partial derivatives given in Eqs. ( 8) and [START_REF] Ibrahim | Friction-induced vibration, chatter, squeal and chaos part ii : Dynamics and modelling[END_REF]. It is important to remember that no other approximation is made in the iterative process of the shooting method by applying the proposed methodology: the determination of the nonlinear solution Z (Z 0 , t) is unaltered.

Initialization of the shooting method

A relevant initial condition is essential for the computational efficiency of the shooting method. Charroyer et al. [START_REF] Charroyer | Self-excited vibrations of a non-smooth contact dynamical system with planar friction based on the shooting method[END_REF] proposed a new approach for the initialization process. In the following, this process is briefly described. More details can be found in [START_REF] Charroyer | Self-excited vibrations of a non-smooth contact dynamical system with planar friction based on the shooting method[END_REF].

It is first admitted that the solution U(t) can be approximated by considering only the contribution of the unstable mode such that:

U(t) = q Φe i (λ)t (26)
where (.) and (.) denote the real and imaginary parts of the variables. λ corresponds to the eigenvalue of the unstable mode that has been previously calculated via the stability analysis. Φ corresponds to the associated eigenvector and q is the modal amplitude. Therefore the velocity U(t) is given by:

U(t) = q i (λ)Φe i (λ)t (27)
Then the initial conditions are defined by

Z 0 0 = q 0 (Φ) + Φ E q 0 (λΦ) (28)
T 0 = 2π/ (λ) [START_REF] Slater | A numerical method for determining nonlinear normal modes[END_REF] where Φ E corresponds to the steady sliding equilibrium of the nonlinear system that has been previously calculated for the stability analysis. The value q 0 is an unknown real value that has to be determined in order to obtain the most relevant and optimal initial conditions for the proposed formulation. The calculation of this value q 0 is performed by considering that the vibrational amplitudes growth until the occurrence of nonlinearities.

The approximated solutions ( 26) and ( 27) are valid until the occurrence of nonlinearities and the stabilization of the solution due to the saturation of the frictional contact reactions and a decrease of the average injected power ratio due to the contact reactions. So the calculation of the optimal value q 0 is performed by considering the power balance of the steady state. It was previously demonstrated in [START_REF] Charroyer | Self-excited vibrations of a non-smooth contact dynamical system with planar friction based on the shooting method[END_REF] that this power balance can be expressed by the following expression: τC = τR [START_REF] Kuether | A numerical approach to directly compute nonlinear normal modes of geometrically nonlinear finite element models[END_REF] where τC corresponds to the average power ratio dissipated by damping and τR defines the average power ratio injected by the frictional contact reactions. Charroyer et al. [START_REF] Charroyer | Self-excited vibrations of a non-smooth contact dynamical system with planar friction based on the shooting method[END_REF] demonstrated that τC is independent of the modal amplitude q whereas τR depends on the modal amplitude q and can be easily calculated by saturating the reactions corresponding to the linearized sliding contact problem.

In conclusion, a relevant initial condition for the shooting method can be obtained by using the approximated linear solution [START_REF] Seydel | Practical Bifurcation and Stability Analysis[END_REF] and by estimating the optimal value q 0 that verifies Eq. ( 30). The initial condition for the unknown period of the solution is given by Eq. ( 29).

Results

In this section, an application of the methodology presented in the previous section is proposed. A stability analysis is first briefly investigated. This allows us to choose a set of parameters for which the brake system exhibits one instability. Then, the shooting method with two specific reductions for the initial guess is carried out for this specific configuration in order to determine the nonlinear self-sustained vibration. Results are discussed and compared with the reference solution obtained by using the full time integration method.

Stability analysis and choice of the configuration for the investigation of nonlinear vibration

A stability analysis is performed for various friction coefficients µ = [0; 0.5]. Only one unstable mode is detected in the frequency range of interest. Fig. 3 shows the mode coupling pattern and the evolution of the associated real parts. Considering these results, the configuration with µ = 0.2 is chosen in the next section in order to predict the nonlinear self-sustained vibration. The associated stability chart is given in Fig. 4. As previously shown in Fig. 3, one unstable mode is depicted at 3382.5 Hz with a value of real part equal to 17 s -1 . Moreover Fig. 5 show the sliding equilibrium of the brake system and the unstable mode shape for µ = 0.2. In this figure, the dotted black lines correspond to the initial position of the brake system. Table 2. Error on the mode shapes based on the 1 -MAC criteria at four specific times over one period

In this section, the stationary regimes obtained by using the reduced shooting method are compared with those provided by the full time integration. The global vibrational behaviour of the brake system during steady state regimes is firstly discussed. Then, the mechanical energy over a period and the associated spectrum is presented. Finally, the limit cycle for a specific node selected on the frictional interface is analysed.

For this application, a simple basis is used in the reduced shooting method. It is only composed of the quasi-static equilibrium Φ E , the real part Φ R and the imaginary part Φ I of the eigenvector of the unstable mode resulting from the stability analysis. This reduced basis is denoted by

T 0 = Φ E Φ R Φ I .
Considering relation [START_REF] Tison | Improvement in the predictivity of squeal simulations: Uncertainty and robustness[END_REF], Q 0 is given by Q 0 = q E q R q I p E p R p I where q E , q R and q I (respectively p E , p R and p I ) correspond to the displacement amplitudes (respectively velocity amplitudes) of the modes Φ E , Φ R and Φ I . This basic choice makes it possible to take into account the potential evolution of the equilibrium position previously established by the vibration stability analysis. Thereby, the static field vector U stat defined in the Eq. ( 12) is null for this application because the proposed reduced basis already contains the equilibrium of the system. In addition, it contains only an initialization on the unstable mode that generates the self-sustaining vibrations. So this choice corresponds to one of the most reduced basis.

Moreover, a relevant initial condition is chosen for the shooting method. Following the process explained in section 4, Q 0 is initialized as follows:

Q 0 = 1 q 0 0 0 (λ) q 0 -(λ) q 0 (31)
where (λ) and (λ) are the real and imaginary parts of eigenvalue for the unstable mode and q 0 corresponds to the optimal value of q that verifies Eq. [START_REF] Kuether | A numerical approach to directly compute nonlinear normal modes of geometrically nonlinear finite element models[END_REF]. First of all, Fig. 6 shows the evolution of the self-sustained steady state vibration of the brake system at different times uniformly distributed over one period. A strong similarity is observed between the results obtained by the proposed shooting process and the classical integration scheme. The fundamental frequency is estimated at 3371.9 Hz by the shooting method, which is very close to that calculated by time integration, which is estimated at 3376.2 Hz (i.e. error of less than 0.13%). Concerning more specifically the convergence of the proposed shooting method, value of the two criteria 1 and 2 has been chosen equal to = 10 -3 . This leads to only 3 iterations in order to obtain the convergence, which demonstrates the relevance of the choice of the initial guess for the prediction of the steady state nonlinear solution. The choice of 10 -3 for the residuals is motivated by the results of a preliminary study on a simplified model [START_REF] Charroyer | Self-excited vibrations of a non-smooth contact dynamical system with planar friction based on the shooting method[END_REF]. In addition, the purpose of the method is to quickly obtain an acceptable approximation of the non-linear regime. The comparisons with the exact solution presented below show that this objective is achieved with this choice of residuals.

In order to better compare the shapes of the nonlinear solution estimated by the shooting method with shapes of the reference solution computed, we use the classical Modal Assurance Criterion (MAC) relative to the mass matrix M given by:

MAC(A, B) = (A T MB) 2 (A T MA)(B T MB) (32)
where A and B are the matrix of the mode shapes (at a specific time) of a reference solution and of the solution computed via the reduced shooting method. Tab. 2 gives the 1 -MAC values for the four different times previously selected for the visualization of nonlinear vibrations of the complete brake system. A low value of the 1 -MAC indicates that the vibration mode shapes between the reference and the approximated solution are well correlated. A high value indicates that mode shapes are orthogonal and uncorrelated so that the shooting method does not work properly to estimate the nonlinear self-sustained vibratory behaviour.

The quantitative error indicators 1 -MAC given in Tab. 2 are less than 7%. This shows that agreement between the reference solution and the proposed shooting process is very good, which demonstrates the ability of this method to estimate the nonlinear self-excited vibration. Fig. 7(a) gives the total mechanical energy and its different harmonic components computed from a Fourier transform of the solution over one period. Fig. 7(b) represents the corresponding relative error between the reference and the nonlinear response via the shooting method. Most of the mechanical energy is distributed between the static and the fundamental components (see Fig. 7(a) for the 0× and 1× harmonic components. These two main contributions are very well estimated by the shooting method with an relative error of less than 0.6% (2.4%, respectively) for the static component (the fundamental frequency, respectively). As shown in 7(a), energy distribution trends on each harmonic components provided by the shooting method (in red color) remain representative of those of the reference (in black color).

Tab. 3 gives the relative variation of mechanical energy over one period ∂E ∂t / Ē as well as the different average power rates: the power rate dissipated by damping τc , the power rate injected by unstable friction τf and the power rate injected during impacts τn . For this last power rate, a negative sign indicates that the energy is in fact dissipated, which is a expected result due to the velocity discontinuities during these inelastic impacts [START_REF] Charroyer | Self-excited vibrations of a non-smooth contact dynamical system with planar friction based on the shooting method[END_REF][START_REF] Acary | Energy conservation and dissipation properties of time-integration methods for nonsmooth elastodynamics with contact[END_REF]. First the relative average power balance is verified (i.e. ∂E ∂t / Ē + τc = τf + τn ). This emphasizes the excellent behaviour of the numerical time integration, which results from the energy conservation properties of the modified θ-method scheme [START_REF] Acary | Energy conservation and dissipation properties of time-integration methods for nonsmooth elastodynamics with contact[END_REF]. Then, it can be observed that the variation of the mechanical energy is almost zero over one period T so that the average power ratio τc dissipated by damping is equal to the average power ratios τf and τn injected by the frictional and contact reactions. This confirms the hypothesis that the nonlinear contact phenomena of separation and impact between the disc and the pads are responsible for the convergence and stabilization towards the nonlinear self-excited solution of the frictional mechanical system. Indeed, the nonlinearities induce a saturation of the frictional contact reactions and a decrease of the average power injected in the system compared with the unstable linear case. In order to better visualize theses nonlinear phenomena such as the separation between the pads and disc at the frictional interface, Fig. 8 shows the evolution of the the upper contact zone at eight specific times over one period. It appears very clearly that separation and contact occur intermittently on several pads. Fig. 9 illustrates the phase diagrams according to the tangential, radial and normal directions for a node of the contact interface (node 1864 as indicated in Fig. 8(h)). It is clearly observed that the results based on the proposed shooting method with a reduced basis for the initial guess are in good agreement with those obtained from the direct time integration. Moreover it is shown that a complex nonlinear behavior such as impact discontinuity can be reproduced (see Fig. 6(c)). The difference in amplitude of the phase diagrams seems to be large with respect to the obtained overall percentage error. However, phase diagrams represent very localized contact responses while energy errors are averaged over the entire structure. It appears that contact responses can present large errors without major effects on the rest of the structure provided that most of the phenomena are well represented. In this case, the phenomena of contact separation allow the system to stabilize on a periodic self-sustaining vibration, even if they are of smaller amplitudes (as shown by the n-limit cycle in Fig. 6(c)). The amplitude and overall shape of the structural vibration is then reached with a good degree of precision.

Finally, the combination of the small number of iterations (i.e. 3) and the small number of modes retained in the reduction basis (3 vs 31974 degrees of freedom of the unreduced system) leads to a large computational gain between the proposed method and the direct time integration (20 min vs 10 h on a standard workstation).

Conclusion

In this paper, a new adaptation of the shooting method is proposed in order to estimate nonlinear selfsustained vibration of large finite element models for which the number of degrees of freedom is important and prevents the application of the conventional shooting method. The originality of the proposed formulation is based on the use of a reduced basis for the determination of the initial conditions. This dramatically reduces the number of variables to be estimated during the iterative process of the shooting method and as a result, the size of the monodromy matrix to be computed, while keeping the resolution of the full original system for the determination of the nonlinear solution. This proposed numerical approach has been tested in the case of a finite element model for the prediction of nonlinear dynamic steady-state solutions of autonomous non-smooth contact systems prone to mono-instability.

All the numerical results indicate that the proposed nonlinear method appears to be efficient for a fairly accurate estimation of the self-sustained vibrations as well as the nonlinear phenomena at the frictional interfaces.

An interesting future study to be considered is the extension of the proposed numerical method in the case of multi-instabilities since experimental results show that squeal noise can have several fundamental frequencies. This case is called quasi-periodic where the evolution of the system is subject to an instability with several unstable modes with incommensurable frequencies. Additional developments to the shooting method proposed in this paper would be necessary in order to approximate the stationary regime of such systems. The definition of a criterion for strict quasi-periodicity in the time domain is not obvious and it should be necessary to use a more flexible criterion (for instance stationarity) or to decide in the frequency domain. In addition, the initialization of the shooting method on a set of unstable mode contributions requires the search for modal amplitudes but also phase shifts between modes. However, preliminary work can consist in seeking periodic solutions, even for a multi-unstable system. Indeed, Loyer et al. [START_REF] Loyer | Study of nonlinear behaviors and modal reductions for friction destabilized systems. application to an elastic layer[END_REF] have shown that in some cases, only one dominant mode can remain in the stationary solution since the contributions of the other modes have been stabilized or even annihilated in the transitional regime. All the proposed methodology is then relevant for this specific objective.

Fig. 1 .Fig. 2 .

 12 Fig. 1. Finite element model of the simple brake system (a) the hub-disc assembly (b) the 12 small cylindrical pins applied on either sides of the disc and the supporting structures modeled as elastically-restrained rigid bodies

Fig. 3 .FrequencyFig. 4 . 2 -Fig. 5 .

 3425 Fig. 3. Evolution of the coupling modes versus the friction coefficient (a) real part (b) frequency

Fig. 6 .Fig. 7 .

 67 Fig.6. Evolution of the self-sustained steady state vibration of the brake system for the reference (a,c,e,g) and the shooting method (b,d,f,h)

Fig. 8 .

 8 Fig. 8. Evolution of the contact and separation states over one period at the upper frictional interface between the disc and the 6 cylindrical pins during the steady state vibrations for the reference case

Fig. 9 .

 9 Fig.9. Phase diagrams of node 1864 provided by the direct time integration (red with a star) and limit cycles given by the shooting method in the (a) tangential (b) radial and (c) normal directions

Table 1 .

 1 Material and geometrical properties of the brake system of the nonlinear vibration.

	F ext	2000 N

8

N/rad Normal force applied on the caliper support

Table 3 .

 3 Variation of the mechanical energy and average power rates on one period of the steady-state vibrations

	Quantity	Value
	Variation of mechanical energy ∂E ∂t / Ē	0.02765 s -1
	Average power rate dissipated by damping τc	363.7 s -1
	Average power rate injected by friction τf	371.62 s -1
	Average power rate injected by impact τn	-7.8962 s -1
	Average power balance	≈ 0

Acknowledgment

This work was performed within the framework of the LABEX CeLyA (ANR-10-LABX-0060) of Université de Lyon, within the program "Investissements d'Avenir" (ANR-11-IDEX-0007) operated by the French National Research Agency (ANR).

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