Parameter Determination of a Minimal Model for Brake Squeal

In the research into the mechanism of brake squeal, minimal models with two degrees of freedom (DoFs) are widely used. Compared with the finite element method, the minimal model is more concise and efficient, making it easier to analyze the effect of parameters. However, how to accurately determine its kinetic parameters is rarely reported in the literature. In this paper, firstly, the finite element model of a disc brake is established and the complex eigenvalue analysis (CEA) is carried out to obtain unstable modes of the brake. Then, an unstable mode with seven nodal diameters predicted by CEA is taken as an example to establish the 2-DoF model. In order that the natural frequency, Hopf bifurcation point and real parts of eigenvalues of the minimal model coincide with that of the unstable mode with seven nodal diameters, the response surface method (RSM) is applied to determine the kinetic parameters of the minimal model. Finally, the parameter-optimized minimal model is achieved. Furthermore, the negative slope of friction-velocity characteristic is introduced into the model, and transient analysis (TA) is used to study the effect of braking velocity on stability of the brake system. The results show that the brake system becomes unstable when braking velocity is lower than a critical value. The lower the velocity is, the worse the stability appears, and the higher the brake squeal propensity is.


Introduction
Squeal noise of brake system is one of the difficult problems that need to be solved by automobile manufacturers.It is generally believed that the disc brake squeal is self-excited vibrations induced by friction forces.However, for the mechanisms of self-excited vibrations, the existing literatures proposed a different hypothesis.The first type of literaturesholds that the variation characteristics of the friction coefficient lead to squeal instability, including mechanisms like stick-slip [1][2][3] and negative friction-velocity relationship [4,5].Nevertheless, the stick-slip mechanism is not sufficient to explain the occurrence of the squeal [6].As pointed out by Ouyang et al. [5], stick-slip occurs only at low speeds of a mass driven across a dry surface, while at a sufficiently high speed the mass will be permanently sliding.This proves that the instability of brake discs can be induced by a negative slope of friction-velocity characteristic even without stick-slip motion.The second type of literature indicates that self-excited vibrations can be induced even at a constant friction coefficient, and demonstrates that the generation of squeal is related to the dynamics of the geometric structure [7,8].The sprag-slip model introduced in early studies by Spurr [7] illustrates the mechanism of squeal in terms of structural geometry, which means that tribological characteristics should not be considered as the only reason for squeal instability [9].Later, many scholars confirmed the existence of mode-coupling phenomenon [10][11][12][13] in brake squeal.Hoffman et al. [11] found that when the friction coefficient is greater than a critical value, namely the Hopf bifurcation point, the oscillation frequencies of two structural modes gradually come closer until they coalesce into a pair of an unstable mode and a stable mode.It is clear that mode coupling explains the squeal phenomenon under a constant friction coefficient instead of a variable friction coefficient.In fact, when mode coupling occurs, the system not only generates unstable oscillations, but simultaneously accompanied by a notable nonlinear variation of friction coefficient.This has an important effect on the dynamic response of the system.Therefore, the stability analysis in this paper is based on the combination of mode coupling and negative slope of friction-velocity characteristics.
The modeling method for brake squeal generally includes two categories.One category is to establish the finite element model of brake, then complex eigenvalue analysis (CEA) or transient analysis (TA) is performed for the finite element model [14][15][16].Ouyang et al. [17] and Kinkaid et al. [18] give a detailed summary of the application of CEA and TA in the study of brake squeal.The other category to investigate squeal is based on minimal models such as the abstract 2-DoF model.Hoffman et al. [11] proposed a 2-DoF model to explain self-excited vibrations induced by mode-coupling mechanism.Shin et al. [19] adopted a 2-DoF model containing lumped masses representing the disc and the pad respectively, through which the effect of negative slope of friction-velocity characteristics on the oscillation instability was investigated.Popp et al. [20] considered the torsional mode of the brake disc in their 2-DoF model.Later, based on the three models above, Wagner [21] put forward a 2-DoF model that further considers vibration of the disc as a rotating elastic plate, so as to associate easily with real disk brake [22].However, the aforementioned literatures have not given detail descriptions about parameters determination method and its process of the model, and the stability of the model depends greatly on the value of parameters.
Since brake components are geometrically complex and are essentially continuous media with infinite DoFs, it is appropriate to use the finite element method with less model simplification and a large number of DoFs [17].The CEA of the finite element method can simultaneously obtain several unstable modes, but it is limited by underestimate or overestimate error due to a linearization assumption [23].Compared with eigenvalue computation, the TA of the finite element method can achieve the dynamic response of components with nonlinear effects included, such as nonlinearity of friction, contact pressure, and lifting off of the pad (e.g., contact separation between the two sliding surfaces), which are referred to as deterministic chaos by Oberst et al. [24].So additional unstable modes may be detected while CEA fails.However, it is not conducive to extensive parameter designs because of excessive time consumption [25], so its engineering application is limited.In contrast with a full-scale brake finite element model, the 2-DoF minimal model is of high efficiency during stability analysis owing to its fewer DoFs, although it is a simplification, and can consider only a certain order mode.It is helpful to analyze the effect of key parameters on brake squeal, and to reveal the squeal mechanism [19].Nevertheless, the parameters in the model are difficult to determine.
Therefore, this paper aims to explore a more accurate method for parameter determination of the 2-DoF minimal model, so as to make up for the shortcomings of existing literatures.The response surface method (RSM) is applied in the parameter determination.Then, TA is performed for the parameter-optimized minimal model to discuss the effect of braking velocity on squeal stability.The results demonstrate that the squeal is more prone to occur under relatively low velocity.

Finite Element Analysis of Brake Squeal
As shown in Figure 1, a simple finite element brake model with a disc, two pads and two back plates is established in ABAQUS.Table 1 shows the material properties of the three brake components [26], where the brake pads are made of orthotropic material, and its nine elastic constants are given.The back plates are constrained using YASYMM boundary condition, namely, constraining the X, Z direction of the displacement and rotation around the Y axis, and the bolt holes on the disc are constrained with PINNED boundary condition, namely, constraining the X, Y, Z direction of the displacement.A uniform pressure of 1 MPa is applied to the back plates to simulate the actual braking hydraulic load.The complex modal analysis of the model is carried out, and the rotation effect is applied to the brake disc in the input file.Finally, eigenvalues of the first 55 complex modes of the model are extracted.Figure 2 shows the frequencies and real parts of the eigenvalues of unstable modes.It can be seen from Figure 2 that the brake system has several unstable modes which may lead to brake squeal.With the increase of the friction coefficient, the number of unstable modes increases; therefore, the propensity of brake squeal increases.For the friction coefficient 9 .0 = μ , 8 unstable The complex modal analysis of the model is carried out, and the rotation effect is applied to the brake disc in the input file.Finally, eigenvalues of the first 55 complex modes of the model are extracted.Figure 2 shows the frequencies and real parts of the eigenvalues of unstable modes.
Appl.Sci.2018, 8, 37 3 of 15 the X, Z direction of the displacement and rotation around the Y axis, and the bolt holes on the disc are constrained with PINNED boundary condition, namely, constraining the X, Y, Z direction of the displacement.A uniform pressure of 1 MPa is applied to the back plates to simulate the actual braking hydraulic load.The complex modal analysis of the model is carried out, and the rotation effect is applied to the brake disc in the input file.Finally, eigenvalues of the first 55 complex modes of the model are extracted.Figure 2 shows the frequencies and real parts of the eigenvalues of unstable modes.It can be seen from Figure 2 that the brake system has several unstable modes which may lead to brake squeal.With the increase of the friction coefficient, the number of unstable modes increases; therefore, the propensity of brake squeal increases.For the friction coefficient 9 .0 = μ , 8 unstable It can be seen from Figure 2 that the brake system has several unstable modes which may lead to brake squeal.With the increase of the friction coefficient, the number of unstable modes increases; therefore, the propensity of brake squeal increases.For the friction coefficient µ = 0.9, 8 unstable modes are predicted in the system.Their frequencies are 4060.2Hz, 5525.5 Hz, 7213.0Hz, 8429.0Hz, 9042. 1 Hz,11,056 Hz,11,198 Hz and 11,496 Hz respectively.
For instance, the unstable mode with seven node diameters is shown in Figure 3. Figure 4a shows the evolution of real parts of the mode with respect to the friction coefficient from 0.2 to 0.9, namely the Hopf bifurcation diagram.The bifurcation point of this mode is observed for µ 0 = 0.4.For friction coefficients greater than 0.4, the real parts of this mode are positive, which signifies an unstable mode.Moreover, the maximum real part value is 713.41 for µ = 0.9. Figure 4b illustrates the evolution of frequencies of this mode versus the friction coefficient.With the increase of friction coefficient, two modal frequencies gradually get closer until coalescence occurs.The coalescence frequency at the bifurcation point is 7168.6Hz.In Figure 4c, a three-dimensional graph including real parts and frequencies of mode coupling is given, where the red line represents the unstable mode, and the blue line represents the stable mode.
Appl.Sci For instance, the unstable mode with seven node diameters is shown in Figure 3. Figure 4a shows the evolution of real parts of the mode with respect to the friction coefficient from 0.2 to 0.9, namely the Hopf bifurcation diagram.The bifurcation point of this mode is observed for 4 .0 0 = μ .For friction coefficients greater than 0.4, the real parts of this mode are positive, which signifies an unstable mode.Moreover, the maximum real part value is 713.41 for 9 .0 = μ .Figure 4b illustrates the evolution of frequencies of this mode versus the friction coefficient.With the increase of friction coefficient, two modal frequencies gradually get closer until coalescence occurs.The coalescence frequency at the bifurcation point is 7168.6Hz.In Figure 4c, a three-dimensional graph including real parts and frequencies of mode coupling is given, where the red line represents the unstable mode, and the blue line represents the stable mode.For instance, the unstable mode with seven node diameters is shown in Figure 3. Figure 4a shows the evolution of real parts of the mode with respect to the friction coefficient from 0.2 to 0.9, namely the Hopf bifurcation diagram.The bifurcation point of this mode is observed for   Then, CEA results of the unstable mode with seven modal diameters will be referred to as a standard in order to establish a minimal model with two DoFs for brake squeal.

2-DoF Minimal Model for Brake Squeal
Referring to [25], a 2-DoF mechanical model is established in Figure 4.In the model, two pads are symmetrically arranged on two sides of the brake disc and have the same motion.1 m denotes the generalized mass of the two pads, and it has a normal mode in Y direction.2 m denotes the generalized mass of the disc, and it has a torsional mode in X direction.To couple the normal mode with the torsional mode, an angle θ is considered between the friction surface and horizontal direction.v is the rotational velocity of the brake disc.Friction coefficient μ of the contact surface is assumed to be a constant.brake F is the braking force applied on the pads. 1 c and 2 c are the damping of the two modes respectively.1 k and 2 k are, respectively, the contact stiffness and torsional stiffness of the disc.Notice that the stiffness is simplified to be linear as an approximation; nevertheless, in reality it should be nonlinear [27].
The dynamic differential equation of the braking system in Figure 5 can be expressed as Then, CEA results of the unstable mode with seven modal diameters will be referred to as a standard in order to establish a minimal model with two DoFs for brake squeal.

2-DoF Minimal Model for Brake Squeal
Referring to [25], a 2-DoF mechanical model is established in Figure 4.In the model, two pads are symmetrically arranged on two sides of the brake disc and have the same motion.m 1 denotes the generalized mass of the two pads, and it has a normal mode in Y direction.m 2 denotes the generalized mass of the disc, and it has a torsional mode in X direction.To couple the normal mode with the torsional mode, an angle θ is considered between the friction surface and horizontal direction.v is the rotational velocity of the brake disc.Friction coefficient µ of the contact surface is assumed to be a constant.F brake is the braking force applied on the pads.c 1 and c 2 are the damping of the two modes respectively.k 1 and k 2 are, respectively, the contact stiffness and torsional stiffness of the disc.Notice that the stiffness is simplified to be linear as an approximation; nevertheless, in reality it should be nonlinear [27].
The dynamic differential equation of the braking system in Figure 5 can be expressed as where N and T are the normal force and Coulomb-type sliding friction force subjected by m 2 , respectively.and y is the vertical component of v. Substituting Equations ( 2) and (3) into Equation ( 1), the force N can be eliminated, so a reduction form with two-dimensional matrices of Equation ( 1) can be written as where the displacement vector is x = (X, Y) T , the mass, damping and stiffness matrices of the brake system are denoted by M, C, K respectively, and F is the braking force vector.Their expressions are deduced as follows: Appl.Sci.2018, 8, 37 6 of 15 where where N and T are the normal force and Coulomb-type sliding friction force subjected by 2 m , respectively.and y is the vertical component of v .Substituting Equations ( 2) and (3) into Equation ( 1), the force N can be eliminated, so a reduction form with two-dimensional matrices of Equation ( 1) can be written as where the displacement vector is , the mass, damping and stiffness matrices of the brake Assuming a small perturbation x = (X, Y) T at the equilibrium point x 0 = (X 0 , Y 0 ) T , which satisfies the following condition: Equation ( 4) is simplified as For convenience, a vector ϕ == (x, .

x)
T is introduced.Then Equation ( 10) can be transformed into state space equation: . ϕ + Aϕ = 0 (11) where The stability of the 2-DoF model can be analyzed by the eigenvalues of A, and the eigenvalues are expressed as where a i is the real part of eigenvalue λ i , and the imaginary part ω i is the frequency of the mode.According to the complex modal theory, the damping ratio of the model at this frequency is [28] Equation ( 14) indicates the propensity of squeal instability of the model.When at least one real part of the eigenvalue satisfies a i > 0, the system is of negative damping and hence unstable.Accordingly, the response becomes divergent with time increases.Conversely, when all real parts of eigenvalues satisfy a i ≤ 0, the system is stable.CEA predicts the stability of systems by the value of real parts.This method has high computational efficiency and can be used in the system with a large number of DoFs, so it is widely applied in engineering.In the above linear stability analysis a small perturbation x is assumed, while large amplitude perturbations happens, stick-slip can arise, which is not the subject of the paper.We refer interested readers to Papangelo et al. [29] for detailed deductions.Lumped models in [29,30] illustrates a bistable zone where steady sliding and stick-slip limit cycle coexists in a certain range of the control parameter.
The bifurcation point and real part are key factors to indicate the squeal propensity of the system.And the unstable modal frequency after modal coupling is regarded as squeal frequency by extensive studies.Therefore, in Figure 4, the coupling frequency 7168.6 Hz, the bifurcation point µ 0 = 0.4, and the maximum real part 713.41 at µ = 0.9 are chosen as three optimization objectives to optimize the kinetic parameters in the 2-DoF model, so that the vibration characteristics of the minimal model coincide with that of the selected finite element unstable mode.

Initialization of Parameters
The natural frequency of the minimal model is mainly affected by the mass and stiffness.First of all, initial values are chosen as m 1 = 1.2 kg, m 2 = 1 kg, θ = 0.15 rad, µ = 0.4.Assuming k 2 = 2.1×10 9 N/m, relative error between the coalescence frequency of the minimal model and the optimization objective, which is 7168.6Hz, are calculated.Figure 6 gives the evolution of error versus k 1 .It can be determined that when k 1 = 2.27×10 9 N/m, the minimum frequency error is only about 0.01%.The Hopf bifurcation point of the minimal model is mainly affected by its slope angle θ.When the kinetic parameters are set as same as the above-mentioned values, as shown in Figure 7, the value of the bifurcation point µ 0 decreases with the increase of θ.It can be determined that when θ = 0.0065 rad, the corresponding bifurcation point is µ 0 = 0.4, which is consistent with the finite element model result.The Hopf bifurcation point of the minimal model is mainly affected by its slope angle θ.When the kinetic parameters are set as same as the above-mentioned values, as shown in Figure 7, the value of the bifurcation point 0 μ decreases with the increase of θ.It can be determined that when θ = 0.0065 rad, the corresponding bifurcation point is , which is consistent with the finite element model result.and θ.

Response Surface Function
To find functional relationships between parameters 1 k , 2 k and the multiple optimization objectives, namely the frequency, real part and θ in the 2-DoF model, response surface methodology (RSM) is utilized for the optimization.The response surface is constructed by the following binary quadratic polynomial: Numerous calculation cases show that the stiffness k 1 , k 2 have the most significant effect on the real part of the model.As the stiffness getting greater, the real part of unstable mode becomes larger, and thus the squeal propensity increase.Variations of the stiffness also result in slight change of the coupling frequency and Hopf bifurcation point, which is mainly affected by angle θ.Therefore, the coupling frequency, bifurcation point and real part of the finite element mode is selected as targets.And multi-parameter and multi-objective optimization is performed to determine k 1 , k 2 and θ.

Response Surface Function
To find functional relationships between parameters k 1 , k 2 and the multiple optimization objectives, namely the frequency, real part and θ in the 2-DoF model, response surface methodology (RSM) is utilized for the optimization.The response surface is constructed by the following binary quadratic polynomial: where k 1 and k 2 are the variables, R represents the response value of the optimization objective, c = c 0 , c 1 , c 2 , c 3 , c 4 , c 5 ) T is a vector composed of 6 coefficients remaining to be estimated, ε is the statistical error between actual values and fitted values.The binary quadratic function has the advantages of simple form and good approximation to actual response.Vector c can be easily solved by the least square method.More than 6 sets of test data are needed for the calculation.Assuming n sets of data are selected, and n > 6.Then the n response surface equations can be written as where R is the n × 1 dimensional response vector, M k is a n × 6 dimensional matrix composed of test variables k 1 and k 2 : . . .
Then the least square solution of coefficient vector c is

Computational Design and Optimization Results
The steepest ascent design was applied to approach the response region of the optimization objectives.Table 2 gives the response results of real part versus k 1 and k 2 .As seen obviously that the response value of computation No. 3 is closest to the optimization objective.Therefore, the following central composite computation is designed with k 1 = 2.37×10 9 N/m, k 2 = 2.1×10 9 N/m as level 0, k 1 + δ, k 2 + δ as level 1 and k 1 − δ, k 2 − δ as level −1, where δ equals 0.03×10 9 N/m.The center composite computation involves 2 variables and 3 levels, so 9 sets of computations are needed.The results are listed out in Table 3.It should be noted that the response values of θ in each set are determined according to the method shown in Figure 7 because θ has a predominant influence on Hopf bifurcation point, whereas the response value of the real part and the coupling frequency is calculated by CEA, as mentioned in Section 2.2.
And the response surface plots are shown in Figure 8. Ignore the effect of fitting error ε , and let Equations ( 19) and ( 20) equal to the optimization objectives of finite element results, respectively.Thus, the following equations are obtained: Values of 1 k and 2 k which meet both optimization objectives of the frequency and real part can be solved by Equation (22).The value of θ can be further obtained from Equation (21).Finally, all kinetic parameters are determined as m1 = 1.2 kg, m2 = 1 kg, k1 = 2.3658 × 10 9 N/m, k2 = 2.0858 × 10 9 N/m, θ = 0.0018 rad.Under these optimized parameters, Figure 9a displays the evolutions of both real parts and frequencies versus friction coefficient.Moreover, a comparison between the CEA results of the minimal model and that of the finite element model is illustrated in Figure 9b,c, where the black line represents the finite element result previously given in Figure 4.It is shown that the Hopf bifurcation point of the 2-DoF model is . In Figure 9b, when the friction coefficient is greater than 0.4, the real parts become positive and continuously increase with the increase of the friction coefficient.
When the friction coefficient is 0.9, the corresponding maximum real part is 711.8.In Figure 9c,  Ignore the effect of fitting error ε, and let Equations ( 19) and ( 20) equal to the optimization objectives of finite element results, respectively.Thus, the following equations are obtained: Values of k 1 and k 2 which meet both optimization objectives of the frequency and real part can be solved by Equation (22).The value of θ can be further obtained from Equation (21).Finally, all kinetic parameters are determined as m 1 = 1.2 kg, m 2 = 1 kg, k 1 = 2.3658 × 10 9 N/m, k 2 = 2.0858 × 10 9 N/m, θ = 0.0018 rad.
Under these optimized parameters, Figure 9a displays the evolutions of both real parts and frequencies versus friction coefficient.Moreover, a comparison between the CEA results of the minimal model and that of the finite element model is illustrated in Figure 9b,c, where the black line represents the finite element result previously given in Figure 4.It is shown that the Hopf bifurcation point of the 2-DoF model is µ 0 = 0.4.In Figure 9b, when the friction coefficient is greater than 0.4, the real parts become positive and continuously increase with the increase of the friction coefficient.When the friction coefficient is 0.9, the corresponding maximum real part is 711.8.In Figure 9c, frequencies of the two close modes tend to approach each other gradually as the friction coefficient increases, and finally coupling occurs at the bifurcation point.The coupling frequency is 7168.6Hz.
After parameter optimization, the coupling frequency, Hopf bifurcation point and the maximum real part of the minimal model are all consistent with the selected finite element mode.The error is within acceptable limits.The optimized 2-DoF model can be more easily associated with the finite element braking model with a disc and two pads.To a certain extent, it represents the coupling properties of the unstable mode with seven nodal diameters.More reliable numerical results of braking stability can be obtained by using the parameters-optimized minimal model.where µ s is the static friction coefficient, µ k is the negative slope coefficient which represents the descent rate of friction coefficient, v is the rotation velocity of the disc which is a constant, .X/ cos θ is the velocity projection of m 2 in horizontal direction.
Subsisting Equation (23) into Equation (3) yields a more complex friction model as the static friction coefficient is introduced.Here a simplification is made for the friction law that the stick phase (e.g., the friction force is T < µ s N) is not considered, so that nonlinearity of discontinuity does not involve in the friction model, and the system can be easily integrated by a standard ODE-solver available in MATLAB [33].Then the time-history response of the equation is solved via the classic fourth-order Runge-Kutta algorithm.The initial conditions are set as X 0 = Y 0 = .X 0 = .Y 0 = 0, and µ s = 0.4625, µ k = 0.025, so that when the velocity of brake disc equals to 2.5 m/s, the initial dynamic friction coefficient is µ = 0.4, just at the bifurcation area.The effect of brake pressure on the stability is not the interest of this paper, and F brake = 100 N in this case.Figure 10a gives the velocity response over time in the Y direction at v = 2.5 m/s, and Figure 10b is the response with timescale zoomed in.Since the dynamic response in the X direction is similar, it is not shown.In Figure 10b, the steady state response in the Y direction is a periodic signal with changing amplitude.The reason for this phenomenon is known from Equation (23), that during oscillation process of the system, the relative velocity v − .
x is constantly changing.This causes the friction coefficient to change repeatedly around the bifurcation point; namely, the system is in constant conversion between stable and unstable.When the friction coefficient exceeds the bifurcation point, the system becomes unstable, so the velocity response increases.At the same time, the friction coefficient reduces below the bifurcation point as the relative velocity increases, then the vibration response decreases.
where s μ is the static friction coefficient, k μ is the negative slope coefficient which represents the descent rate of friction coefficient, v is the rotation velocity of the disc which is a constant, θ cos / X  is the velocity projection of 2 m in horizontal direction.
Subsisting Equation (23) into Equation (3) yields a more complex friction model as the static friction coefficient is introduced.Here a simplification is made for the friction law that the stick phase (e.g., the friction force is ) is not considered, so that nonlinearity of discontinuity does not involve in the friction model, and the system can be easily integrated by a standard ODE-solver available in MATLAB [33].Then the time-history response of the equation is solved via the classic fourth-order Runge-Kutta algorithm.The initial conditions are set as , and , so that when the velocity of brake disc equals to 2.5 m/s, the initial dynamic friction coefficient is 4 .0 = μ , just at the bifurcation area.The effect of brake pressure on the stability is not the interest of this paper, and N 100 brake = F in this case.Figure 10a gives the velocity response over time in the Y direction at m/s 5 .

= v
, and Figure 10b is the response with timescale zoomed in.
Since the dynamic response in the X direction is similar, it is not shown.In Figure 10b, the steady state response in the Y direction is a periodic signal with changing amplitude.The reason for this phenomenon is known from Equation (23), that during oscillation process of the system, the relative velocity is constantly changing.This causes the friction coefficient to change repeatedly around the bifurcation point; namely, the system is in constant conversion between stable and unstable.When the friction coefficient exceeds the bifurcation point, the system becomes unstable, so the velocity response increases.At the same time, the friction coefficient reduces below the bifurcation point as the relative velocity increases, then the vibration response decreases.In addition, to study the effect of different velocities on oscillations, the steady-state responses for 3.5 m/s, 2.5 m/s, 1.5 m/s and 0.5 m/s are calculated, respectively.The limit cycles are shown in Figure 11a,b.It can be seen that when v equals 3.5 m/s, the system is stable and the limit cycle amplitude is very small.When v equals 2.5 m/s, the system is unstable, and the limit cycle amplitude is greater than that of 1.5 m/s and 0.5 m/s.With the decrease of disc velocity, the limit cycle amplitude becomes larger, that is to say, the instability increases.However, the increase of the amplitude is not proportional to the decrease of v .When v continues to decrease to 0.5 m/s, the limit cycle amplitude In addition, to study the effect of different velocities on oscillations, the steady-state responses for 3.5 m/s, 2.5 m/s, 1.5 m/s and 0.5 m/s are calculated, respectively.The limit cycles are shown in Figure 11a,b.It can be seen that when v equals 3.5 m/s, the system is stable and the limit cycle amplitude is very small.When v equals 2.5 m/s, the system is unstable, and the limit cycle amplitude is greater than that of 1.5 m/s and 0.5 m/s.With the decrease of disc velocity, the limit cycle amplitude becomes larger, that is to say, the instability increases.However, the increase of the amplitude is not proportional to the decrease of v.When v continues to decrease to 0.5 m/s, the limit cycle amplitude no longer increases significantly.It is also found that the limit cycle of 2.5 m/s presents as a thick ring due to its amplitude changing (as seen in Figure 10), whereas the limit cycles amplitudes of 1.5 m/s and 0.5 m/s are constants.A large number of calculations show that when the rotational velocity of brake disc is greater than a critical value about 3.1 m/s, the system is stable.The brake system becomes unstable for a braking velocity lower than 3.1 m/s.Moreover, the lower the velocity is, the worse the stability appears.So the brake squeal tends to increase.
no longer increases significantly.It is also found that the limit cycle of 2.5 m/s presents as a thick ring due to its amplitude changing (as seen in Figure 10), whereas the limit cycles amplitudes of 1.5 m/s and 0.5 m/s are constants.A large number of calculations show that when the rotational velocity of brake disc is greater than a critical value about 3.1 m/s, the system is stable.The brake system becomes unstable for a braking velocity lower than 3.1 m/s.Moreover, the lower the velocity is, the worse the stability appears.So the brake squeal tends to increase.
Figure 11.Limit cycles of different velocity: (a) X -direction limit cycles of 3.5 m/s; (b) X -direction limit cycles of 2.5 m/s; (c) X -direction limit cycles of 1.5 m/s; (d) X -direction limit cycles of 0.5 m/s; (e) Y-direction limit cycles of 3.5 m/s; (f) Y-direction limit cycles of 2.5 m/s; (g) Y-direction limit cycles of 1.5 m/s; (h) Y-direction limit cycles of 0.5 m/s.
In conclusion, in the friction characteristic adopted in this paper, the friction coefficient is mainly affected by the value of . When k μ is a constant, the friction coefficient decreases with the increase of v .A relatively low friction coefficient can reduce the vibration amplitude of displacement and velocity of the system.If the friction coefficient is smaller than the bifurcation point, the modecoupling phenomenon disappears.Thus, the occurrence of brake squeal is suppressed.It should be noted that the friction coefficient, as a key parameter affecting braking performance, cannot be designed too small from a safety perspective.On the contrary, a sufficient small value of v can make the friction coefficient higher than the bifurcation point, and exacerbates the instability of the system.This is consistent with the fact that the brake squeal usually occurs in relatively low-speed braking conditions.

Conclusions
A large number of literatures use a minimal model to study the mechanism of brake squeal, but an accurate determination method of parameters in minimal models is rarely found.In the light of the CEA results of a finite element brake model, this paper established a more representative 2-DoF model.Taking the natural frequency, Hopf bifurcation point and real part of the eigenvalue of a finite element unstable mode as optimization objectives, the response surface method is applied to determine the kinetic parameters of the minimal model.So the minimal model can feature the modecoupling phenomenon of the unstable mode with seven nodal diameters effectively.Further analysis for brake squeal based on this model is reliable.
Based on the optimized parameters, the negative slope friction characteristic is further introduced, and TA is performed for the dynamic equations of the 2-DoF model.The results show that nonlinear variation characteristics of the friction coefficient with braking speed can affect modal coupling of the system, and therefore lead to squeal instability.The squeal is more likely to occur In conclusion, in the friction characteristic adopted in this paper, the friction coefficient is mainly affected by the value of µ k •v.When µ k is a constant, the friction coefficient decreases with the increase of v.A relatively low friction coefficient can reduce the vibration amplitude of displacement and velocity of the system.If the friction coefficient is smaller than the bifurcation point, the mode-coupling phenomenon disappears.Thus, the occurrence of brake squeal is suppressed.It should be noted that the friction coefficient, as a key parameter affecting braking performance, cannot be designed too small from a safety perspective.On the contrary, a sufficient small value of v can make the friction coefficient higher than the bifurcation point, and exacerbates the instability of the system.This is consistent with the fact that the brake squeal usually occurs in relatively low-speed braking conditions.

Conclusions
A large number of literatures use a minimal model to study the mechanism of brake squeal, but an accurate determination method of parameters in minimal models is rarely found.In the light of the CEA results of a finite element brake model, this paper established a more representative 2-DoF model.Taking the natural frequency, Hopf bifurcation point and real part of the eigenvalue of a finite element unstable mode as optimization objectives, the response surface method is applied to determine the kinetic parameters of the minimal model.So the minimal model can feature the mode-coupling phenomenon of the unstable mode with seven nodal diameters effectively.Further analysis for brake squeal based on this model is reliable.
Based on the optimized parameters, the negative slope friction characteristic is further introduced, and TA is performed for the dynamic equations of the 2-DoF model.The results show that nonlinear variation characteristics of the friction coefficient with braking speed can affect modal coupling of the system, and therefore lead to squeal instability.The squeal is more likely to occur under relatively low velocity, which is consistent with objective facts.While the relative velocity is higher than a critical value, the system is stable.Therefore, the uncertainty of braking speed is closely related to the unpredictability of the squeal.
The minimal model presented in this paper is intuitive.The model parameter determination method based on the response surface optimization fill in a gap in the literature.Moreover, compared with the FEA which consumes considerable computation time and data resources, the 2-DoF model enjoys the superiority of high efficiency in solving the dynamic response, making it easier to analyze the effect of the nonlinear variation of parameters such as friction coefficient on the limit cycle of the system.Therefore, the parameter determination method and parameter analysis results of the minimal model proposed in this paper can provide a reference for an improved design of the brake system.

Figure 1 .
Figure 1.Finite element model of brake.

Figure 2 .
Figure 2. Unstable modes of the brake.

Figure 1 .
Figure 1.Finite element model of brake.

Figure 1 .
Figure 1.Finite element model of brake.

Figure 2 .
Figure 2. Unstable modes of the brake.

Figure 2 .
Figure 2. Unstable modes of the brake.

.
coefficients greater than 0.4, the real parts of this mode are positive, which signifies an unstable mode.Moreover, the maximum real part value is 713.41 for 9 Figure4billustrates the evolution of frequencies of this mode versus the friction coefficient.With the increase of friction coefficient, two modal frequencies gradually get closer until coalescence occurs.The coalescence frequency at the bifurcation point is 7168.6Hz.In Figure4c, a three-dimensional graph including real parts and frequencies of mode coupling is given, where the red line represents the unstable mode, and the blue line represents the stable mode.

Figure 4 .
Figure 4. Complex eigenvalue analysis (CEA) results of the unstable mode: (a) the real part of eigenvalues; (b) the frequency of eigenvalues; (c) 3D plot of the coalescence patterns of the finite element mode at 7168.6 Hz.

Figure 4 .
Figure 4. Complex eigenvalue analysis (CEA) results of the unstable mode: (a) the real part of eigenvalues; (b) the frequency of eigenvalues; (c) 3D plot of the coalescence patterns of the finite element mode at 7168.6 Hz.

Figure 6 .
Figure 6.Evolution of frequency error versus stiffness 1 k .

Figure 11 .
Figure 11.Limit cycles of different velocity: (a) X-direction limit cycles of 3.5 m/s; (b) X-direction limit cycles of 2.5 m/s; (c) X-direction limit cycles of 1.5 m/s; (d) X-direction limit cycles of 0.5 m/s; (e) Y-direction limit cycles of 3.5 m/s; (f) Y-direction limit cycles of 2.5 m/s; (g) Y-direction limit cycles of 1.5 m/s; (h) Y-direction limit cycles of 0.5 m/s.

Table 1 .
Material properties of brake components.Z direction of the displacement and rotation around the Y axis, and the bolt holes on the disc are constrained with PINNED boundary condition, namely, constraining the X, Y, Z direction of the displacement.A uniform pressure of 1 MPa is applied to the back plates to simulate the actual braking hydraulic load.

Table 1 .
Material properties of brake components.

Table 1 .
Material properties of brake components.

Table 2 .
Steepest ascent computational design and measured responses.

Table 3 .
Central composite computational design and the measured response.