Prediction of Squeal Instabilities of a Finite Element Model Automotive Brake With Uncertain Structural and Environmental Parameters With a Hybrid Surrogate Model


 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.


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 [10][11][12][13]. 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.

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.

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.
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).

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: where M, C, and K are the structural mass, damping and stiffness matrices of the FEM, respectively. F p is the external forces applied on the system (i.e., the pressure on the piston) and F nl 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 F nl . The nonlinear sliding equilibrium position U S of the system is the solution of KU S + F nl (U S ) = F p . To study the stability of the solution, a small perturbation of the equilibrium position is considered and the corresponding eigenvalue problem is: where J nl is the Jacobian matrix of the non-linear forces F nl at the point U S . The eigenvalues λ j = a j + iω 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.

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.

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 L 2 sense. For numerical reasons, this series is truncated at the order P. It writes: where (Φ k ) is the polynomial basis and (α k ) is the unknown deterministic coefficient. They can be computed with different nonintrusive 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 3.1.2 Kriging. Let's consider a function α that depends on the vector x of variables x M1 and x M2 , the kriging approximates it with the following [33,34]: 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 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 = ‖x − x ′ ‖ 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.

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]: 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 Ξ = 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 M (m) a 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 M (m) a0 and is built by completing artificially the experimental design of M (m) a 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 M (m) a or M (m) ω at a new point (x 0 , ξ 0 ), the test-meta-model M (m) a0 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.

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 M 1 and M 2 can vary between 0 g and 250 g. A first experimental design composed of 50 couples (M 1 ,M 2 ) 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 Transactions of the ASME 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.

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 (M 1 ,M 2 ), 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 (M 1 ,M 2 ), 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.
Once the PCE are constructed, for each PCE coefficient, a metamodel 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 M a0 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 metamodels have been validated by comparison of the hybrid metamodel 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 k 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.

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 M 1 and M 2 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.
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 M 1 and M 2 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 M 1 = M 2 = 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 M 1 = 50 g. Then, a more or less large stability area is observed for low values of M 2 and large values of M 1 . 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 M 1 + M 2 < 250 g). For low values of M 2 and large values of M 1 , 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 M 1 = 50 g and M 2 > 50 g. So, if in average the mode has a higher real part at M 1 = M 2 = 50 g than at M 1 = M 2 = 200 g (and could be considered as "more unstable"), it is also more fleeting considering the variance as the variance is extremely low at M 1 = M 2 = 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 M 1 and M 2 , 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 M 1 null and M 2 equal to 250 g gives the lowest number of instabilities, whereas for a B (2,8) law, the minimum is reached for both M 1 null and M 2 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.
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 M 1 = 50 g and M 2 = 40 g, whereas for a uniform law, it is reached for M 1 = 75 g and M 2 = 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 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.

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.