Abstract

This study focuses on the prediction of the stability behavior of an industrial automotive brake system under structural and environmental uncertainties. Uncertainties are modeled with a random distribution or an interval and are propagated with a hybrid surrogate model associating polynomial chaos and kriging. The objective is to create a surrogate model of each eigenvalue computed with the complex eigenvalue analysis (CEA). As the modes can be tracked only when unstable, the effective size of the training sets can become extremely small. Despite this limitation, it is shown the hybrid meta-model is still able to predict the stability of the brake system. Moreover, the hybrid meta-model gives a direct access to the mean and variance of the eigenvalues with respect to the design parameters without any additional Monte Carlo simulations (MCS). By considering different probability density function for the friction coefficient, it is shown it has a high influence on the stability and the latter should be accurately estimated.

1 Introduction

Friction-induced vibrations are a complex and rich topic of research since several decades for both industrials and academics [1,2]. It is well-known that the phenomenon is fugitive and depends on numerous parameters, such as the lining materials [3], the wear [4], the environment [5], the geometry of the contact surface [6], etc., making the study and the modeling complex. Thus, the consideration of the variability or the uncertainty associated to some parameters emerged a few years ago. Indeed, in 2009, Culla and Massi [7] considered some parameters as uncertain and performed Monte Carlo Samplings (MCS) on a numerical model to improve the squeal prediction. In 2010, Butlin and Woodhouse [8] investigated the variability of squeal related parameters and laid the foundations of the consideration of uncertainty for improving squeal predictions. In 2014, Tison et al. [9] demonstrated that considering uncertainty improves the squeal prediction and confirmed their results with experiments. Since, numerous studies were dedicated to the consideration and the modeling of uncertainties in brake systems from a numerical and/or experimental point of view [1013]. These studies remain complex as the different types of uncertainties must be identified and characterized, before using advanced and costly tools to take them into account.

Meta-modeling is now a classical method to propagate uncertainty for an expensive numerical model as it requires a small number of evaluations for its construction. From a few simulations, a mathematical model that surrogates the real model is constructed. Different methods exist to do so, chosen depending on the description of the uncertainty. For random parameters, the polynomial chaos expansion (PCE) [14,15], the perturbation method [16], or the FORM method [17] can be employed. The fuzzy logic may also be considered [18,19]. An other approach consists in the creation of a surface response function with regression [10,20], the kriging method [11,13,21], support vector machine, neural networks, etc.

However, different types of uncertainties might be present in the same system and the need for hybrid methods that allow different descriptions of the uncertainty in the same surrogate model emerges. Hence, Lu and Yu associate the interval approach with the response surface method (RSM) [22] for optimizing a brake system regarding squeal. In Ref. [19], they associate the RSM method with the fuzzy logic. In both cases, a RSM model of one eigenvalue of the system is created. This model depends on the interval (resp. fuzzy) variables, and then, the uncertainty is propagated through the RSM. In Ref. [23], Denimal et al. presented a new hybrid surrogate model associating kriging with PCE. It enables to deal with both random and parametric variables. The method was tested and validated on a 4-dof phenomenological model to predict its stability. If promising results were obtained, some limitations were raised. First, a high number of learning points were required to build the hybrid meta-model (about 3000 for three parameters). Thus, its application on a full finite element model (FEM) remains an open question as the number of points is limited due to large simulation time. Additionally, when dealing with real FEM of industrial brake systems, one observes a high modal density together with numerous unstable modes and complex coupling phenomena. These make the modal identification complex, and modes can only be tracked when unstable [13]. This is a major limitation as it means the functions (here the eigenvalues) to be surrogated are only partially defined and sometimes on a small portion of the design space. More concretely, on the total number of available points, only a part of them can be used to create the surrogate model of an eigenvalue, reducing drastically the number of training points. It raises the questions of the ability of the hybrid surrogate model to deal with partially defined problems with a low number of points in the experimental design. So, this work can be seen as an extension of Ref. [23] not only to demonstrate the ability and validity of the hybrid surrogate model to predict the stability of industrial FEM of brake systems but also to discuss its limitations.

The automotive brake system is subjected to three uncertain parameters, namely the friction coefficient at the pad/disc interface and two mass additions on the caliper. The first one varies from one braking event to an other and is largely related to the environmental conditions. For this reason, it is commonly described by a probability density function (PDF) [1,24,25]. Different PDFs will be considered to demonstrate its high influence and the necessity to identify it correctly. A uniform distribution but also beta distributions will be considered as they are more likely to fit a large variety of experimentally measured distributions.

The two mass additions on the caliper correspond to what one can call design parameters and hence are modeled with an interval. From a practical point of view, in a design process one will look for the best choice of mass addition considering the uncertainty on the friction coefficient. The hybrid surrogate method associating the PCE and the kriging is adopted to propagate both uncertainties simultaneously and to predict each unstable eigenvalue of the system. From the hybrid formulation, a direct access to the stochastic moments of each eigenvalue with regard to the design parameters is assessed without any additional computation. Different probability laws are considered for the friction coefficient, and in each case, the hybrid meta-modeling method is able to predict correctly the evolution of each unstable mode. The results clearly show the strong sensitivity of the brake system stability to the friction coefficient probability law.

2 Mechanical System Under Study

In this section, the mechanical system as well as the uncertain parameters are presented. Then, the numerical strategy used to predict the squeal propensity is briefly reminded.

2.1 Finite Element Model and Uncertain Parameters.

The mechanical system studied is a FEM of an industrial automotive brake system with a floating caliper technology. It is represented in Fig. 1(a). During a breaking action, a hydraulic pressure makes move the piston which pushes the inner pad against the disc. The reaction forces bring back the caliper and the outer pad, which also tackles against the disc. The contact between the pads and the disc creates friction. The contact is modeled with a linear penalty enforcement law and with a Coulomb’s law for the friction. For the sake of consistency, the implementation of the FEM and the influence of the internal contacts as well as its validation are not presented here, but the interested reader is invited to refer to Refs. [13,26] for more details.

Fig. 1
(a) Exploded view of the brake and (b) mass modification on the caliper
Fig. 1
(a) Exploded view of the brake and (b) mass modification on the caliper
Close modal

The friction coefficient μ at the pad/disc interface has a high influence on the stability of the brake. Moreover, this parameter is often not well known and varies from one braking action to another one and largely depends on the environment (temperature, humidity, surface, etc.) [1,2,8,9,27,28]. For this reason, the latter is often considered as variable to indirectly take into consideration these environmental variabilities. Moreover, results shown in Refs. [24,25] demonstrate that the friction coefficient has a statistical distribution and can be described by a probability density function. When a brake exhibits a too high squeal propensity, structural modifications are investigated to improve the brake with regard to the squeal. From an industrial point of view, a classical structural modification consists in the addition of mass at the extremities of the bracket as it is illustrated Fig. 1(b). These modifications have a high impact on the squeal propensity of the brake system [13] as they directly modify the modes that are involved in the mode coupling instabilities.

In the present study, these three parameters are considered as uncertain and will be modeled with the hybrid surrogate model. First of all, experimental data on the distribution or evolution of the friction coefficient was not accessible for the example of automotive brake system proposed in our study. This explains why the friction coefficient is taken as random. From the authors’ knowledge, the number of available experimental studies in the literature about the distribution law of the friction coefficient in industrial brakes is low. From Ref. [29], the range of variation is taken in [0.3, 0.8]. To illustrate the importance of identifying correctly the friction distribution, three laws are considered here, illustrated in Fig. 2, namely a uniform distribution, which can be seen as “neutral” and two beta laws of shape parameters (2,8) and (8,2) are considered. The choice of using beta distributions is motivated by the idea that they are more likely to fit experimental distribution than a Gaussian one. Indeed, in Refs. [1,24,25], the distributions observed for the friction coefficient are non-Gaussian. The shape parameters are chosen arbitrarily in this study for illustration. Once the random parameter is generated, a linear transformation on [0.3, 0.8] is applied to get the corresponding physical friction coefficient. On the other hand, two small masses are added on each part of the bracket and can vary between 0 g and 250 g. They correspond to design parameters, and so they are not described by a probability law. They are modeled in the FEM by a modification of the density of a few elements of the caliper, represented in Fig. 1(b).

Fig. 2
Considered probability laws for the random parameter associated with the friction coefficient μ: (a) uniform, (b) B(8,2), and (c) B(2,8)
Fig. 2
Considered probability laws for the random parameter associated with the friction coefficient μ: (a) uniform, (b) B(8,2), and (c) B(2,8)
Close modal

2.2 Complex Eigenvalue Analysis.

The squeal propensity of a brake system is determined from the complex eigenvalue analysis (CEA). It consists in the analysis of the stability of the non-linear equilibrium position.

First of all, the complete nonlinear equation of the automotive brake system can be written as follows:
(1)
where M, C, and K are the structural mass, damping and stiffness matrices of the FEM, respectively. Fp is the external forces applied on the system (i.e., the pressure on the piston) and Fnl is the non-linear forces coming from the contact and friction. It is worth emphasizing here that the mass modifications affect directly the mass matrix M, whereas the uncertainty on the friction coefficients impacts the non-linear forces Fnl.
The nonlinear sliding equilibrium position US of the system is the solution of KUS + Fnl(US) = Fp. To study the stability of the solution, a small perturbation of the equilibrium position is considered and the corresponding eigenvalue problem is:
(2)
where Jnl is the Jacobian matrix of the non-linear forces Fnl at the point US. The eigenvalues λj = aj + j are complex and ωj corresponds to the pulsation of the mode. The mode shape is given by Ψj. The equilibrium position is considered as stable if all the real parts are negative and is considered as unstable if at least one real part is positive.

3 Hybrid Surrogate Method

In this section, the hybrid method used is presented. The objective of the surrogate model is to predict the evolution of the eigenvalues with respect to the different uncertain parameters. The influence of the friction coefficient, which is described by a PDF, is modeled with the PCE, whereas the two masses are taken in consideration with a kriging surrogate model.

3.1 Mathematical Descriptions.

3.1.1 Polynomial Chaos Expansion.

Let’s consider that an eigenvalue Λ is a random variable that depends on the random friction coefficient ξμ. Then, the Generalized Polynomial Chaos (GPC) allows to approximate Λ with a polynomial expansion where the polynomials are orthogonal with respect to the PDF of ξμ [30,31]. These polynomials are given by the Askey scheme [32]. In the case of a uniform law, Legendre polynomials are employed whereas for the case of a Beta distribution, Jacobi polynomials are used. This series is convergent to L2 sense. For numerical reasons, this series is truncated at the order P. It writes:
(3)
where (Φk) is the polynomial basis and (αk) is the unknown deterministic coefficient. They can be computed with different non-intrusive techniques such as the regression method [31]. Once the coefficients (αk) are determined, the estimation of the stochastic moments is straightforward. Indeed, the mean of Λ is given by E[Λ]=α0 and the variance by σΛ2=k=0Pαk2ΦkL22.

3.1.2 Kriging.

Let’s consider a function α that depends on the vector x of variables xM1 and xM2, the kriging approximates it with the following [33,34]:
(4)
where the first part is the regression part and corresponds to the trend of the kriging meta-model. It corresponds to the projection on q functions that are often polynomials of low order. The second part, namely Z(x), is a zero-mean Gaussian process characterized by its covariance Cov[Z(x),Z(x)]=σ2k(x,x), where σ2 is the process variance and k(x, x′) is the correlation function between two input points x and x′. k(x, x′) depends on the Euclidean distance h = ‖xx′‖2 and on a vector of hyperparameters θ to be optimized. If the input dimension is superior to 1, the correlation function is constructed with a product of univariate correlation functions.

3.1.3 Hybrid Formulation.

Each eigenvalue of the stability analysis depends on the friction coefficient μ and on the caliper mass modification. The hybrid surrogate model starts from the Eq. (3) and considers that the PCE coefficients depend on the mass modification, represented by the variable x. It would be numerically impossible to create a PCE for all the mass modifications; hence, this epistemic uncertainty is modeled with a kriging meta-model. For the m-th eigenvalue, it comes [23]:
(5)

To create this meta-model, an experimental design (ED) is generated. An ED corresponds to a set of input and output, it is used to compute the parameters characterizing a meta-model (i.e., the regression coefficient for the PCE, and the regression coefficient as well as the hyper-parameter for the kriging). For the input space, two sets of points are generated. The first one is denoted Ξ={ξμ(1),,ξμ(N)} and is composed of N values of the random friction coefficient ξμ. The second one is denoted A={x(1),x(M)} and is composed of M mass couple values (M1, M2). The final ED is obtained by combining these two sets in a full factorial experiment (i.e., N × M configurations). The output set is obtained by evaluating the function to be approximated. This is done by performing a CEA for the N × M system configurations. It is worth emphasizing here that the formulation of the hybrid meta-model makes impossible to generate one ED with the three parameters together as for each x, it is necessary to construct a PCE from N evaluations.

The objective is to create two surrogate models for each eigenvalue m: one for the real part, denoted Ma(m) and one for the imaginary part, denoted Mω(m). It is thus necessary to track each eigenvalue over the variations of the system configurations. This is done by using a MAC criterion to track the associated mode shapes [11,13,23,35]. However, as demonstrated in Ref. [13], the MAC criterion is not enough to track the mode at the bifurcation point. Indeed, large variations of the mode shapes are observed at the bifurcation point, making the tracking difficult if the number of points is limited. Moreover, the two complex modes involved in a mode coupling converge to the same real mode at the bifurcation point. It is then impossible, based on a MAC criterion, to differentiate the two modes before the Hopf bifurcation.

Hence, the meta-models are created only when the real part of their eigenvalue is non-zero and are valid only on this space. This has a substantial impact on the ED as on the total number of available points, only a few will be available for the surrogate model creation. Moreover, each mode has its own experimental design that corresponds to the part of the total experimental design where the mode is (un)stable (i.e., with a strictly positive or negative real part). Because it does not correspond to a pavement of the input space, a third meta-model is built and acts as a test function to check if the real part is equal to zero or not. The latter is denoted Ma0(m) and is built by completing artificially the experimental design of Ma(m) with a zero output at the points where the mode was not tracked [13]. This third meta-model takes zero values when the mode has an eigenvalue with a zero real part, and a non-zero value otherwise (equal to the real part of the eigenvalue), and depending on its value, one can know if the considered mode is stable or not, and so if the surrogate models are in their validity area or not. So, before evaluating a meta-model Ma(m) or Mω(m) at a new point (x0, ξ0), the test-meta-model Ma0(m) is evaluated to ensure that the two meta-models are in their definition space. If not, one can consider the mode is stable and has a zero-real part. A small threshold equal to 5 is taken in practice to avoid any oscillatory effect.

Finally, for each mode m and for each mass couple (i.e., for each x(k)A), a PCE is constructed based on the N model evaluations for which the mode is unstable. The PCE coefficients are the (αj(k),(m))j{0,P}. Once all the PCE are constructed, for each PCE coefficient, a kriging meta-model of the PCE coefficient (αj(m)) is build based on the available data (i.e., the M values where the mode is unstable). The general workflow is illustrated in Fig. 3.

Fig. 3
General workflow of the PCE–Kriging construction
Fig. 3
General workflow of the PCE–Kriging construction
Close modal

4 Application for Squeal Prediction

The section here is devoted to the use of the hybrid surrogate model for the prediction of the squeal propensity of the brake system. In a first time, the construction of the experimental designs is explained, then the meta-models are created and finally the hybrid meta-model is used for the prediction of squeal.

4.1 Construction of the Experimental Designs.

As explained previously, three uncertain parameters are considered in the present study: the friction coefficient and two masses added to the bracket. Since the friction coefficient is unknown, it is described by a random law. Three laws are considered here, namely the uniform law and two beta laws of shape parameters (2,8) and (8,2) represented in Fig. 2. The two masses M1 and M2 can vary between 0 g and 250 g. A first experimental design composed of 50 couples (M1,M2) is generated with a latin hypercube sampling (LHS) for the kriging, and a second experimental design composed of 15 points is generated by LHS for the friction coefficient for the PCE. Finally, 50 × 15 = 750 input points are used for the hybrid-meta-model creation and the corresponding CEA are performed. For each considered friction law, an experimental design is generated.

Figure 4 shows the stability analysis for the three different laws. Each mode is tracked with an MAC criterion and identified by one color in the complex plans. A large number of unstable modes appear for each configuration with potentially complex crossover phenomena between the different unstable modes. These results show without any ambiguity the complexity of the phenomena involved and consequently the potential difficulty to create an accurate hybrid surrogate model for predicting the stability of a FEM automotive brake system. The influence of the law on the squeal propensity is obvious, and the unstable modes as well as their evolution ranges are also impacted. For example, the mode at 4050 Hz is unstable only when the friction coefficient follows a B(2,8) law. Or, the mode at 4300 Hz has a real part that can reach about 2000 when the friction coefficient follows a uniform law or a B(8,2), but reaches only 1100 for a B(2,8) law. In the following, only the modes with a large enough experimental design are considered for the creation of the hybrid meta-model (i.e., more than 50 points available). This is motivated by two things: first, to have a surrogate model of good quality one needs enough training points; second, if the mode is unstable only on a small part of the experimental design, then one may neglect its unstable behavior in a design process. The threshold of 50 points is arbitrary here, to ensure enough points in the training set. The different modes retained are summarized in Table 1, where for each mode one gives its frequency, its nature, the maximum of the real part (minimum if stable) and the number of points in the experimental design for each law. Again, one can see that the uncertainty related to the friction coefficient has a large impact on the system stability: some modes are unstable only for one PDF, the stability area depends on the PDF (directly related to the number of points in the ED) and variation of the maximum (or minimum) value of the real part of the eigenvalue.

Fig. 4
Experimental design from the mode tracking for the different laws: (a) uniform law, (b) B(8,2) law, and (c) B(2,8) law
Fig. 4
Experimental design from the mode tracking for the different laws: (a) uniform law, (b) B(8,2) law, and (c) B(2,8) law
Close modal
Table 1

Modes involved in mode coupling for the different experimental designs with their frequency, real part, and the number of points in each experimental design for each law

Mode no.Frequency (Hz)NatureUniform lawB(8,2) lawB(2,8) law
Real partNo. pts EDReal partNo. pts EDReal partNo. pts ED
1879.7Unstable27.511222.9100
21855.1Unstable248.3403248.3475217.0324
32853.2Unstable131.1168131.122089.580
43460.8Unstable124.1601387.869696.5459
53863.5Unstable670.2741774.2750406.4678
64395.7Unstable1948.97391948.97501194.9674
74879.5Unstable249.940982.6275
85541.1Unstable54.9144
9879.7Stable−27.5112−22.9100
101855.1Stable−248.3403−248.3475−217.0324
112853.2Stable−131.1168−131.1220−89.580
123460.8Stable−124.1601−387.8696−96.5459
133863.5Stable−670.2741−774.2750−406.4678
144395.7Stable−1948.9739−1948.9750−1194.9674
154879.5Stable−249.9409−82.6275
165541.1Stable−54.9144
Mode no.Frequency (Hz)NatureUniform lawB(8,2) lawB(2,8) law
Real partNo. pts EDReal partNo. pts EDReal partNo. pts ED
1879.7Unstable27.511222.9100
21855.1Unstable248.3403248.3475217.0324
32853.2Unstable131.1168131.122089.580
43460.8Unstable124.1601387.869696.5459
53863.5Unstable670.2741774.2750406.4678
64395.7Unstable1948.97391948.97501194.9674
74879.5Unstable249.940982.6275
85541.1Unstable54.9144
9879.7Stable−27.5112−22.9100
101855.1Stable−248.3403−248.3475−217.0324
112853.2Stable−131.1168−131.1220−89.580
123460.8Stable−124.1601−387.8696−96.5459
133863.5Stable−670.2741−774.2750−406.4678
144395.7Stable−1948.9739−1948.9750−1194.9674
154879.5Stable−249.9409−82.6275
165541.1Stable−54.9144

4.2 Construction and Validation of the Meta-Models.

Once the experimental designs are generated, the hybrid meta-models are constructed in two steps. First, the PCE is used to estimate the evolution of the frequencies and the real parts of the eigenvalues with respect to the friction coefficient μ. Thus, for each couple (M1,M2), a PCE is constructed for each retained mode. For the uniform law and the beta laws, Legendre and Jacobi polynomials are used, respectively. The degree of each chaos is chosen to be the same for all the modes in one case (same law and same surrogate model) and to obtain the best quality of prediction. To do so, the average relative error between the PCE predictions and the exact values at the points of the ED is minimized. These degrees are summarized in Table 2. For each couple (M1,M2), the number of available points in the experimental design related to the PCE might be small, especially if the mode is mostly stable on the considered interval which explains the low degrees retained for the PCE. As an illustration, the PCE prediction for the real part of the fifth eigenvalue is given in Fig. 5(a). Red points correspond to the training set and the black line to the PCE prediction. One can see the good accuracy of the prediction. Similar results are obtained for the other modes.

Fig. 5
Real part of the eigenvalue of the mode 5: (a) comparison of the PCE prediction (black line) to reference points (red dot) for (M1, M2) = (212.1, 137.3) g and (b) comparison of the PCE-Kriging prediction (black cross) to different validation points (red circle) for different values of μ, M1, and M2 (Color version online.)
Fig. 5
Real part of the eigenvalue of the mode 5: (a) comparison of the PCE prediction (black line) to reference points (red dot) for (M1, M2) = (212.1, 137.3) g and (b) comparison of the PCE-Kriging prediction (black cross) to different validation points (red circle) for different values of μ, M1, and M2 (Color version online.)
Close modal
Table 2

Polynomial chaos order for the different laws and the different surrogate models

LawMωMaMa0
Uniform457
B(8,2)553
B(2,8)333
LawMωMaMa0
Uniform457
B(8,2)553
B(2,8)333
Once the PCE are constructed, for each PCE coefficient, a meta-model with the ordinary kriging method is built with a Matérn 5/2 correlation kernel. The kriging predictions of the first two coefficients of the PCE of the Ma0 surrogate model of the mode 3 for the three laws are given in Fig. 6. The red points are the learning points, and the surface is the kriging prediction. The influence of the chosen law is quite obvious, the surfaces are complex and vary from one law to another. In the case of the coefficient 0, which is directly related to the mean, the surfaces have similar shapes but the numerical values are different. Indeed, the maximum is equal to 90, 110, and 82 for the uniform law, the B(8,2) law and the B(2,8) law, respectively. For the coefficient 1, both the shapes and the values are strongly influenced by the law. From this, its comes that the mean of the real part of the eigenvalue of the mode at 2850 Hz has a similar evolution (but more or less higher values) for different PDF of the friction coefficient, but that the variance will be strongly impacted. The different meta-models have been validated by comparison of the hybrid meta-model prediction to reference values. To do so, a set of 45 reference points has been generated based on additional simulations. For each mode, the PCE-kriging predictions are compared to these reference values. The number of validation points depends not only on the considered mode but also on the probability law as the surrogate models are not valid on the full design space. As an illustration, the comparison between the predictions for the real part of the eigenvalue of the fifth mode and the reference values are given in Fig. 5(b) where red circles are the reference values and black crosses the predictions. The accuracy of the surrogate model is obvious. This comparison is summarized for each real and imaginary parts of each eigenvalue and for each law in Table 3, where the average relative error e¯ for each mode is given by
(6)
where yk(pred) is the prediction of the considered meta-model at the kth validation point and yk(ref) the reference value. To be noted that the quantities, yk(pred) and yk(ref) are associated with the real or imaginary parts of the eigenvalue for one specific mode. n corresponds to the number of points in the experimental design for which the real part of the associated mode is superior to zero. Considering more specifically the average relative error e¯ for imaginary parts, all the results indicate unambiguously the relevance and validity of the PCE-kriging predictions. Considering results for real parts, the level of error e¯ is always low giving confidence in the different hybrid surrogate models except for some limited specific cases. These specific cases for which an important error e¯ is detected correspond to modes for which we have both a very small number of points n to calculate e¯ as well as positive real parts very close to zero which has the consequence of increasing the relative error e¯.
Fig. 6
Kriging prediction of the first two PCE coefficients for the real part of the mode 3 for the three laws with respect to the masses M1 and M2—ED (): (a,d) uniform law, (b,e) B(8,2) law, and (c,f) B(2,8) law (Color version online.)
Fig. 6
Kriging prediction of the first two PCE coefficients for the real part of the mode 3 for the three laws with respect to the masses M1 and M2—ED (): (a,d) uniform law, (b,e) B(8,2) law, and (c,f) B(2,8) law (Color version online.)
Close modal
Table 3

Average relative error e¯ of the PCE-kriging predictions at validation points on the eigenvalues estimation for the different laws

Error real partError imaginary part
UniformB(8,2)B(2,8)UniformB(8,2)B(2,8)
Mode 18.51E-21.042.57E-45.07E-4
Mode 24.45E-21.95E-12.96E-11.01E-39.53E-32.34E-3
Mode 37.35E-20.42E-22.97E-16.38E-42.52E-42.1E-3
Mode 42.63E-11.31E-14.82E-11.12E-38.12E-47.00E-4
Mode 51.42E-22.59E-23.75E-21.74E-47.82E-44.17E-4
Mode 67.65E-26.28E-31.11E-11.32E-31.14E-31.33E-3
Mode 74.88E-17.65E-24.10E-48.28E-5
Mode 84.48E-12.87E-4
Mode 98.51E-21.042.57E-45.07E-4
Mode 104.51E-22.55E-12.96E-11.04E-39.53E-32.39E-3
Mode 117.35E-22.97E-16.38E-42.14E-3
Mode 123.61E-12.73E-14.69E-15.36E-21.66E-37.26E-4
Mode 131.42E-32.59E-23.44E-21.74E-47.82E-44.34E-4
Mode 147.65E-26.28E-31.1E-11.32E-31.15E-31.33E-3
Mode 154.88E-17.65E-24.10E-48.28E-5
Mode 164.47E-12.87E-4
Error real partError imaginary part
UniformB(8,2)B(2,8)UniformB(8,2)B(2,8)
Mode 18.51E-21.042.57E-45.07E-4
Mode 24.45E-21.95E-12.96E-11.01E-39.53E-32.34E-3
Mode 37.35E-20.42E-22.97E-16.38E-42.52E-42.1E-3
Mode 42.63E-11.31E-14.82E-11.12E-38.12E-47.00E-4
Mode 51.42E-22.59E-23.75E-21.74E-47.82E-44.17E-4
Mode 67.65E-26.28E-31.11E-11.32E-31.14E-31.33E-3
Mode 74.88E-17.65E-24.10E-48.28E-5
Mode 84.48E-12.87E-4
Mode 98.51E-21.042.57E-45.07E-4
Mode 104.51E-22.55E-12.96E-11.04E-39.53E-32.39E-3
Mode 117.35E-22.97E-16.38E-42.14E-3
Mode 123.61E-12.73E-14.69E-15.36E-21.66E-37.26E-4
Mode 131.42E-32.59E-23.44E-21.74E-47.82E-44.34E-4
Mode 147.65E-26.28E-31.1E-11.32E-31.15E-31.33E-3
Mode 154.88E-17.65E-24.10E-48.28E-5
Mode 164.47E-12.87E-4

4.3 Squeal Prediction.

Once the meta-models are constructed, it is possible to predict the evolution of the different eigenvalues. The main interest of the present methodology is the estimation of the PCE coefficients with the kriging. Since the PCE coefficients are directly related to the stochastic moments, it is possible from the kriging to predict straightforwardly the evolution of the mean and the variance of the eigenvalues without any MCS as it would have been the case if a unique kriging meta-model was built for the three parameters [13]. Thus, means and variances with respect to the masses M1 and M2 can be directly evaluated. For example, the evolution of the mean and the variance of the real part of the modes 3 and 4 are given in Figs. 7 and 8 for the three different laws.

Fig. 7
Evolution of the mean of the real part with regard to the masses M1 and M2 for the modes 3 (top) and 4 (bottom) for the ((a,d)) uniform law, (b,e) B(8,2) law, and (c,f) B(2,8) law
Fig. 7
Evolution of the mean of the real part with regard to the masses M1 and M2 for the modes 3 (top) and 4 (bottom) for the ((a,d)) uniform law, (b,e) B(8,2) law, and (c,f) B(2,8) law
Close modal
Fig. 8
Evolution of the variance of the real part with regard to the masses M1 and M2 for the modes 3 (top) and 4 (bottom) for the (a,d) uniform law, (b,e) B(8,2) law, and (c,f) B(2,8) law
Fig. 8
Evolution of the variance of the real part with regard to the masses M1 and M2 for the modes 3 (top) and 4 (bottom) for the (a,d) uniform law, (b,e) B(8,2) law, and (c,f) B(2,8) law
Close modal

Considering the stability behavior of the mode 3, looking at the average of the real part of the eigenvalue (see in Fig. 7 upper line), the behavior and the evolution is similar over the three different PDF of the friction coefficient. Indeed, the mode is mostly stable, and when M1 and M2 increase simultaneously, the mode becomes unstable. This unstable area is slightly smaller for a B(2,8) law. If the evolution is similar, the maximum mean value of the real part is impacted by the friction PDF and larger mean values are expected with a B(8,2) law for example. On the opposite, the behavior of the variance is strongly impacted by the PDF (see in Fig. 8 upper line). For a uniform law a large variance is observed on the whole unstable part, for a B(8,2) law the variance is low and the maximum is observed at the bifurcation area around M1 = M2 = 50 g and for a B(2,8) law a large variance is observed only at the most unstable area. The extremely large variance for the uniform law is of importance, as it might imply a change in the stability behavior of the mode depending on the friction coefficient value, and so depending on the environmental conditions a squeal instability at this frequency can be observed or not. But if the friction coefficient is described by a B(2,8) or a B(8,2) law, the variance on the unstable part is low and the mode will always be unstable whatever the friction coefficient value and the environmental conditions.

Looking at the stability behavior of the mode 4, the evolution of the mean of the real part is similar but it impacts the stability behavior of the mode (see in Fig. 7 bottom line). First, the average real part is always much higher when M1 = 50 g. Then, a more or less large stability area is observed for low values of M2 and large values of M1. This stability area is almost non-existent for a uniform friction law, and much larger when the friction coefficient follows a B(2, 8) law. When the two masses have higher values, the mode is also unstable but with lower values of the real part in average for all cases. In this case, the PDF that describes the friction coefficient has an impact on the average stability behavior of the mode. As for mode 3, the variance is strongly impacted by the PDF of the friction law (see in Fig. 8 bottom line). For the uniform law, a large variance is observed in the bottom left part (area where M1 + M2 < 250 g). For low values of M2 and large values of M1, it means that the mode is sometimes stable and so that squeal event might take place, or not, depending on the environment. For the B(8,2) case, the variance is always very low and the stability behavior of the mode is not impacted. For the B(2,8) case, the variance is always low except at M1 = 50 g and M2 > 50 g. So, if in average the mode has a higher real part at M1 = M2 = 50 g than at M1 = M2 = 200 g (and could be considered as “more unstable”), it is also more fleeting considering the variance as the variance is extremely low at M1 = M2 = 200 g.

Finally, it is possible to access directly the mean number of instability and the sum of the mean of the positive real parts with regards to the masses M1 and M2, they are given in Fig. 9 for each law. The influence of the law for the friction coefficient is obvious. Indeed, when a uniform law is considered, from 2.5 to 5 instabilities are observed in mean, whereas for a B(8,2) law, up to six instabilities can be observed in average. Moreover, the location of the maximum and minimum number of instabilities is different depending on the law. Indeed, for a uniform law, a mass M1 null and M2 equal to 250 g gives the lowest number of instabilities, whereas for a B(2,8) law, the minimum is reached for both M1 null and M2 inferior to 80 g. Also, the shape of the limits between the average number of instabilities (see black lines in Fig. 9) is completely different from one case to another.

Fig. 9
Mean number of unstable modes (top) and sum of the mean of the positive real parts (bottom) with regard to the masses M1 and M2 for (a,d) the uniform law, (b,e) the law B(8,2), and (c,f) the law B(2,8). Black lines in (a,b,c) denote the contours of integer numbers of instabilities.
Fig. 9
Mean number of unstable modes (top) and sum of the mean of the positive real parts (bottom) with regard to the masses M1 and M2 for (a,d) the uniform law, (b,e) the law B(8,2), and (c,f) the law B(2,8). Black lines in (a,b,c) denote the contours of integer numbers of instabilities.
Close modal

A similar analysis can be done for the sum of the mean of the positive real parts, given in Figs. 9(d), 9(e), and 9(f). Indeed, the mean values are highly impacted by the distribution: see the different ranges of variations in Figs. 9(b), 9(d), and 9(f). For a B(2,8) law, the highest value is equal to 700, whereas for a B(8,2), it is equal to 2600. The locations of the maximum and minimum are also impacted. Indeed, for a B(2,8) law, the minimum is reached for M1 = 50 g and M2 = 40 g, whereas for a uniform law, it is reached for M1 = 75 g and M2 = 240 g.

Results obtained here demonstrate the complex influence of the design parameters on the system stability which illustrates the necessity to conduct deep studies to end in a correct understanding of squeal phenomenon. The evolution of the different modes is not trivial which complicates the identification of the best set of design parameters for an engineer. Moreover, the identification of the probability laws of the random parameters must be done beforehand and precisely as it has a high impact on the system stability.

4.4 Discussion.

The results presented give a concrete use of a hybrid surrogate model for the prediction of stability of an industrial brake system, but also raises some limitations:

  • – size of the ED (here 750): it is quite high and is definitely a limitation related to the use of this method, but this high number is related to different factors. The first reason is the full factorial construction of the ED imposed by the formulation of the hybrid surrogate model as the uncertainties are propagated one after another, and work should be conducted to lift this limitation. However, this high number of points is mostly due to the physical problem under consideration. Indeed, due to the high modal density and the important variation of the mode shapes over the design space, the ED must be dense enough to ensure a robust mode tracking, which implies a minimal number of points. Moreover, as modes can be tracked only when the eigenvalues real parts are non-zero, the ED must be dense enough to ensure that each instability is indeed detected, but it also implies that on the total number of points, only a limited number is finally available to create a surrogate model for a mode. So, when 750 simulations are run, in some cases only 80 points are really available (see mode 3) at the end for the construction of the hybrid meta-model. So finally, the performances of the hybrid meta-model should be considered with the effective number of available points. And, finding a criterion more robust than a MAC criterion for the mode tracking and able to pass the bifurcation would be the solution to lift this difficulty.

  • – comparison with a kriging with all the parameters: one could naturally consider using one kriging meta-model for the three parameters [13] and then perform MCS for the propagation of the random law related to the friction coefficient. Using this strategy represents a reduction of the initial number of CEA to perform. However, the MCS to be performed at the end takes time, several minutes as many modes are present and as the validity of the surrogate models must be tested each time. By using the hybrid surrogate model, this simulation time is removed as the stochastic moments are directly obtained from the PCE coefficients. In some cases, it might be more convenient to have a more important simulation time initially and get directly statistical moments at the end, than having to run again MCS after the surrogate model creation. One last advantage, not illustrated here as only one parameter is random, is that the Sobol indices can also be directly accessed from the PCE coefficients.

  • – computational time: except the simulations related to the CEA, the most computationally expensive part of the method is the mode tracking part, with the computation of many MAC matrices over the ED and represents a few minutes in the present case. The construction of the hybrid meta-models represents only 1.5 s on average, showing the efficiency of the method. Once the surrogate models are constructed, obtaining the evolution of the mean and variance on a 50 × 50 grid (Fig. 7) takes only 0.13 s for a mode, when a few minutes would be required with a MCS on a single kriging surrogate model.

5 Conclusion

The present study proposes the application of a hybrid surrogate method for the prediction of the stability behavior of an industrial brake system subjected to many instabilities. The main objective is to demonstrate the validity and the interest of hybrid methods when a system is put through different types of uncertainties. Indeed, for many application epistemic and aleatory uncertainties are present, and modeling them at the same time in the same surrogate model is not necessarily the best representation.

The CEA methodology is employed for the prediction of the stability behavior of the finite element model of the brake system and three uncertain parameters are considered. First, the friction coefficient is described by three different probability laws, and two design parameters corresponding to two small masses that might be added to the caliper are described by an interval of possibility. The friction coefficient is modeled via a PCE whereas the two design parameters are modeled through the kriging method. Once the meta-models are constructed and validated, the analysis of the stability of the brake system can be performed. The work demonstrates the applicability of such an approach for large models where computational times are important and where the construction of the learning sets is limited due to the mode tracking, which limits drastically the number of available points to construct the meta-models.

From the PCE coefficients, it is possible to access without any additional computations (i.e., no MCS) to the evolution of the mean and the variance of each eigenvalue regarding the design parameters allowing for a stochastic description of the system. This would not be possible if a unique surrogate model for all the uncertain parameters was created, since MCS simulations would be performed to propagate the uncertainty related to the random parameters. The approach used here may require more simulations to create the surrogate models, but no cost is added then for the UQ part. Considering the global process, this strategy might be of interest in some cases when stochastic properties are required.

Moreover, the strong impact of the probability law chosen for the friction coefficient on the system stability is also observed and demonstrates the necessity to identify it in advance.

Acknowledgment

This work was achieved within Stellantis Stellab program-OpenLab Vibro-Acoustic-Tribology@Lyon (VAT@Lyon). J.-J. Sinou acknowledges the support of the Institut Universitaire de France.

Conflict of Interest

There are no conflicts of interest.

References

1.
Ibrahim
,
R. A.
,
1994
, “
Friction-Induced Vibration, Chatter, Squeal, and Chaos Part I: Mechanics of Contact and Friction
,”
ASME Appl. Mech. Rev.
,
47
(
7
), pp.
209
226
.
2.
Ibrahim
,
R. A.
,
1994
, “
Friction-Induced Vibration, Chatter, Squeal, and Chaos Part II: Dynamics and Modeling
,”
ASME Appl. Mech. Rev.
,
47
(
7
), pp.
227
253
.
3.
Eriksson
,
M.
,
Bergman
,
F.
, and
Jacobson
,
S.
,
1999
, “
Surface Characterisation of Brake Pads After Running Under Silent and Squealing Conditions
,”
Wear
,
232
(
2
), pp.
163
167
.
4.
Bergman
,
F.
,
Eriksson
,
M.
, and
Jacobson
,
S.
,
1999
, “
Influence of Disc Topography on Generation of Brake Squeal
,”
Wear
,
225
, pp.
621
628
.
5.
Eriksson
,
M.
,
Lundqvist
,
A.
, and
Jacobson
,
S.
,
2001
, “
A Study of the Influence of Humidity on the Friction and Squeal Generation of Automotive Brake Pads
,”
Proc. Inst. Mech. Eng., Part D: J. Auto. Eng.
,
215
(
3
), pp.
329
342
.
6.
Bergman
,
F.
,
Eriksson
,
M.
, and
Jacobson
,
S.
,
2000
, “
The Effect of Reduced Contact Area on the Occurrence of Disc Brake Squeals for An Automotive Brake Pad
,”
Proc. Inst. Mech. Eng., Part D: J. Auto. Eng.
,
214
(
5
), pp.
561
568
.
7.
Culla
,
A.
, and
Massi
,
F.
,
2009
, “
Uncertainty Model for Contact Instability Prediction
,”
J. Acoust. Soc. Am.
,
126
(
3
), pp.
1111
1119
.
8.
Butlin
,
T.
, and
Woodhouse
,
J.
,
2010
, “
Friction-Induced Vibration: Quantifying Sensitivity and Uncertainty
,”
J. Sound. Vib.
,
329
(
5
), pp.
509
526
.
9.
Tison
,
T.
,
Heussaff
,
A.
,
Massa
,
F.
,
Turpin
,
I.
, and
Nunes
,
R. F.
,
2014
, “
Improvement in the Predictivity of Squeal Simulations: Uncertainty and Robustness
,”
J. Sound. Vib.
,
333
(
15
), pp.
3394
3412
.
10.
,
H.
, and
Yu
,
D.
,
2014
, “
Brake Squeal Reduction of Vehicle Disc Brake System With Interval Parameters by Uncertain Optimization
,”
J. Sound. Vib.
,
333
(
26
), pp.
7313
7325
.
11.
Nobari
,
A.
,
Ouyang
,
H.
, and
Bannister
,
P.
,
2015
, “
Uncertainty Quantification of Squeal Instability Via Surrogate Modelling
,”
Mech. Syst. Signal. Process.
,
60
, pp.
887
908
.
12.
Renault
,
A.
,
Massa
,
F.
,
Lallemand
,
B.
, and
Tison
,
T.
,
2016
, “
Experimental Investigations for Uncertainty Quantification in Brake Squeal Analysis
,”
J. Sound. Vib.
,
367
, pp.
37
55
.
13.
Denimal
,
E.
,
Sinou
,
J-J.
, and
Nacivet
,
S.
,
2019
, “
Influence of Structural Modifications of Automotive Brake Systems for Squeal Events With Kriging Meta-Modelling Method
,”
J. Sound. Vib.
,
463
, p.
114938
.
14.
Sarrouy
,
E.
,
Dessombz
,
O.
, and
Sinou
,
J.-J.
,
2013
, “
Piecewise Polynomial Chaos Expansion With An Application to Brake Squeal of a Linear Brake System
,”
J. Sound. Vib.
,
332
(
3
), pp.
577
594
.
15.
Nechak
,
L.
,
Besset
,
S.
, and
Sinou
,
J.-J.
,
2018
, “
Robustness of Stochastic Expansions for the Stability of Uncertain Nonlinear Dynamical Systems–Application to Brake Squeal
,”
Mech. Syst. Signal. Process.
,
111
, pp.
194
209
.
16.
Nobari
,
A.
,
Ouyang
,
H.
, and
Bannister
,
P.
,
2015
, “
Statistics of Complex Eigenvalues in Friction-Induced Vibration
,”
J. Sound. Vib.
,
338
, pp.
169
183
.
17.
,
H.
, and
Yu
,
D.
,
2015
, “
Stability Analysis and Improvement of Uncertain Disk Brake Systems with Random and Interval Parameters for Squeal Reduction
,”
ASME J. Vib. Acoust.
,
137
(
5
), p.
051003
.
18.
Do
,
H. Q.
,
Massa
,
F.
,
Tison
,
T.
, and
Lallemand
,
B.
,
2017
, “
A Global Strategy for the Stability Analysis of Friction Induced Vibration Problem With Parameter Variations
,”
Mech. Syst. Signal. Process.
,
84
, pp.
346
364
.
19.
,
H.
,
Shangguan
,
W-B.
, and
Yu
,
D.
,
2017
, “
Uncertainty Quantification of Squeal Instability Under Two Fuzzy-Interval Cases
,”
Fuzzy Sets Syst.
,
328
, pp.
70
82
.
20.
Treimer
,
M.
,
Allert
,
B.
,
Dylla
,
K.
, and
Müller
,
G.
,
2017
, “
Uncertainty Quantification Applied to the Mode Coupling Phenomenon
,”
J. Sound. Vib.
,
388
, pp.
171
187
.
21.
Nechak
,
L.
,
Gillot
,
F.
,
Besset
,
S.
, and
Sinou
,
J-J.
,
2015
, “
Sensitivity Analysis and Kriging Based Models for Robust Stability Analysis of Brake Systems
,”
Mech. Res. Commun.
,
69
, pp.
136
145
.
22.
,
H.
, and
Yu
,
D.
,
2016
, “
Optimization Design of a Disc Brake System With Hybrid Uncertainties
,”
Adv. Eng. Soft.
,
98
, pp.
112
122
.
23.
Denimal
,
E.
,
Nechak
,
L.
,
Sinou
,
J.-J.
, and
Nacivet
,
S.
,
2018
, “
A Novel Hybrid Surrogate Model and Its Application on a Mechanical System Subjected to Friction-Induced Vibration
,”
J. Sound. Vib.
,
434
, pp.
456
474
.
24.
Soom
,
A.
, and
Kim
,
C.
,
1983
, “
Interactions Between Dynamic Normal and Frictional Forces During Unlubricated Sliding
,”
ASME J. Lubr. Tech.
,
105
(
2
), pp.
221
229
.
25.
Oberst
,
S.
, and
Lai
,
J. C. S.
,
2011
, “
Statistical Analysis of Brake Squeal Noise
,”
J. Sound. Vib.
,
330
(
12
), pp.
2978
2994
.
26.
Denimal
,
E.
,
Sinou
,
J.-J.
,
Nacivet
,
S.
, and
Nechak
,
L.
,
2019
, “
Squeal Analysis Based on the Effect and Determination of the Most Influential Contacts Between the Different Components of an Automotive Brake System
,”
Int. J. Mech. Sci.
,
151
, pp.
192
213
.
27.
Kinkaid
,
N. M.
,
O’Reilly
,
O. M
, and
Papadopoulos
,
P.
,
2003
, “
Automotive Disc Brake Squeal
,”
J. Sound. Vib.
,
267
(
1
), pp.
105
166
.
28.
Ouyang
,
H.
,
Nack
,
W.
,
Yuan
,
Y.
, and
Chen
,
F.
,
2005
, “
Numerical Analysis of Automotive Disc Brake Squeal: A Review
,”
Int. J. Vehicle Noise Vib.
,
1
(
3–4
), pp.
207
231
.
29.
Bergman
,
F.
,
Eriksson
,
M.
, and
Jacobson
,
S.
,
1999
, “
Influence of Disc Topography on Generation of Brake Squeal
,”
Wear
,
225
, pp.
621
628
.
30.
Wiener
,
N.
,
1938
, “
The Homogeneous Chaos
,”
Am. J. Math.
,
60
(
4
), pp.
897
936
.
31.
Sudret
,
B.
,
2008
, “
Global Sensitivity Analysis Using Polynomial Chaos Expansions
,”
Reliab. Eng. Syst. Saf.
,
93
(
7
), pp.
964
979
.
32.
Xiu
,
D.
, and
Karniadakis
,
G. E.
,
2003
, “
Modeling Uncertainty in Flow Simulations Via Generalized Polynomial Chaos
,”
J. Comput. Phys.
,
187
(
1
), pp.
137
167
.
33.
Forrester
,
A.
,
Sobester
,
A.
, and
Keane
,
A.
,
2008
,
Engineering Design Via Surrogate Modelling: a Practical Guide
,
John Wiley & Sons
,
Chichester, UK
.
34.
Roustant
,
O.
,
Ginsbourger
,
D.
, and
Deville
,
Y.
,
2012
,
Dicekriging, Diceoptim: Two R Packages for the Analysis of Computer Experiments by Kriging-Based Metamodeling and Optimization
.
35.
Denimal
,
E.
,
Nechak
,
L.
,
Sinou
,
J-J.
, and
Nacivet
,
S.
,
2016
, “
Kriging Surrogate Models for Predicting the Complex Eigenvalues of Mechanical Systems Subjected to Friction-Induced Vibration
,”
Shock Vib.
,
2016
, p.
3586230
.