N Coudeyras

J-J Sinou

S Nacivet

A new treatment for predicting the self-excited vibrations of nonlinear systems with frictional interfaces: the Constrained Harmonic Balance Method (CHBM). Application to disc brake squeal

Noise brake squeal is still an issue since it generates high warranty costs for the automotive industry and irritation for customers. Key parameters must be known in order to reduce it. Stability analysis is a common way of studying nonlinear phenomena and has been widely used by the scientific and the engineering communities for solving disc brake squeal problems. This type of analysis provides areas of stability versus instability for driven parameters, thereby making it possible to define design criteria. Nevertheless, this technique does not permit obtaining the vibrating state of the brake system and nonlinear methods have to be employed. Temporal integration is a well known method for computing the dynamic solution but as it is time consuming, nonlinear methods such the Harmonic Balance Method are preferred. This paper presents a novel nonlinear method called the Constrained Harmonic Balance Method (CHBM) that works for nonlinear systems subject to flutter instability. An additional constraining based condition is proposed that omits the static equilibrium point (i.e. the trivial static solution of the nonlinear problem that would be obtained by applying the classical Harmonic Balance Method), and therefore focus on predicting both the Fourier coefficients and the fundamental frequency of the stationary nonlinear system. The effectiveness of the proposed nonlinear approach is illustrated by an analysis of disc brake squeal. The brake system under consideration is a reduced finite element model of a pad and a disc. Both stability and nonlinear analyses are performed and the results are compared with a classical variable order solver integration algorithm. Therefore the objectives of the following paper are to present not only an extension of the Harmonic Balance Method (the Constrained Harmonic Balance Method -CHBM), but also to demonstrate an application to the specific problem of disc brake squeal with extensively parametric studies that investigate the effects of the friction coefficient, piston pressure,nonlinear stiffness and structural damping.

Introduction

Disc brake squeal is a complex phenomenon and has been a challenging issue for many researchers for a long time. There is no precise definition of brake squeal [START_REF] Kinkaid | Automotive disc brake squeal[END_REF], but it could be defined as a monoharmonic 1 sound emitted at a frequency over a range of 1 kHz to 20 kHz during a braking event. Squeal is a fugitive noise, i.e. at each braking action a brake system may or may not be noisy. Although brake squeal is undesirable, it is the sign that the brake is reliable [START_REF] North | Disc brake squeal, a theoretical model[END_REF]. Ouyang et al. [START_REF] Ouyang | Friction-induced parametric resonances in discs: effect of a negative friction-velocity relationship[END_REF] demonstrated that parametric resonances can occur when an elastic element is rotated around an annular disc with friction having a negative slope with velocity. They also indicated an elastic system can oscillate in the stickslip mode (that is a low sliding speed phenomenon caused when the static friction coefficient is higher than the dynamic coefficient) in the plane of a disc due to nonsmooth friction nonlinearity [START_REF] Ouyang | Friction-induced vibration of an elastic slider on a vibrating disc[END_REF]. Spurr [START_REF] Spurr | A theory of brake squeal[END_REF] analyzed squeal as a sprag-slip phenomenon. He indicated that the tribological property cannot be considered as the only reason for brake squeal, and that vibration could occur when the friction coefficient remains fairly constant with speed. Later, the sprag slip phenomenon was generalized as a geometrically induced or kinematic constraint instability. For example, Jarvis and Mills [START_REF] Jarvis | Vibration induced by dry friction[END_REF] and Millner [START_REF] Millner | Analysis of disc brake squeal[END_REF] worked on mass-spring models and showed that autonomous vibrations are due to friction that couples two modes together, thus instability occurs even if the friction coefficient is constant. Brake squeal has been identified by Oden and Martins [START_REF] Oden | Model and computational methods for dynamic friction phenomena[END_REF] as a result of friction-induced vibration. Hence, if friction force couples two degrees of freedom (dof), unstable modes could merge and generate squeal. Liles [START_REF] Liles | Analysis of disc brake squeal using finite element methods[END_REF] studied a large finite element model and confirmed that brake squeal is due to the friction coupling effect, leading to mode coalescence. The friction coefficient appears to be the essential parameter for detecting squeal phenomena. Moreover, Ouyang et al. [START_REF] Ouyang | Methodology for the determination of dynamic instabilities in a car disc brake[END_REF] proposed to study the stability analysis of a car disc brake system (with pads, calliper and mounting) by considering a combined analytical and numerical method that uses the finite element method.

Stability analysis is a classical method for studying the brake squeal phenomenon [START_REF] Liles | Analysis of disc brake squeal using finite element methods[END_REF][START_REF] Moirot | Etude de la stabilité d'un équilibre en présence de frottement de coulomb, application au crissement des freins à disque[END_REF][START_REF] Nack | Brake squeal analysis by finite elements[END_REF][START_REF] Kung | Complex eigenvalues analysis for reducing low frequency brake squeal[END_REF][START_REF] Lorang | Tgv disc brake squeal[END_REF]]. An analytical finite element model with nonlinear contents such as contact and frictional elements is considered. A complex eigenvalue computation of the respective linearized system is then performed, followed by a study of the corresponding real parts. A positive real part indicates that the corresponding eigenmode is unstable and squeal may occur. Parametric studies are carried out and several design criteria are derived. However, as mentioned by Ouyang et al. [START_REF] Ouyang | Numerical analysis of automotive disc brake squeal : a review[END_REF], eigenvalue calculation is insufficient due to linearization which provides valid results only close to the steady sliding state. The real part of an eigenvalue indicates the growth rate of oscillations; however it does not provide information on the amplitude of the dynamic response [START_REF] Liles | Analysis of disc brake squeal using finite element methods[END_REF]. Moreover, eigenvalues analysis overestimates the number of unstable modes and they cannot all be observed in experiments [START_REF] Massi | Linear and non-linear numerical approaches[END_REF]. Finally, the starting vibration mechanism is unknown. As a result, transient analysis is the natural second step in studying brake squeal. Contrary to eigenvalue analysis, transient analysis can include nonlinear aspects of the model. The models can be refined and the use of time-dependent loading, sophisticated friction laws and so on is possible. Better qualitative and quantitative results are derived, considerably contributing to the improvement of brake systems. A large number of transient analyses of finite element models have been performed in the past. Nagy et al. [START_REF] Nagy | A new method development to predict squeal occurrence[END_REF] was one of the first to perform a numerical integration in a finite element disc brake. Chargin et al [START_REF] Dune | Nonlinear dynamics of brake squeal[END_REF] carried out a transient computation on a very simple brake system by using an implicit integration scheme. Mahajan et al [START_REF] Mahajan | Vehicle disc brake squeal simulation and experiences[END_REF] ran both complex and temporal approaches and found that both methods are useful for design modifications. Hu et al [START_REF] Hu | Squeal doe using non-linear transient analysis[END_REF] performed an explicit time integration analysis combined with Taguchi's method and found that the characteristics of friction materials are an important factor, along with rotor thickness, pad chamfer and pad slot. More recently, Massi et al [START_REF] Massi | Linear and non-linear numerical approaches[END_REF] performed a dynamic transient computation on a large dofs disc brake model with an in-house finite element code and correlated the modal complex analysis with the time simulation in the sense that the vibrating steady state matched one of the unstable modes found in the complex analysis. The major drawback of the works mentioned above concerns the excessive calculation time required to obtain the stationary state of oscillations, which penalizes design modifications. Shorter computation times can be achieved with these methods at the cost of over-simplified finite element models. AbuBakar and Ouyang [START_REF] Abubakar | Complex eigenvalue analysis and dynamic transient analysis in predicting disc brake squeal[END_REF] performed a transient analysis for three different contact regimes in ABAQUS and found only one that gives acceptable results, i.e. the oscillation frequency is similar to one of the results obtained by complex analysis. Computation took more than 24 hours before calculations diverged. One way of tackling these drawbacks is to carry out the solution in the frequency domain. A very well-known method for solving nonlinear problems is the Harmonic Balance Method (HBM) in combination with the Alternate Frequency Time Domain Method proposed by Cameron and Griffin [START_REF] Cameron | An alternating frequency/time domain method for calculating the steady-state response of nonlinear dynamic systems[END_REF]. This method has been used by many authors to solve nonlinear problems [START_REF] Pierre | Multi-harmonic analysis of dry friction damped systems using an incremental harmonic balance method[END_REF][START_REF] Cheung | Application of the incremental harmonic balance method to cubic nonlinearity systems[END_REF][START_REF] Groll | The harmonic balance method with arclength continuation in rotor/stator contact problems[END_REF][START_REF] Sinou | A non-linear study of a cracked rotor[END_REF]. The key factor of HBM is the computation of the steady-state solution without the transient part. HBM is well designed for systems under periodic excitations. It is less time consuming and requires less disc storage. In the particular case of self-excited systems subjected to Hopf bifurcations, it is a little more complicated to apply this method since the uniqueness of the solution is lost [START_REF] Moirot | Etude de la stabilité d'un équilibre en présence de frottement de coulomb, application au crissement des freins à disque[END_REF] and both static and dynamic solutions coexist. Hence the system is driven only by initial conditions and leads to a unique final solution. In an optimization domain, that corresponds to two local minima the solver computes either the static solution or the dynamical one without any control, but the static solution is always reached whatever the initial conditions. This reflects the fact that the static solution corresponds to an "exact solution" of the nonlinear system in the Fourier domain (the solution is only composed by the static Fourier coefficients), whereas the dynamic solution will be an approximation of the nonlinear system due to the truncated Fourier series. This is a major drawback that has been tackled in this study. In this paper, we propose a novel nonlinear approach, called the Constrained Harmonic Balance Method (CHBM) that works for nonlinear systems subjected to flutter instability. An additional constraining based condition is proposed for predicting both the Fourier coefficients and the fundamental frequency of the stationary nonlinear dynamic system amplitudes, called the "limit cycles amplitudes". This paper is divided into four sections. The first one deals with the presentation of the brake system under study. Secondly, the stability analysis of the nonlinear brake system is performed. The third part concerns the nonlinear analysis with the CHBM and the results are compared with a classical temporal integration scheme. The advantages and drawbacks of both methods are discussed. The last one is devoted to parameter analyses where the advantages of the new CHBM are illustrated.

Finite Element Model of the brake system

Figure 1 shows the finite element model of the car front brake under consideration, developed using the ABAQUS finite element software package. The model consists of the two main components contributing to squeal: the disc (see Figure 1(a)) and the pad (see Figure 1(b)). There are about 60000 nodes and tennode quadratic tetrahedron elements are used. They are very useful for meshing complex shapes and are second order elements which provide accurate results without requiring very fine meshing [27].

Model reduction

As seen previously in Figure 1, finite element models of the two brake components need a large number of dofs to represent geometrical details of the brake system. One of the first classical processes is to reduce the finite element models of the pad and disc by using a Craig and Bampton technique [START_REF] Craig | Coupling of substructures for dynamic analyses[END_REF] for keeping certain contact nodes and generalized dofs.

It consists in building a projection basis combining constraint modes T C and a truncated basis T ϕ,N of normal modes computed with a fixed interface. Hence the relationship between physical boundaries (or interface dofs) u B , interior dofs u I , modal coordinates of constraint modes q B and modal coordinates of normal modes q ϕ is given by:

u = u B u I = T C T ϕ,N q B q ϕ (1)
The constraint modes T C are computed by assuming unit displacements onto u B .

K BB K BI K IB K II < u B > u I = R B < 0 > (2)
where <> denotes prescribed quantities. Then the corresponding under-space for the static condensation is written as:

T C = I -K II -1 K IB ( 3
)
where I defines the identity matrix. T ϕ is achieved by resolving an eigenvalue problem with fixed interface dofs u B .

-

ω 2 M BB M BI M IB M II + K BB K BI K IB K II < 0 > ϕ I = R ϕ < 0 > (4)
where <> denotes prescribed quantities. Finally T ϕ,N is deduced by retaining N modes computed in Equation ( 4).

T ϕ,N = 0 ϕ 1:N,I (5)
Finally, stiffness matrix K and mass matrix M are given by

K = T C T ϕ,N T K BB K BI K IB K II T C T ϕ,N = KCC 0 0 KNN (6)
M = T C T ϕ,N T M BB M BI M IB M II T C T ϕ,N = MCC MCN MNC MNN (7)
where 0 defines the zero matrix. For the sake of convenience, in the following of the paper the hat above the reduced matrices is deleted. Nine contact nodes are kept on each structure at the disc/pad interface and four extra nodes are retained on the back-pad where piston force is applied. Moreover, the first fifty modes of each structure are held. Boundary conditions are achieved by embedding disc in the hub while the pad is only free to translate in the normal contact direction. The resulting model is a 158 dofs system.

Non-linear system

Many contact definitions could be used to model the contact between structures in finite element models but the simplest is the the penalty method mentioned by Ouyang in his review [START_REF] Ouyang | Numerical analysis of automotive disc brake squeal : a review[END_REF]. It consists in adding contact stiffness at the disc/pad interface. The frictional material is a mixture of many components and is about one thousand times less stiff the disc, thus its deformation under loading is greater and nonlinearity behavior occurs. Thus contact stiffnesses values are chosen to fit the first and the third order of the pad compression curves obtained in the tests. The computation of the nonlinear contact force takes the form:

F contact = k l (U i -U j ) + k nl (U i -U j ) 3 if (U i -U j ) > 0 0 otherwise (8)
with k l and k nl linear and nonlinear contact stiffnesses respectively and U i and U j displacements of interfaces nodes i and j that are related to the pad and the disc respectively. The friction forces are deduced from the contact forces by using the classical Coulomb law. A permanent sliding state is considered and a constant friction coefficient is assumed:

F friction = µF contact sgn(v r ) (9)
with µ being the friction coefficient and v r the relative velocity between both bodies. Thus the vectors of the nonlinear forces take the form:

F nl (x) = F contact (x) + F friction (x) (10)
Hence the reduced final model of the brake system is written as:

Mẍ + C ẋ + Kx + F nl (x) = F piston ( 11
)
where x is the response displacement of dofs, the dot denotes the derivation with respect to time, M, C, K are respectively the mass, damping and stiffness matrices of the system. F piston is the piston pressure force and vector F nl (x) corresponds to nonlinear forces. C is built by projecting the modal damping matrix D onto the undamped, non-frictional inverse modal basis Φ -1 of the reduced model:

C = Φ -1 T DΦ -1 (12)
The modal damping matrix D is built so that modal damping is added on both the modes involved in the instability. Initially, equal damping distribution is considered.

D = diag (0 . . . 0 D1 D2 0 . . . 0) (13)
with D1 = D2 = 1;

Stability Analysis

As mentioned before, the stability analysis is the first step for studying nonlinear systems subjected to instability phenomena. For a given set of parameters, a static equilibrium position can become unstable and stationary periodic oscillations, called limit cycles, occur. This analysis is performed in two steps [START_REF] Sinou | Friction-induced vibration for an aircraft brake system-part 1: Experimental approach and stability analysis[END_REF]. The first one consists in performing the static loading of the system defined in Equation [START_REF] Moirot | Etude de la stabilité d'un équilibre en présence de frottement de coulomb, application au crissement des freins à disque[END_REF] corresponding in an action by the driver on the brake pedal. Pressure is applied onto the piston which acts on the pads entering into contact with the disc. The corresponding nonlinear static equation is written as:

Kx 0 + F nl (x 0 ) = F piston (14)
where x 0 corresponds to the static equilibrium of the nonlinear brake system. Note that the static equilibrium is achieved with a non-zero rotational disc speed involving friction forces.

Then the system ( 11) is linearized about the static equilibrium position x 0 by using the perturbation technique. Such a perturbation:

x = x 0 + x (15)
Replacing Equation [START_REF] Ouyang | Numerical analysis of automotive disc brake squeal : a review[END_REF] into Equation [START_REF] Moirot | Etude de la stabilité d'un équilibre en présence de frottement de coulomb, application au crissement des freins à disque[END_REF] leads to:

M ẍ + C ẋ + K (x 0 + x) + F nl (x 0 + x) = F piston (16)
Supposing that F nl belong to the C 1 class, developing the nonlinear force as a Taylor series and keeping the first order leads to:

F nl (x 0 + x) ≈ F nl (x 0 ) + K L NL x (17)
K L NL corresponds to the linearized nonlinear forces Jacobian matrix which is composed by the following elementary matrix for each contact element:

K L NL ij =     ∂F nl i ∂x i x 0 ∂F nl i ∂x j x 0 ∂F nl j ∂x i x 0 ∂F nl j ∂x j x 0     (18)
And writing k l i = ∂F nl i ∂x i x 0 , the elementary stiffness matrix takes the following form:

K L NL ij =         0 -µk l i 0 0 µk l i 0 0 k l i 0 0 -k l i 0 0 0 0 0 0 0 0 µk l i 0 0 -µk l i 0 0 -k l i 0 0 k l i 0 0 0 0 0 0 0         (19)
It should be pointed out that the friction force is only written in the x direction that corresponds to the longitudinal direction of the vehicle. Hence the nonlinear system ( 11) is approximated at the equilibrium position by the following linearized one: The previous model ( 20) is then written in the following state-space and complex eigenvalues are derived :

M ẍ + C ẋ + K + K L NL x = 0 (20)
A = 0 I -M -1 K + K L NL -M -1 C (21)
Since the stiffness matrix ( 19) is asymmetrical due to the contribution of friction forces, the computed eigenvalues are complex and are written as:

λ = a + iω ( 22
)
where a is the real part of the eigenvalue that corresponds to the growth rate of the amplitude and ω is the imaginary part of the eigenvalue that corresponds to the pulsation of the mode. A negative real part indicates that the corresponding mode is stable. In other words, a perturbation about the static equilibrium sliding state will not modify the equilibrium position of the system. A positive real part equivalent to a negative damping leads to an unstable mode. Thus modifying one of the parameters will induce growing oscillations about the static equilibrium position of the system until the dynamical steady-state is reached. Figure 2 shows evolutions of normalized real parts and normalized frequencies of the associated eigenvalue versus µ, which is normalized with respect to the Hopf bifurcation point µ 0 . As seen in Figure 2(a), the real parts curves split into two branches near the Hopf bifurcation point µ 0 . One goes towards the positive real part half-space and become positive whereas the other branch decreases and remains in the negative real part half-space. Figure 2(b) shows the typical lock-in phenomenon between both modes of the system with µ increasing. It can be seen that coalescence between the two modes is perfect since they are equally damped. The effects of equally and non-equally damped modes will be illustrated in the last part of this paper. It is then possible to define stable areas versus unstable areas of the linearized system for a given set of parameters. In the following, the paper is devoted to nonlinear dynamic computation.

Non-linear dynamic and self-excited limit cycles

The classical approach of nonlinear analysis consists in using temporal integration schemes to compute nonlinear dynamic solutions. However, it can be observed that this kind of method is costly in terms of computation time and resources for large finite element models. Enhanced nonlinear methods have to be employed to save time. Extensive reviews on this topic have been given in [START_REF] Sinou | Methods to reduce non-linear mechanical systems for instability computation[END_REF]. In the following part of the paper, an original adaptation of the harmonic balance method for self-excited nonlinear systems will be introduced and discussed.

The Constrained Harmonic Balance Method

In this section, we propose to introduce an extension of the harmonic balance method, called the Constrained Harmonic Balance Method (CHBM), for approximating stationary nonlinear responses of selfexcited systems subjected to flutter instabilities. Traditional harmonic balance methods are well-known numerical methods that have been commonly used to solve nonlinear problems in engineering [START_REF] Sinou | Methods to reduce non-linear mechanical systems for instability computation[END_REF]. However, they do not permit obtaining the stationary nonlinear vibrational responses of self-excited systems due to the fact that the static nonlinear solution corresponds to the trivial solution of the problem. In this section, the classical harmonic balance method with a condensation procedure on the nonlinear dofs is presented first. Then the additional constraining condition allowing the determination of the limit cycle amplitudes and an optimized initial condition process are presented and discussed.

The Harmonic Balance Method with a condensation procedure

Considering harmonic balance methods, a nonlinear solution is assumed to be a truncated Fourier series and the exact nonlinear periodic solution X (t) is replaced as:

X app (t) = N h k=0 U C k cos (kωt) + N h k=1 U S k sin (kωt) (23)
where U C k and U S k are vectors of Fourier coefficients and ω defines the final pulsation of the nonlinear limit cycles. It can be seen that ω is an unknown parameter in this study since we are in the presence of a self-excited system and the frequency of the stability analysis differs slightly from that of the nonlinear steady-state solution. Thus it cannot be used as a fixed parameter. N h is the number of the harmonic coefficients retained for the approximated nonlinear stationary solution. Velocities and accelerations are obtained by derivation of Equation [START_REF] Pierre | Multi-harmonic analysis of dry friction damped systems using an incremental harmonic balance method[END_REF] with respect to the time. The advantage of harmonic balance method is that it allows keeping only the first terms of Equation ( 23) where, generally, a preponderant energy part of the signal is concentrated. Replacing the approximated solution X app (t) into Equation ( 11) leads to

R N h (t) = N h k=0 K -(kω) 2 M U C k + (kωC) U S k cos (kωt) + N h k=1 K -(kω) 2 M U S k -(kωC) U C k sin (kωt) + F nl U C,S k -F piston (24)
Projecting the residue on sine and cosine orthonormal bases, and writing the multi-harmonics vector Ũ such that:

Ũ = U C 0 T U C 1 T U S 1 T • • • U C N h T U S N h T T (25)
leads to the following approximated equation:

Λ Ũ + Fnl Ũ = Fout (26)
with

Λ =           K 0 0 0 0 0 0 Λ h,1 0 0 0 0 0 0 . . . 0 0 0 0 0 0 Λ h,k 0 0 0 0 0 0 . . . 0 0 0 0 0 0 Λ h,N h           (27)
and

Λ h,k = -(kω) 2 M + K kωC -kωC -(kω) 2 M + K for k = 1 : N h (28)
Λ h,k is the dynamical stiffness matrix associated with the k th harmonic and Fout are the external forces. Equation ( 25) gathers the Fourier coefficients that have to be balanced to obtain the periodic solution of the nonlinear system. Non-linear force Fourier coefficients depend on Ũ and their determination can be fastidious analytically because of the size of the system and the number of harmonics. The Alternate Frequency Time Domain Method proposed by Cameron and Griffin [START_REF] Cameron | An alternating frequency/time domain method for calculating the steady-state response of nonlinear dynamic systems[END_REF] permits omitting this issue as outlined below:

Ũ FFT -1 -→ X (t) -→ F nl (X(t)) FFT -→ Fnl Ũ (29)
When a nonlinear system has a significant number of dofs but only a few of then are related to nonlinear components, it is possible to reduce system (26) on the nonlinear dofs without loss of accuracy [START_REF] Hanh | Harmonic balance analysis of general squeeze film damped multidegree-of-freedom rotor bearing systems[END_REF][START_REF] Villa | Stability and vibration analysis of a complex flexible rotor bearing system[END_REF]. Linear and nonlinear nodes are separated (i.e. the new vector is such that the nonlinear dofs are stored at the vector's end). Equation ( 26) may be rewritten in the following form

Λ ln,ln Λ ln,nl Λ nl,ln Λ nl,nl Ũln Ũnl + 0 Fnl = Fout,ln Fout,nl (30)
where U ln and U nl define the linear dofs and nonlinear dofs, respectively. F out,ln and F out,nl are the associated linear and nonlinear external forces. In the current model, only external linear force is available, i.e. the piston force. The purpose of the condensation aims at solving the algebraic nonlinear system of equations only for nonlinear dofs, leaving the other linear ones to be determined later by a linear transformation. Hence rewriting [START_REF] Sinou | Methods to reduce non-linear mechanical systems for instability computation[END_REF] on nonlinear dofs leads to:

Λ eq Ũnl + Fnl (U nl ) = Feq (31)
with

Λ eq = Λ nl,nl -Λ nl,ln (Λ ln,ln ) -1 Λ ln,nl (32)
and

Feq = Fout,nl -Λ nl,ln (Λ ln,ln ) -1 Fout,ln (33)
Thus system [START_REF] Hanh | Harmonic balance analysis of general squeeze film damped multidegree-of-freedom rotor bearing systems[END_REF] has the size of the number of nonlinear dofs and is lighter than [START_REF] Sinou | Methods to reduce non-linear mechanical systems for instability computation[END_REF]. Equation ( 31) is rewritten in the following form to be solved:

f Ũnl = Λ eq Ũnl + Fnl Ũnl -Feq (34)
When optimization is finished and Ũnl is known, linear displacements are obtained with:

Ũln = Λ -1 ln,ln Fout,ln -Λ ln,nl Ũnl (35)

The additional constraining condition

Equation ( 34) is a cost function that has a minimum when Ũnl is a solution of the system and can be solved by nonlinear least-square algorithms such as those of the Gauss-Newton and Leveberg-Marquardt methods. As stated before, the uniqueness of the solution is lost for systems at the Hopf bifurcation point [START_REF] Moirot | Etude de la stabilité d'un équilibre en présence de frottement de coulomb, application au crissement des freins à disque[END_REF]: the exact and trivial solution of Equation ( 34) corresponds to the static equilibrium point which is unstable. If the classical Harmonic balance Method is used, the only solution that will be found will be this static solution due to the fact that the residue of Equation ( 34) will be equal to zero for the static equilibrium point. So, in order to reject this trivial static solution and obtain the stationary nonlinear dynamical oscillations that correspond to the limit cycle amplitudes, it is necessary to add a constraint to the Harmonic Balance Method to reach only the minimum of Equation ( 34) which corresponds to the stationary nonlinear periodic motion. The constraining condition will be outlined in this paragraph. By considering the nonlinear autonomous system [START_REF] Moirot | Etude de la stabilité d'un équilibre en présence de frottement de coulomb, application au crissement des freins à disque[END_REF] and writing it in the state space gives:

Ẏ = AY + F P + F NL (Y) (36)
with

Y = X Ẋ , F P = 0 M -1 F piston , F NL = 0 -M -1 F nl (37)
and

A = 0 I -M -1 K -M -1 C (38)
A nonlinear periodic solution Y e (t) of Equation ( 36) is such that a real T exists so that:

Y e (t + T ) = Y e (t) and Y e (t + T ) = Y e (t) for 0 < T < T (39)
T is the period of the solution. It may be noted that T is an unknown parameter due to the absence of forced excitations and the difference between the stability frequency analysis and that of the nonlinear system.

By disturbing a solution Y e (t) with a perturbation ǫ(t) we obtain:

Y(t) = Y e (t) + ǫ(t) (40)
And by substituting Equation (40) in Equation [START_REF] Fritz | Investigation of the relationship between damping and mode-coupling patterns in case of brake squeal[END_REF], we obtain:

Ẏe + ǫ = A (Y e + ǫ) + F P + F NL (Y e + ǫ) (41)
By supposing that F NL is C 1 class, its development in the Taylor series at Y e at the first order gives

Ẏe + ǫ ≈ A (Y e + ǫ) + F P + F NL (Y e ) + J NL ǫ (42)
with J NL the Jacobian matrix of the first derivatives of the nonlinear forces F NL respect to the periodic solution Y e (t).

Since Y e (t) is the solution of Equation ( 36), Equation ( 42) is under the form

ǫ ≈ Aǫ + J NL ǫ = Jǫ (43) with J = A + J NL (44)
J is the Jacobian matrix of the nonlinear system [START_REF] Sinou | The influence of damping on the limit cycles for a self-exciting mechanism[END_REF] and depends on the dynamical solution Y e (t).

Thereby the eigenvalues of J define the evolution of the limit cycles amplitudes of the nonlinear system. If one or more eigenvalues are positive, the approximated nonlinear solution of the system is increasing and is governed by the unstable modes. If all the eigenvalues are negative, the nonlinear solution is decreasing (i.e. we are still in transient motion). If one eigenvalue is equal to zero whereas all the others are negative, the nonlinear approximated solution defines the stationary motion of the nonlinear system subjected to flutter instability. The replacement of the nonlinear contributions F NL by a linear approximation

J NL is done to mini- mize the difference ζ ζ = F NL (Y e (t)) -J NL Y e (t) (45)
This kind of transformation refers to the equivalent linearization concept proposed by Iwan [START_REF] Iwan | A generalization of the concept of equivalent linearization[END_REF]. Finally, ζ can be minimized by using a least square method.

It may be noted that J NL is the Jacobian matrix of the periodic nonlinear forces and does not depend on time. The eigenvalues of J are clearly related to the evolution of the nonlinear periodic solution Y e (t). Indeed, the real part of the corresponding unstable mode becomes equal to zero while the other real parts are negative when the computed solution has reached dynamical steady-state.

Hence the unstable real part is used as an extra equation for the root finding algorithm. In such a case, it will converge towards the steady-state solution where the dynamical equation and the real part are minimized. In conclusion, the final set of equations that has to be minimized is given by the two following functions f 1 Ũnl , ω and f 2 Ũnl , ω

f 1 Ũnl , ω = Λ eq (ω) Ũnl + Fnl Ũnl , ω -Feq < ǫ 1 (46) f 2 Ũnl , ω = |Re (λ) | < ǫ 2 (47)
where λ defines the eigenvalue of J that has the maximum real part. ǫ 1 and ǫ 2 are chosen residual coefficients.

The complete procedure and description of the Constrained Harmonic Balance Method is given in Figure 3.

The optimized additional initial conditions

As explained previously, the unknown parameters that have to be determined are the Fourier coefficients Ũnl and the frequency ω of the stationary periodic signal.

Firstly, when employing the static equilibrium position as the initial condition, computation can be very difficult and expensive. Hence they are too far from the final stationary nonlinear dynamical solution and so another initial estimation must be determined to save time and improve the computation procedure.

In this part of the paper, optimized additional initial conditions based on the complex nonlinear modal analysis [START_REF] Sinou | Stability analysis and non-linear behaviour of structural systems using the complex non-linear analysis (cnlma)[END_REF] are introduced and discussed. As explained previously by Sinou et al. [START_REF] Sinou | Stability analysis and non-linear behaviour of structural systems using the complex non-linear analysis (cnlma)[END_REF], starting from the hypothesis that the nonlinear unstable mode drives the dynamical solution, the evolution of the approximated solution curve, defined by considering only the contribution of the unstable mode, is given by:

Y 0 (t, p, λ) = p Ψe λt + Ψe λt (48)
where Ψ defines the nonlinear unstable mode and Ψ its conjugate. λ is the eigenvalue that corresponds to the unstable mode and p is an arbitrary chosen coefficient.

The optimized additional initial conditions are defined as the decomposition into Fourier coefficients of the previous expression of Y 0 (t, p, λ). In this case, Ψ is the eigenvector of the nonlinear unstable mode that has been obtained from the stability analysis. These optimized initial conditions work for quite a wide range of p and lead to the convergence of the harmonic balance method for the first calculation. Afterwards it is easy to compute solutions for given sets of parameters starting from previous results. Secondly, when using the harmonic balance method it is necessary to know the frequency of the periodic signal since the dynamical matrices are frequency-dependent and thus convergence may be laborious or even impossible if the chosen frequency is too approximative. The initial frequency value is selected to be equal to the unstable mode frequency computed by stability analysis. Figure 3 displays the algorithm procedure of the Constrained Harmonic Balance Method.

Results

In this part of the paper, the effectiveness of the Constrained Harmonic Balance Method will be illustrated for the nonlinear brake system presented in Section 2.

Non-linear stationary solutions

Limit cycles are computed for µ = µ 0 for cases with 1, 2, 3 and 10 harmonics. Results are compared with those obtained by applying temporal integration. The results shown in Figures 4 (a-b) correspond to a physical interface node where previous static reductions have been performed together with a modal displacement (c-d). First of all, the case with only one harmonic is discussed. Considering Figure 4, it clearly appears that the approximated solution obtained with one harmonic does not converge exactly with the final temporal solution. However, this solution gives a initial approximation of the limit cycle amplitudes, indicating that the first harmonic is one of the most significant components for the complete nonlinear solution. Hence estimations of the approximated nonlinear solutions computed with 2 or 3 harmonics are considered. The nonlinear limit cycles show a good fit with the temporal integration results in Figure 4. Moreover, the computed unknown frequency from the CHBM is very close to that of the temporal integration since the difference is less than 0.03% for 1,2, 3 and 10 harmonics. To make understanding easier for the reader, the evolutions of the residue of the nonlinear equation are given in Figure 7 for each case. Computational tests have been performed by assuming that the frequency of the limit cycles is different from the frequency obtained from the stability analysis. For the reader comprehension, neglecting changes in frequency of the self-excited nonlinear vibrations does not allow good estimation of the limit cycle amplitudes even if the difference between the frequency obtained from the stability analysis and the fundamental frequency of the non-linear oscillations appears to be very c-d), it is obvious that more complex solutions are better estimated when augmenting the number of harmonics. Slight differences appear between 2 and 3 harmonics because the third order becomes no more negligible and has to be taken into account for matching curves from the numerical integration well. This is even true for higher nonlinearities which generally involve coupling between Fourier coefficients and thus higher order responses. All the limit cycles computed with 10 harmonics match numerical integration correctly. Table 1 summarizes the results and relative errors for the three cases studied. Finally, Figure 5 shows the power spectrum ratios of the limit cycles computed at µ = µ 0 for each harmonic. Computation is done by summing power of each dof for a given harmonic: where P j is the power of the j th harmonic, a i and b i are the cosine and the sine coefficients of the i th dofs.

P j = 1 2 N dof i=1 a 2 i,j + b 2 i,j (49)
Then, the power ratio is computed. This consists in dividing each power harmonic by the total power of limit cycles computed by the temporal integration:

R j = P j P temporal ( 50
)
where R j defines the power ratio of the j th harmonic and P j is the power of the j th harmonic of the CHBM and P temporal is the total power of the temporal integration solution. Most of the energy appears to be concentrated in the first terms of the Fourier series. Table 1 displays the mean error between the numerical integration and the CHBM for the four different cases. As expected, the error decreases as the number of harmonics increases. Calculations with 10 harmonics have been carried out to ensure the convergence of computed solutions towards numerical integration and negligible differences have been detected. Moreover, limit cycles with 3 and 10 harmonics are close and the mean error of both methods is almost equal. Even when using only 3 harmonics for computing, the steady-state solution seems to be adequate in the present study. It can be seen that it is easy to consider a large number of harmonics with the proposed method. Furthermore, considering the limit cycle frequency as an unknown is fundamental in CHBM, otherwise the computation of the dynamical solution fails.

Convergence and Time Computation

This section concerns the convergence of the Constrained Harmonic Balance Method. Figure 6 displays the evolutions of real parts during the optimization computations. Whatever the number of harmonics, the real part converges to zero at the end of optimization, indicating the nonlinear stationary self-excited vibration of the brake system. Figure 7 shows the residue norm for every iteration. Although using several harmonics generally involves many iterations, only 11 iterations were needed for 2 and 3 harmonics in this case. Using 10 harmonics requires the convergence of 5 extra harmonics. Computation with 1 harmonic is not significant since the final approximated solution does not exactly match with the exact nonlinear solution. Thus it appears that the algorithm has difficulties in finding a root and seeks a convergence path that is costly in terms of iterations. Table 2 displays the computation time needed by temporal integration and CHBM to reach the dynamical solution. Both methods were run on MATLAB.

Even though the number of iterations for both 2 and 3 harmonics is equal, computation time is somewhat different with respectively 485 and 753 seconds required. This is obviously due to the fact that more harmonics mean more unknowns and thus more Jacobian evaluations in the optimization process. Computation time can be improved by using continuation methods, but this is not the aim of the paper. Using 10 harmonics is costlier and needs about 1 hour CPU time. However, CHBM computation time is still very low compared to the 200 hours needed for temporal integration. Moreover, massive disc storage is needed for temporal integration with more than one gigabyte of stored data (each dof is computed at a b each time t and stored on the disc) compared to a few kilobytes used by the Fourier coefficients of the CHBM. In the following part of this paper 3 harmonics will be used, regarding noticeably low relative errors both on displacement and on velocity. Moreover it offers a good compromise between accuracy and computation time. In conclusion, the proposed Constrained Harmonic Balance Method is well designed for a self-excited system because the results are accurate regarding the temporal approach and it is cheaper in terms of time consumption and disc storage.

Parametric studies : Interest of the CHBM

In this section of the paper, parametric studies will be undertaken for both the stability analysis and the limit cycles amplitudes.

Friction coefficient

The friction coefficient is generally considered as one of the most important parameters in brake systems. Good brake performances often signify a high friction coefficient yielding a high sound pressure level. Limit cycles are computed for µ = µ 0 , µ = 1.2µ 0 , µ = 1.4µ 0 , µ = 1.6µ 0 and µ = 2µ 0 and displayed In the rest of the study, the real parts as well as the frequencies computed and displayed in tables are normalized in relation to those in the nominal model (i.e. P 0 , k s , µ 0 and D1 = D2 = 1). To demonstrate the interest of considering the frequency as an unknown in the nonlinear method proposed in this paper, Table 3 gives the difference ∆f between the initial frequency of the unstable mode that has been obtained via the stability analysis for the nominal parameters and the final frequency of the self-excited vibration that has been obtained via the nonlinear method.

Considering Figures 8, it clearly appears that increasing the friction coefficient involves the higher vibration amplitudes for both pad and disc. Moreover, evolutions of the equilibrium point are observed, as indicated in Figure 8(a).

It can be seen that the Constrained Harmonic Balance Method allows determining the limit cycle amplitudes not only in the vicinity of the Hopf bifurcation point but also far from the Hopf bifurcation point µ 0 . In Table 3, a difference between the frequency of the limit cycle amplitudes and the frequency of the unstable mode can be observed.

Piston pressure

A variation in piston pressure has an effect on the pressure distribution at the disc/pad interface and the stability analysis of the brake system may be affected. Eigenvalues are computed for three different pressures 0.8P 0 , P 0 and 1.2P 0 where P 0 is the operational piston pressure. The evolutions of the real parts and the coalescences of the unstable and stable modes are illustrated in Figures 9(a-b). Table 4 gives the evolution of the Hopf bifurcation point. It appears that the pressure has an important effect of the stability. Basically, a higher piston pressure increases the degree of instability by moving the Hopf bifurcation point towards lower values without modifying the pattern shapes of the evolution of the real part or frequency coalescences. However, the frequency coalescence point is affected by a change in piston pressure, typically a high piston pressure results in a higher coupled frequency. Figures 9(c-d) illustrate the limit cycle amplitudes for the pad and the disc. To facilitate comprehension, the nonlinear vibrations are obtained at a fixed normalized friction coefficient µ = 1.05µ 0 . Hence pad interface deformation is equivalent for the three cases although the static positions are more affected for nodes at the pad/piston interface (see Figures 9(c)). This fact can be clearly explained by considering that higher piston pressure results in greater pad compression. The piston/pad interface is moved towards pad/disc interface. Disc deformation interface (Figure 9 4: Hopf bifurcation points and frequencies for variable piston pressure vibration amplitudes as the piston pressure increases.

Contact Stiffness

This part is devoted to the analysis of contact stiffness applied to the problem of brake squeal phenomenon. This parameter is very dependent on contact body stiffness and contact surface shapes. Since the disc is about a thousand times stiffer than the pad, contact stiffness usually depends on frictional materials. Both analyses are performed for three contact stiffnesses corresponding to a variation of the pad friction material' properties. For the sake of simplicity, k s is considered as a function of linear k l and nonlinear k nl contact stiffness springs, k s = f (k l , k nl ). The evolution of the real parts and frequencies of the stable and unstable modes are shown in Figures 10(a-b). Table 5 gives the evolution of the Hopf bifurcation point. As for the piston pressure case, an increase in the contact stiffness destabilizes the nonlinear system by decreasing the Hopf bifurcation point. Nevertheless, the patterns look alike and only a translation is observed for varying k s . Frequency lock-in is changed and a higher contact stiffness results in augmenting the coalescence frequency. This is quite logical since high contact stiffness tends to rigidify the whole system and thus increase the resonance frequencies.

When looking at limit cycles on Figure 10(c-d) which are computed for µ = 1.03µ 0 , it can be seen that increasing the contact stiffness decreases the nonlinear vibrations. Nevertheless, amplitudes do not change with the same ratio compared to the variable stiffness. A 30% increase in stiffness only changes the amplitudes by a factor of 1.1 while a 30% drop in contact stiffness results in lowering the limit cycle amplitudes by a mean factor of 3.3. Moreover, the static position of the equilibrium point can change drastically, as illustrated in Figure 10(c). Hence contact stiffness seems to be a key parameter in reducing brake squeal noise; nevertheless it has to satisfy other specifications which concern brake efficiency that are considered as the most important factors.

Damping

Damping seems to be a major parameter in reducing disc brake squeal. Squeal problems are resolved by applying damping shims on back plates, which works well sometimes. However, damping influences self-excited systems in which the coalescence of two modes is not yet fully understood and many studies have been performed to study this phenomenon. Hoffmann and Gaul [START_REF] Hoffmann | Effects of damping on mode-coupling instability in friction induced oscillations[END_REF] performed a stability analysis of a two dof system and showed that adding damping without precaution can lead to a paradoxical effect i.e. it can destabilize the system. Sinou and Jézéquel [START_REF] Sinou | The influence of damping on the limit cycles for a self-exciting mechanism[END_REF] studied the impact of modal damping on both stable and nonlinear systems and found that non-equally damped modes lead to a destabilization of the system and increase the amplitudes of limit cycles. Shin [START_REF] Shin | Analysis of disc brake noise using a two-degree-of-freedom model[END_REF] worked on a two dof system representing both pad and disc modes and shown that the equally damped case stabilizes the system by diminishing limit cycles when higher damping is applied but this is no longer the case when adding damping on only one mode. While the amplitude of the more highly damped system decreases, the amplitude of the other one increases. More recently, Fritz et al. [START_REF] Fritz | Effects of damping on brake squeal coalescence patterns-application on a finite element model[END_REF][START_REF] Fritz | Investigation of the relationship between damping and mode-coupling patterns in case of brake squeal[END_REF] performed stability analyses on a complete finite element model brake system. He confirmed that the ratio of the damping of the two modes involved in squeal is an essential key for controlling the stability of systems. For equally damped modes, the stability curves are lower and thus instability occurs for a higher friction coefficient. Nevertheless, in the case of large non-equally damped modes a smoothing effect occurs and pushes the Hopf bifurcation point towards lower values, thus instability appears for a lower friction coefficient compared to the equally damped case.

In the following, we investigate the effects of modal damping on stability and its impact on limit cycle amplitudes. Both the cases of equally and non-equally damped modes will be considered. Firstly, Figures 11(a-b) illustrate the evolutions of the real parts and the frequencies of the unstable and stable equally damped modes. As explained previously by Fritz et al. [START_REF] Fritz | Effects of damping on brake squeal coalescence patterns-application on a finite element model[END_REF][START_REF] Fritz | Investigation of the relationship between damping and mode-coupling patterns in case of brake squeal[END_REF], increasing equal modal damping causes a lowering effect and by consequence a stabilizing effect, by lowering the branches of the real part in the stable area. The critical Hopf bifurcation point is moved towards higher values for higher equal modal damping. Evolutions of the Hopf bifurcation point are given in Table 6. Moreover, it appears that the frequency lock-in phenomenon remains identical in the three cases. Mode 1 turns out to be the stable mode while Mode 2 is the unstable one. Secondly, Figures 12 (a)-(b) display nonlinear behavior and the limit cycle amplitudes for one DOF of the pad and one DOF of the disc computed for µ = 1.07µ 0 . The effects of damping appear to be very complex. Not only limit cycle amplitudes but also changes in the static position of the limit cycles are Case D1 D2 D1/D2 µ c /µ H ∆f (Hz) ∆f (Hz) at µ = 1.4µ 0 12(a). Surprisingly, the highly damped case does not necessarily induce low vibration amplitudes. Although the vibration amplitudes of the pad are slightly higher for the lower damped case, the highest dynamical response of the disc is found for the highest damped case, with an amplitude ratio of almost 4 in relation to the lowest damped case.

Figures 12 (c)-(d) show limit cycles computed by increasing the friction coefficient (µ = 1.4µ 0 in the current case). As explained previously in section 6.1, a higher friction coefficient involves higher limit cycle amplitudes. It can be seen that the effects of equally damped modes cannot be neglected and that the combined effect of the friction coefficient and damping is not trivial. For example, it appears that the influence of damping is weaker for µ = 1.4µ 0 compared to the case in the vicinity of the Hopf bifurcation point at µ = 1.07µ 0 . For the disc, the amplitude is still highest for the larger damped case, but the amplitude ratio in relation to the lowest damped case is less significant, with a value of 1.3. Moreover, it can be noted that the nonlinear amplitudes of the limit cycles do not follow proportionally the growth rate of the positive real part and the commonly held belief that the added damping would result in lower vibrations is not necessarily true.

Non-equally damped modes

Now we investigate the influence of non-equally damped modes for two cases that are D1/D2 = 0.5 and D1/D2 = 0.2. Figures 11(c-d) illustrate the associated evolutions of the real parts and frequencies of the unstable and stable equally damped modes. The reference remains the equally damped case where D1/D2 = 1. It should be noted that the damping of unstable Mode 2 turns out to be higher. The lowering effect due to high damping remains but its stabilizing effect is counterbalanced by the wellknown smoothing effect occurring in the vicinity of the Hopf bifurcation point (see Figure 11(c)), as mentioned by [START_REF] Hoffmann | Effects of damping on mode-coupling instability in friction induced oscillations[END_REF][START_REF] Sinou | The influence of damping on the limit cycles for a self-exciting mechanism[END_REF][START_REF] Shin | Analysis of disc brake noise using a two-degree-of-freedom model[END_REF][START_REF] Fritz | Effects of damping on brake squeal coalescence patterns-application on a finite element model[END_REF][START_REF] Fritz | Investigation of the relationship between damping and mode-coupling patterns in case of brake squeal[END_REF]. The real part branches of non-equally damped modes split with a smoother slope and become positive at a lower friction coefficient than the real part branch of the equally damped mode. This effect is stronger for higher asymmetrical modal damping cases, for example when D1/D2 = 0.2. Figure 13 (a-b) illustrates the nonlinear limit cycles that are computed at µ = µ 0 . The largest limit cycle amplitudes appear for the lowest damping ratio, D1/D2 = 0.2 although the unstable mode is the more highly damped one. More complex behavior is found for pad deformation when D1/D2 = 0.2. With a higher friction coefficient, the real part curves cross each other at around µ = 1.03µ 0 and case D1/D2 = 1 has the larger positive real part, contrary to the D1/D2 = 0.2 case which becomes the most stable beyond µ = 1.03µ 0 . Figure 13 (c-d Figure 11: Evolution of real parts and frequency coalescences (a-b) equally damped: -1 st case, --2 nd case, ... 3 rd case (c-d) non-equally damped: -1 st case, --6 th case, ... 7 th case Large pad amplitudes are obtained for low damping ratios and are considered as the most stable cases at µ = 1.4µ 0 by the stability analysis. Nevertheless, the highest amplitudes for the disc are derived from the most unstable case, which is the equally damped one. It can be clearly observed that it is not possible to establish a link between the values of the real parts and corresponding vibrating states since they can be higher or lower depending on the different effects of the physical parameters on both stability and nonlinear behavior. Moreover, it is noted that the static equilibrium point changes with the variation of non-equally damped modes, as indicated in Figure 13 (a).

To further investigate the influence of modal damping, we propose to invert damping ratios D1/D2 = 2, D1/D2 = 5 and perform a stability analysis and determine the nonlinear limit cycles. It should be borne in mind henceforth that stable Mode 1 is the more highly damped one, whereas the damping of unstable Mode 2 is decreased. The evolutions of the real parts and frequencies of the stable and unstable modes are similar to the previous case thus the conclusions on stability are identical. Figure 14 (a-b) and (c-d) illustrates the nonlinear limit cycles for µ = 1.07µ 0 and µ = 1.4µ 0 , respec- tively. In the of the Hopf bifurcation point (i.e. µ = 1.07µ 0 ), the highest amplitude is still obtained for the largest difference in damping, i.e. D1/D2 = 5, but with an amplitude ratio of almost 57 compared to an amplitude ratio of 37 for the preceding case. When considering a case far from the Hopf bifurcation point (at µ = 1.4µ 0 for example), the evolution of the limit cycles are complex and it appears more difficult to give a general rule on the effects of damping. For example, although the real part of the unstable mode for the more asymmetrical damping case (i.e. D1/D2 = 5) is the lowest one, this configuration involves the higher vibrating state for both pad and disc. It must be compared to the previous case study where the more asymmetrical damping case (i.e. D1/D2 = 0.2) exhibited higher amplitudes only for the pad. Here, the unstable mode is the more weakly damped one and even if the corresponding real part is below that of case 1, the corresponding amplitudes are higher.

Inverting modal damping has a weak influence in the vicinity of the Hopf bifurcation point, as seen in Figure 13 (a-b) and Figure 14 the friction coefficient. When the unstable mode is the more strongly damped one, the limit cycles are smaller than the equally damped case by a factor of 2 (Figure 13 (d)), but the conclusions are totally different when inverting the damping ratio: The corresponding limit cycles are about 3 times larger than for the equally damped case (Figure 14 (d)) .

This example shows that inverting modal damping distribution has considerable effects not only on the stability of the system but also on the nonlinear amplitudes of the limit cycles. For example, case D1/D2 = 0.2 displays the smallest disc oscillations while case D1/D2 = 5 has the largest disc oscillations and the conclusions for finding the best model are totally different. Considering the previous results, it appears that structural damping is a key factor when dealing with nonlinear autonomous systems, but nevertheless it is a complex phenomenon and it has to be considered with care to ensure good silent brake system design. Not only the quantity but also the distribution of damping have to be taken into account thoroughly to avoid unexpected results.

Conclusion

In this paper proposed a novel nonlinear method called the Constrained Harmonic Balance Method. This original approach allows the determination of the stationary nonlinear periodic solution of a nonlinear mechanical system subject to flutter instability, by the addition of an extra-constraint in the classical Harmonic Balance Method. This additional constraint allows eliminating the static equilibrium point (i.e. the trivial static solution of the nonlinear problem that would be obtained by applying the classical Harmonic Balance Method)) and gives only the stationary nonlinear oscillations. Moreover, the frequency is added as an unknown since the frequency of a self-excited system is not known a priori and may change for varying parameters. Also, the dynamical solution cannot be computed if only using the frequency resulting from the stability analysis. An application to disc brake squeal was performed to illustrate the effectiveness of the nonlinear method.

Numerical results correlated well with a classical time domain algorithm in terms of both amplitude and frequency. The results of the CHBM is highly dependent on the number of harmonics. A power ratio computation shows that the major part of the energy is concentrated in the first harmonic, but retaining only the latter does not lead to the steady-state solution. In order to adapt to the complex behaviors of the solutions, more harmonics are required in the Fourier series. The computation time of the new method is very short compared to that of a classical temporal integration algorithm and thus is well designed for intensive computation in the case of parameter-dependent systems. The effectiveness of this method's application to a disc brake system is emphasized in the last part of the paper, which describes the parametric studies performed. Fast limit cycle computations were achieved for a large number of operational parameters and conclusions were obtained. The complementarity between the stability analysis and the complex nonlinear vibrational behavior appears to be essential for carrying out a complete design study of a brake system. Moreover, it was shown that not only the friction coefficient, but also piston pressure, nonlinear stiffness and structural damping are important factors to take into account to avoid poor design.

Figure 1 :

 1 Figure 1: Finite element models of the brake system (a) Pad (b) Disc

Figure 2 :

 2 Figure 2: Stability analysis (a) Evolution of the real parts of the stable and unstable modes (b) Coalescence of the two corresponding eigenvalues

Figure 3 :

 3 Figure 3: Algorithm procedure

1 )Figure 4 :

 14 Figure 4: Limit cycles using the classical temporal integration and the modified HBM at µ 0 ; -temporal integration, -1 harmonic, --2 harmonics, -.-3 harmonics, ... 10 harmonics; (b,d) zoom

Figure 5 :

 5 Figure 5: Power ratio for the first 10 harmonics

Figure 6 :

 6 Figure 6: Evolution of the real part; -1 harmonic, --2 harmonics, -.-3 harmonics, ... 10 harmonics (a) all cases (b) zoom on 2, 3 and 10 harmonics

Figure 7 :

 7 Figure 7: Evolution of residues during optimization; -1 harmonic, --2 harmonics, -.-3 harmonics, ... 10 harmonics (a) all cases (b) zoom on 2, 3 and 10 harmonics

Figure 8 :

 8 Figure 8: Limit cycles for pad and disc nodes; -µ 0 , -1.2µ 0 , --1.4µ 0 , ... 1.6µ 0 , -.-. 2µ 0 (a) Pad, (b) Disc

Figure 9 :

 9 Figure 9: Influence of piston pressure; -P 0 , --0.8P 0 , ... 1.2P 0 (a) Evolution of real parts (b) Frequency coalescences; Limit cycles for (c) Pad, (d) Disc

1 )Figure 10 :

 110 Figure 10: Influence of contact stiffness; -k s , --0.7k s , ... 1.3k s (a) Evolution of real parts (b) Frequency coalescences; Limit cycles for (c) Pad, (d) Disc

Figure 12 :

 12 Figure 12: Limit cycles for equally damped modes and plotted for two different µ (a-b) Limit cycles plotted at µ = 1.07µ 0 -1 st case, --2 nd case, ... 3 rd case (c-d) Limit cycles plotted at µ = 1.4µ 0 -1 st case, --2 nd case, ... 3 rd case; (a-c) Pad, (b-d) Disc

Figure 13 :

 13 Figure 13: Limit cycles for non-equally damped modes and plotted for two different µ (a-b) Limit cycles plotted at µ = µ 0 -1 st case, --4 th case, ... 5 th case (c-d) Limit cycles plotted at µ = 1.4µ 0 -1 st case, --4 th case, ... 5 th case; (a-c) Pad, (b-d) Disc

1 )Figure 14 :

 114 Figure 14: Limit cycles for non-equally damped modes and plotted for two different µ (a-b)Limit cycles plotted at µ = µ 0 -1 st case, --6 th case, ... 7 th case (c-d) Limit cycles plotted at µ = 1.4µ 0 ; -1 st case, --6 th case, ... 7 th case; (a-c) Pad (b-d) Disc

Table 1 :

 1 Relative error for 1, 2, 3 and 10 harmonics

	1 H 2 H 3 H 10 H

Table 3 :

 3 Hopf bifurcation points and frequencies for variable friction coefficient

Table 5 :

 5 Hopf bifurcation points and frequencies for variable stiffness

	4µ 0

Table 6 :

 6 Hopf bifurcation points and frequencies for damping parameters observed, as indicated in Figure

		1	1	1	1	0.21	0.72
	2	5	5	1	1.01	0.18	0.70
	3	10 10	1	1.07	0.10	0.63
	4	1	2	0.5	0.97	0.68	1.28
	5	1	5	0.2	0.84	1.64	1.98
	6	2	1	2	0.96	-0.50	0.24
	7	5	1	5	0.84	-0.87	-0.18

  ) presents the limit cycles of both the nodes considered previously, computed at µ = 1.4µ 0 .

	a										b
		0.8										1.0015
		0.6
												1.001
		0.4
	Normalized Real Part	-0.2 0 0.2									Normalized Frequency	0.9995 1 1.0005
		-0.4
												0.999
		-0.6
		0.7 -0.8	0.8	0.9	1	1.1	1.2	1.3	1.4	1.5	1.6	0.7 0.9985	0.8	0.9	1	1.1	1.2	1.3	1.4	1.5	1.6
					Normalized Friction Coefficient						Normalized Friction Coefficient
	c										d
		0.8										1.0015
		0.6
												1.001
		0.4
	Normalized Real Part	-0.2 0 0.2									Normalized Frequency	0.9995 1 1.0005
		-0.4
												0.999
		-0.6
		0.7 -0.8	0.8	0.9	1	1.1	1.2	1.3	1.4	1.5	1.6	0.7 0.9985	0.8	0.9	1	1.1	1.2	1.3	1.4	1.5	1.6
					Normalized Friction Coefficient						Normalized Friction Coefficient

Acknowledge

Jean-Jacques Sinou gratefully acknowledges the financial support of the french National Research Agency through the program Young researcher ANR-07-JCJC-0059-01-CSD 2.