Kriging Surrogate Models for Predicting the Complex Eigenvalues of Mechanical Systems Subjected to Friction-Induced Vibration

This study focuses on the kriging based metamodeling for the prediction of parameter-dependent mode coupling instabilities. The high cost of the currently used parameter-dependent Complex Eigenvalue Analysis (CEA) has induced a growing need for alternative methods. Hence, this study investigates capabilities of kriging metamodels to be a suitable alternative. For this aim, kriging metamodels are proposed to predict the stability behavior of a four-degree-of-freedom mechanical system submitted to friction-induced vibrations. This system is considered under two configurations defining two stability behaviors with coalescence patterns of different complexities. Efficiency of kriging is then assessed on both configurations. In this framework, the proposed kriging surrogate approach includes a mode tracking method based on the Modal Assurance Criterion (MAC) in order to follow the physical modes of the mechanical system. Based on the numerical simulations, it is demonstrated by a comparison with the reference parameter-dependent CEA that the proposed kriging surrogate model can provide efficient and reliable predictions of mode coupling instabilities with different complex patterns.


Introduction
Studies of mechanical systems subjected to friction-induced vibrations benefit from a growing interest due to the large amount of applications in the field of mechanical engineering.The different and complex mechanisms that can be responsible for undesirable dynamic characteristics and appearance of instabilities in many mechanical systems have been extensively studied in the last decades [1][2][3][4][5].There are typically two different analyses and categories of mechanisms available for defining the origin of friction-induced system instability: the first one is mainly due to tribological properties whereas the second one relies on geometrical conditions.While the variation of the friction coefficient is considered as one of the most important factors for the emergence of instability in the first category (i.e., in the case of a tribological approach), the origin of friction-induced vibrations is rather related to kinematic constraints or sprag-slip phenomenon [6] and modal coupling in the second case (i.e., in the case of a structural dynamics approach based on geometrical conditions).In this last case, the emergence of instability can be detected even with a constant friction coefficient.In the present study, this last approach that is based on structural coupling mechanism will be discussed.
Nowadays, two kinds of analysis are classically used to undertake numerical studies of friction-induced vibrations and dynamic instabilities on mechanical systems: the Complex Eigenvalue Analysis (CEA) to detect unstable frequencies [7,8] and time analysis to determine self-excited vibrations [9,10].As explained in previous papers [9,11,12], both approaches have their pros and cons.However CEA based methods and the calculations of self-excited vibrations may become too costly when parametric analysis and/or uncertainty propagation are needed for engineering design problems [13].In these cases, it may be worthwhile to work towards the development of sophisticated methods based on surrogate models in order to perform design optimization or design space approximation (i.e., emulation).The main aim is to substitute any complex model by a suitable surrogate model which offers a convenient compromise between the accuracy of its predictions and the cost related to its implementation.In the present study, one is interested in estimating the occurrence of instability in a predefined design space approximation.In this context, the main purpose of the surrogate modeling is the generation of a surrogate that is as accurate as possible for the prediction of the occurrence of instabilities in the complete design space of interest, using as few simulation evaluations as possible.Such approximation models, known as metamodels or emulators, mimic the behavior of the simulation model (i.e., estimation of all the real and imaginary parts of eigenvalues in our case) as closely as possible while being computationally cheaper to evaluate.It may be noted that the accuracy of the surrogate depends on the number and location of samples in the design space of interest required for its implementation.Moreover, surrogate models are characterized by some tuning parameters that control their accuracy.
In the field of friction-induced vibrations, numerous formalisms have been developed to define surrogate models for the prediction of mode coupling instabilities.Surrogate models that are based on the Generalized Polynomial Chaos (GPC) formalism [14] have been proposed this last decade to deal with the stability of mechanical systems subjected to friction-induced vibrations under uncertainties [15][16][17][18].This approach has been proposed for propagating uncertainties described by probability density functions in systems submitted to friction-induced instabilities, a task which is prohibitive when performed by using the Monte Carlo method.The latter was exploited for estimating of the probability of squeal occurrence in [19] and in several other studies as a reference method [15][16][17].So the main idea governing the GPC formalism consists of expressing the system's degrees of freedom or eigenvalues within a functional space built from polynomials that are orthogonal with respect to probabilistic measures associated with the system's design parameters.The chaos order is the most important tuning parameter which is fixed to a suitable value from a convergence study.This probabilistic surrogate model has shown an interesting efficiency in propagating and quantifying uncertainties on the stability behavior of such systems.However, it may present some limits when the number of uncertain parameters is relatively high and/or when high chaos orders are required in particular for functions that are strongly nonlinear in the random space.In this last case, surrogate models based on the multielement GPC can be useful [20].The response surface methodology (RMS) is also proposed in [21] to deal with the stability, reliability, and sensitivity analysis of brake systems submitted to interval and random uncertainties.The proposed surrogate model consists of using basis functions (defined by monomials in the uncertain parameters) to express the system eigenvalues.
The same surrogate model is proposed in [22] for the optimization design of brake system under interval and random uncertainties.Another approach consists of constructing surrogate models based on the perturbation principle [23].In this specific case, the main principle consists in expressing system's eigenvalues by means of Taylor expansions near the mean value of uncertain parameters.For example, the firstorder perturbation method has been proposed by Butlin and Woodhouse [24] to quantify sensitivity of frictioninduced instabilities to the design parameters.Moreover the recent study of Nobari et al. [25] proposes a second-order Taylor expansion to estimate statistical properties of eigenvalues characterizing mode coupling instabilities.Despite its efficiency, this approach has limitations, especially when standard deviations of parameters are important.
Another type of surrogate models can be constructed based on the kriging method [26,27].This approach exploits spatial correlations between a small number of function values at some samples generated from a suitable experience plan to predict unknown values of the function within its design space of interest.The kriging based model consists in two main parts; the regression model roughly represents the global tendency of the analyzed function while the second part is defined by a stochastic process representing the spatial correlations in the design space of interest.Extensive reviews of kriging metamodeling in simulation and other applications as in the sensitivity analysis and optimization in design process can be found in [28,29].
In the recent past, the prediction of squeal instability in brake systems via surrogate modeling has been introduced by Nobari et al. [30] and Nechak et al. [13] in order to construct a predictor of squeal instability.However, these two interesting studies were limited to the estimation of unstable frequencies without considering the prediction of all the stable modes or the behavior of the system before the emergence of squeal instability.Extension of the prediction of both stable and unstable behaviors (i.e., estimation of all the real parts and imaginary parts of eigenvalues) for mechanical systems subjected to friction-induced vibrations via surrogate modeling is proposed in the present study.Furthermore, the construction of a surrogate model for each separated mode is carried out with a careful selection of output data, due to the evolution of the complex eigenvalues and the order of the modes when some specific parameters change.If the surrogate model is not constructed by paying attention to this point, errors due to an improper surrogate model will appear.To overcome this difficulty, a tracking process based on the Modal Assurance Criterion (MAC) is proposed in conjunction with the generation of kriging surrogate models characterized by their tuning parameters, namely, the order of the regression model, the spatial correlation model, and the size of the experience plan.The two previous studies [13,30] have not analyzed exhaustively the effects of these tuning parameters on the accuracy of kriging for instability predictions of mechanical system subjected to friction-induced vibrations.
So the main objective and originality of the present paper lies in the analysis of performances (in terms of accuracy and cost) of kriging surrogate models with respect to their tuning parameters, when dealing with the prediction of not only mode coupling instabilities but also the complete approximation of the real and imaginary parts of both stable and unstable modes.To undertake such a study and to validate the methodology of constructing kriging surrogate models with a tracking process, numerical simulations will be performed on a minimal four-degree-of-freedom model.Two specific numerical cases will be investigated: the first one will be a classical baseline with "a simple mode coupling mechanism" (i.e., coalescence of two modes, one being unstable and the other unstable).The second case deals with more complex modes coupling mechanisms with the successive appearances and disappearances of instabilities, and the crossing phenomenon between modes.This study is organized as follows.First, the mechanical system under study is presented and the classical stability analysis (i.e., CEA) is briefly discussed.Then, the methodology for constructing kriging surrogate models is developed.The proposed procedure allows not only the prediction of mode coupling instabilities but also the prediction of all the complex eigenvalues of the mechanical system (i.e., the real parts and imaginary parts for both stable and unstable modes) by performing a selection and arrangement of the modes via a suitable MAC criterion.The last part of the present study is devoted to the presentation and discussion of the numerical results.

Mechanical System and Stability Analysis
2.1.Description of the Phenomenological Model. Figure 1 shows the minimal four-degree-of-freedom model to be used in the present study.This phenomenological model has its origins in the previous two-degree-of-freedom model proposed by Hultén [31,32] and was used to point out the role of the damping and the destabilization paradox [33] and for the prediction of mode coupling instabilities submitted to parameter uncertainties [15][16][17]34].In the present study, an extension of this minimal Hultén model is proposed in order to investigate the case of multi-instabilities.The model consists of two masses  1 and  2 held against moving bands disposed as shown in Figure 1.Contacts between masses and bands are modeled by plates supported by springs and damping.The masses  1 and  2 are linearly coupled by a spring (i.e., stiffness   ) and the associated damping   .
For the sake of simplicity, it is assumed that the two masses and the three band surfaces are always in contact due to a preload applied to the mechanical system.Considering the friction forces between the four plates and the three bands, the coefficient of friction  is assumed to be constant.Its value is the same for all contacts and the classical Coulomb law is applied.The velocity of the moving bands is considered as constant.Moreover, the relative velocities between the band speed and the displacements of masses are assumed to be positive so that the direction of the tangential friction force   does not change.According to the Coulomb law, the tangential friction force is assumed to be proportional to the normal force   (i.e.,   =   , where  is the friction coefficient).
Equation of motion for the four-degree-of-freedom model can be expressed as where X is the displacement vector defined by X = ( 1  2  3  4 )  .Ẍ and Ẋ are the associated acceleration and velocity vectors, respectively.M, C, and K are, respectively, the mass, damping, and structural stiffness matrices of the mechanical system.The matrix K  corresponds to the frictional contributions between the four plates and the three bands.Expressions of the mass matrix M, the damping matrix C, the structural stiffness matrix K, and the stiffness matrix K  due to frictional forces are given by ] ] . (2)

Prediction of Instabilities Based on the Complex Eigenvalues Analysis.
From the literature, it is admitted that there are typically two different and complementary methodologies in order to predict instabilities of mechanical systems subjected to friction-induced vibrations: the CEA and the dynamic transient and stationary analysis.Both methodologies have their advantages and disadvantages and they can be performed separately.As previously explained by Sinou et al. [9,10], the calculation of the nonlinear transient or/and stationary self-excited vibrations and the estimation of the acoustic noise [35,36] define the most relevant process to detect the presence or absence of instabilities in mechanical systems subjected to friction-induced vibrations.However, the approach is often unfortunately too expensive computationally.From this point of view, methods based on the CEA are often used in order to predict instabilities in a given frequency range.Even if it may lead to an underestimation or an overestimation of the unstable modes observed in the nonlinear time simulation, it is nowadays of common use in industry.One of the main limitations of CEA lies in the fact that the predictions of the onset of instability are valid just locally (i.e., in the neighborhood of a given static equilibrium point).
The prediction of friction-induced instabilities using the well known CEA is based on the numerical analysis of the system's eigenvalues.Assuming a solution of the form X() = Φ  (where  corresponds to the complex eigenvalues of the system and Φ is the associated eigenvector), system (1) becomes So stability of the system can be investigated by performing an eigenvalue analysis of the characteristic equation Due to the contribution of the nonsymmetric stiffness matrix K  , the mechanical system may become unstable.The stability analysis based on the CEA is then based on the system's eigenvalues given by the complex roots   =   +   of the characteristic polynomial.Indeed, according to the Lyapunov theory, the asymptotic stability of system (1) is stated if all eigenvalues are with strict negative real parts while the system is said to be unstable (i.e., appearance of one or more instabilities) if at least one eigenvalue has a positive real part   .The corresponding imaginary part   defines the pulsation of the unstable mode.The number of unstable modes (i.e., instabilities) is related to the number of eigenvalues with a positive real part.

Kriging Models for Stability Analysis
In this section, the use of surrogate models in the context of the prediction of complex eigenvalues for mechanical systems subjected to friction-induced vibration will be presented.
The construction of a surrogate model will enable us to reproduce all the outputs of expensive simulation (i.e., the prediction of all the eigenvalues based on the CEA in the present case) while requiring a limited number of simulations and thereby avoiding considerable resources in terms of both computation time and data storage.
In the next parts of this section, the mathematical formulation of kriging for the prediction of eigenvalues will be first presented.Secondly, the need to implement a tracking process based on a MAC criterion for all the complex eigenvalues will be discussed.

Mathematical Formulation of Kriging.
The kriging, also named Gaussian process, is an interpolation method that permits estimating a response surface of a function just from a relatively small number of simulations performed at sample points generated randomly or pseudorandomly from the parameter space of the function.The exploited principle is based on the correction of a rough approximation given by a linear or nonlinear regression model by using a zeromean Gaussian process characterized by a spatial correlation function which estimates the similarity of two points in the parameter space.The essential ideas of kriging are presented in the sequel through the considering of estimating of the response surface of eigenvalues.In this perspective, let us consider  sample points p = (p () ) =1,..., in the dimensional space parameter and the  associated eigenvalues  () .Each output  () is obtained from the solution of the parameter-dependent characteristic equation: where M () , C () , K () , and K ()  are the system's matrices evaluated at the th sample of parameter p.
Based on the kriging theory [26,27,37], a parameterdependent eigenvalue can be expressed as where the first term is  basis polynomial functions   that are weighted by the regression parameters   .It describes the global tendency of the approximated function against parameters p.The second term  is a realization of a Gaussian process that is assumed to have a zero-mean value and a covariance function given by where  2 is the process variance and R(, s, p) ∈ [0, 1] is the spatial correlation function with the scaling parameter , [⋅] being the expectation operator.The correlation function is a monotone function depending on the distance between p and s.This function is build such that two identical points have a unitary correlation when two infinitely separated points have zero correlation.It has the following form: Basic functions are given in Table 1.
Three parameters have to be determined to define the kriging model (6), namely, the regression parameter , the process variance , and the scaling parameter .Since  and  are -dependent, the latter has to be first estimated.This is performed by solving the following maximum likelihood function: where |R| is the determinant of the correlation matrix R whose entries are given by   = R(, p () , p () ) with  = 1, . . .,  and  = 1, . . ., .Then, by considering the  ×  training matrix F defined from the evaluation of regression functions at the generated  samples and one of the entries defined by   =   (p () ) (with  = 1, . . .,  and  = 1, . . ., ), the regression parameter vector  is defined by as the least square solution of the following regression problem: where Λ comprises the  eigenvalue samples obtained from the  solutions of the parameter-dependent characteristic equation ( 5) associated with the generated samples p () with  = 1, . . ., .The associated variance  2 is given by Hence, according to [26], the best linear unbiased predictor of an eigenvalue can be obtained as follows: where r(p) = [R(, p (1) , p) ⋅ ⋅ ⋅ R(, p () , p)]  is the vector containing values of the correlation between each of the  input sample points p () with  = 1, . . .,  and the parameter variable p.
In the present study, the surrogate kriging model is constructed by using the toolbox Matlab DACE developed by Sondergaard et al. [27].

Efficient Surrogate Modelling for Squeal Instability Based
on Mode Shape Criterion.Because the evolution of the complex eigenvalues and the order of the associated modes can change when the parameters change, it is necessary to follow the evolution of all separated eigenvalues and the associated eigenvectors computed via the CEA in order to construct a reliable surrogate model on complex eigenvalues.
One of the most efficient ways to allow an accurate tracking of the mode evolution is to compare two sets of modes and perform the pairing between these sets of modal vectors using a MAC.The MAC that corresponds to a measure of the degree of linearity or consistency between one modal and another reference modal vector is given by where * designate the conjugate transpose of a complex vector.Φ 1 and Φ 2 define the mode shape (i.e., eigenvector) of the baseline system and the mode shape of the modified system.This criterion can be used for both real-valued and complex-valued vector and is insensitive to the phase and to the norm of Φ 1 and Φ 2 .It is also appropriate to the case of monophased vectors.Due to the fact that the MAC criterion is based upon the minimization of the squared error between two vector spaces, one of the main drawbacks on the use of the MAC criterion is its sensibility to components of large value in the mode shapes (i.e., eigenvectors).Moreover, its low sensibility when decreasing the number of components in the eigenvectors is another potential limitation to performing pairing between two sets.This classical MAC has been previously used successfully by Nobari et al. [30] in order to be sure that the same unstable modes are used for constructing a predictor.The authors proposed to correlate the modified unstable mode shapes of the randomized inputs (due to randomizing the material properties of the reference system) with the baseline deterministic design.

Numerical Results
This section discusses performances of kriging metamodels with respect to their tuning parameters (i.e., the order of the regression function, the spatial correlation function, and the size of the experience plan) in the prediction of the parameter-dependent stability and instability regions in the design space.In this perspective, two configurations are considered with the set of parameters presented in Table 2.
The two studied configurations differ from each other with respect to the complexity of the phenomenon observed in the coupling modes that is mainly due to the coupling degree defined by the stiffness parameter   between the two masses  1 and  2 .The first configuration may be considered as a classical baseline with "a simple mode coupling mechanism": it corresponds to two decoupled coalescences of two modes (for each coalescence one mode is unstable and the other is stable).The second configuration gives rise to a more complex coalescence pattern including more complex modes coupling mechanisms with the successive appearances and disappearances of instabilities and the crossing phenomenon between stable and unstable modes.These behaviors are mainly due to the "high coupling" between the two masses  1 and  2 .The main aim in the consideration of the two configurations is to observe the effect of the growing complexity of the proposed system on the kriging performances.
Moreover, two-and three-dimensional parameter spaces are considered in order to analyze the influence of the dimension of the parameter space on the kriging performance.So two studies with p = {, 11 } and p = {, 11 ,  22 } are performed, respectively.
For the reader comprehension, the following proposed studies and the associated results are conducted as follows: (i) The "reference" is computed without any surrogate model.A classical parametric study in connection with a stability analysis based on the CEA method is performed with both configurations and both sets.The obtained results constitute the database used to assess the validity and the performances of the constructed kriging metamodels.They are obtained by considering the friction coefficient  ∈ [0,1] and the stiffness coefficients  11 and  22 around their nominal values with dispersions of 10%.
(ii) Calculations are performed using kriging metamodels for the predicting of the complex eigenvalues.In this case, the selected samples for the construction of kriging models are taken homogeneously in the design space.This choice is not necessarily optimal but is made to study more specifically different regression models in term of order (zero, one, and two) and different spatial correlation models (linear, Gaussian, exponential, cubic, and spherical).
(iii) Calculations based on kriging surrogate models using Latin Hypercube Sampling (LHS) are conducted.The main aim of this study is to analyze efficiency of kriging metamodels combined with a reduced sample data set.

Reference Model.
As mentioned previously, stability analysis based on the CEA method is performed by sampling the -dimensional space parameter (i.e., generating the set of samples p () for  = 1, . . .,  and by solving, for each sample, the characteristic equation ( 5)).The main drawback of this method is related to its high cost induced by the high number  of solutions required to cover the whole parameter space and thus to ensure confident stability analysis.It must also be noted that this cost is much more important when dealing with large scale systems (see [13]).First of all, evolutions of the system's eigenvalues (both real and imaginary parts) against the two parameters  and  11 are plotted in Figure 2 for the two configurations.Two different behaviors of the mechanical system are clearly observed.As shown in Figures 2(a) and 2(c), the first one gives rise to two independent classical mode coupling phenomena (i.e., two decoupled single mass Hultén systems).Two independent pairs of eigenvalues can be observed.For each pair, the corresponding frequencies approximate each other when the friction coefficient increases until coalescence.Near this point the real parts of the eigenvalues separate.First instability is detected for  > 0.3 whereas the second instability appears for  > 0.58.The associated frequencies of the unstable modes are 65.3 Hz and 41.7 Hz, respectively.Each coalescence pattern appears to be not affected by the other mode coupling.It is also observed that increasing the value of the stiffness parameter  11 brings forward the apparition of the coalescence point.Figure 3(e) shows the evolutions of the real and imaginary parts versus  and  11 in the complex plan.For the reader comprehension, it is worth remembering that all the numerical results obtained via the CEA do not use the tracking process based on the MAC criterion for these first case studies.However, we chose to use the MAC criterion in order to draw four different color maps for both the real and imaginary parts.In this framework, the MAC criterion is used only for the sake of graphical representation of results.
Then, Figures 2(b) and 2(d) illustrate the coalescence patterns for the second case.Evolutions of the real and imaginary parts versus  and  11 in the complex plan are also shown in Figure 3(f).Although a "classical simple mode coupling phenomenon" was observed for the first configuration, this second case study presents more complex coalescence patterns.Three mode couplings are detected: the first one appears for  = 0.13 with the frequency of the unstable mode equal to 59.3 Hz; the second and third mode couplings are detected for  > 0.48 and  > 0.81, respectively.The frequencies of the associated unstable modes are equal to 66.6 Hz and 45.1 Hz, respectively.Moreover, it is observed that the modes are involved in several successive coupling coalescences (see, e.g., the two modes represented by the orange and yellow maps).At last, a crossing phenomenon between stable and unstable modes is also observed around  = 0.55 (see the two modes represented by the orange and yellow maps in Figure 2(d)).
In order to investigate the efficiency of kriging surrogate models and to offer easy viewing of metamodels for both real and imaginary parts in comparison with the reference based on CEA, we propose to draw envelopes (i.e., the minimum and maximum values) for each real and imaginary part of eigenvalues while keeping the representation of variations according to the control parameter .For example, Figures 3(a)-3(d) show the representation of each envelope for both real and imaginary parts for the two previous cases.Moreover, evolutions of all eigenvalues in the complex plane are used to compare results obtained by kriging surrogate models and the reference.These results for the reference model are plotted in Figures 3(e) and 3(f) while the number of unstable modes characterizing the stable and unstable design spaces is also displayed versus the two design parameters in Figure 4 for the two considered configurations.Moreover, a numerical study by tracking the system's eigenvalues against three design parameters is performed in order to investigate the potential of constructing kriging surrogate models: the friction coefficient  in [0, 1] and the stiffness  11 and  22 of 10% around their mean values are considered.Results displayed in Figures 5(e) and 5(f) represent the evolution of eigenvalues in the complex plan assuming the variation of ,  11 , and  22 .Evolutions of the envelopes (i.e., the minimum and maximum values) versus the control parameter  are shown in Figures 5(a)-5(d).All these results that are performed via the CEA will be used as the reference data in order to assess the accuracy of the constructed kriging surrogate models in the following studies.

Kriging Surrogate Model.
Performances of kriging based metamodels in the prediction of mode coupling instabilities are assessed with respect to their tuning parameters, namely, the regression order, the spatial correlation model, and the size of the training data.Both configurations with the two parameter sets used previously are also considered.Moreover, two experience plans are used to generate the training data.The first one is a linear grid.The objective of this first part is to conduct a study rather dedicated to the impact of the tuning parameters by considering homogeneous input data on the design space.The second experience plan is based on the LHS that is a widely used method to generate controlled random samples.The final goal here is to analyze the capacity of kriging metamodels to deal with different complexity levels induced by the dimension of the input parameters space and various coalescence patterns.

Study with p = {𝜇, 𝑘 11 }: Training Data from a Linear
Grid.A 2D linear sample grid (20 × 50) is generated to cover the 2D parameter space defined by the friction coefficient  and the linear stiffness  11 (i.e., p = {,  11 }).The generated samples are used to build the training matrix F and the correlation matrix R associated with different regression models in term of order (zero, one, and two) and different spatial correlation models (linear, Gaussian, exponential, cubic, and spherical).We are interested by the evolutions of envelops of the system's eigenvalues (real and imaginary parts) versus the friction coefficient .
The quadratic mean of error on the imaginary and real parts of each eigenvalues and the associated maximum of absolute error for different regression and correlation kriging functions are given in Tables 3, 4, 5, and 6 for the first configuration (in Tables 7, 8, 9, and 10 for the second configuration).The quadratic mean of error on the imaginary parts denoted by   QM as well as the quadratic mean of error on the real parts denoted by  real QM are calculated as the square root of the mean of the squares as follows: where  defines the number of samples. ref and  kriging correspond to the pulsation of the reference model and the kriging surrogate model, respectively. ref and  kriging correspond to the real part of one eigenvalue for the reference model and the kriging surrogate model, respectively.Results show influences of both regression and correlation functions on the accuracy of kriging based predictions.Moreover, they are different with respect to the considered configuration.Each eigenvalue requires a kriging metamodel with specific regression and correlation function (not necessarily the same for all eigenvalues) to ensure suitable levels of accuracy.Indeed, for the first configuration as example, it can be observed that the accuracy of kriging predictions  measured by the quadratic mean of absolute errors is relatively suitable for all the built metamodels.However, for the same eigenvalue (real and imaginary parts), the average accuracy is weakly sensitive to the tuning parameters (the regression order and the correlation function type).A relative variation of the accuracy of predictions is rather observed between kriging metamodels of different eigenvalues.The maximum of absolute error indicates that the kriging metamodel constructed with the general exponential function gives predictions with the best and less sensitive accuracies on all eigenvalues (real and imaginary part).Indeed, the sum of maximum errors recorded for all eigenvalues is stationary (2.39% for the real part and 2.33% for the imaginary part).Kriging based predictions (by using the general exponential function with a second-order regression function) and reference predictions of envelopes of eigenvalues (real and imaginary parts) induced by the variability of parameter  11 are plotted versus the friction coefficient  in Figures 6(a) and 6(c).Moreover, Figure 6(e) shows evolutions of eigenvalues in the complex plan for the first configuration.It is clearly shown that kriging based predictions globally show the same tendencies as the reference so the whole stability behavior characterized by the coalescence patterns within the considered parameter intervals is faithfully represented.Indeed, the kriging predictor, built by using the general exponential correlation function and the second-order regression function, has suitably modeled the influence of  11 parameter on the first instability in which the envelope in the -dimension is suitably estimated at [0.29, 0.39].Otherwise, the kriging based metamodel has also well predicted the zeroimpact of  11 on the second instability located at  = 0.57 and on the associated coupling modes.Secondly, evolutions of all eigenvalues (real and imaginary parts) are plotted for the second configuration by using the general exponential function with a second-order regression function.Comparisons of the reference model and the chosen kriging surrogate model are illustrated in Figure 6.More specifically, Figures 6(b) and 6(d) give envelopes of imaginary and real parts of eigenvalues versus p = {,  11 }. Figure 6(f) gives evolutions of eigenvalues in the complex plan.It is clearly illustrated that the complex coalescence patterns are very well reproduced with the three mode couplings: the modes are involved in several successive coupling coalescences with a crossing phenomenon between stable and unstable modes.Envelopes of the corresponding instabilities in the -dimension are obtained within the intervals [0.12, 0.27] and [0.81, 0.82], respectively.Moreover, the complex plane for the reference and the kriging surrogate model are very similar.So increasing the system's complexity (as proposed for the second configuration) by increasing the coupling stiffness   has not strongly impacted the average accuracy of kriging based predictions (see Tables 7 and  9).Finally, it is also noted that the general exponential correlation function has also allowed a less sensitive accuracy towards eigenvalues (real and imaginary parts) and the regression order.This is indicated by the stationary sum of maximum absolute errors on kriging predictions (2.01% for the real part and 2.33% for the imaginary part), as given in Tables 8 and 10.Otherwise, it is noted from Figure 7 that the number of unstable modes versus design parameters  and  11 is well predicted by the kriging models constructed for the two considered configurations.Indeed, the comparison with the reference deterministic results shows the suitable accuracy of the stability and instability zones estimated by the constructed kriging metamodels.From all these results, the efficiency of the proposed kriging surrogate model is clearly demonstrated.The minimum and maximum envelopes are suitably approximated for all the real parts and frequencies.So the kriging model is able to reproduce not only the classical baseline of the first configuration (i.e., coalescence of two modes, one being unstable and the other stable) but also the more complex modes coupling mechanism with the successive appearances of instabilities and the crossing phenomenon between modes (as proposed for the second configuration).

Study with p = {𝜇,𝑘 11 }:
Training Data from a Latin Hypercube Sampling.In the previous part, the training process of the generated kriging metamodels is performed by using data generated from linear grilling of the parameter space.In fact, linear grid is not necessarily the optimal plan for an efficient training giving a suitable compromise between    the accuracy and the cost (size of the training data).So, in this part, the pseudorandom sampling based on the LHS is used.
It is a widely used method to generate controlled random samples.It is based on the basic idea to make sampling point distribution close to probability density function.Opposite to linear grids, the performance of the training process depends not only on the size of the LHS plan but also on the location of the generated samples in the parameter space as samples are generated pseudorandomly.Thus, in order to analyze performances of kriging with respect to the size of LHS plan by taking into account the location of samples in the parameter space, we have considered the convergence study of kriging from the statistical view point.This is performed by generating 50 times the LHS plan for different size.The mean value and the standard deviation of relative errors associated with the predicted envelopes of eigenvalues (real and imaginary parts) are determined versus the size of the used LHS plans.Results are displayed in Figure 8.  Results show a clear decrease of the mean relative error versus the size of LHS plan.An increasing number of points induces a diminution of the mean relative error as well as the associated standard deviation.For example for the first eigenvalue, stationary regimes for the mean relative error and the standard deviation are reached from  = 300 points.At this level, the accuracy of kriging prediction is weakly dependent on the repartition of training data within the parameter space.The reached stationary error level is 8.7% for the real part while it is lower for the imaginary part 0.23%.Otherwise, it can be remarked that the size  required for a convergent training with an LHS plan is strongly smaller than the one required with the linear grid ( = 1000).Then Figure 9 shows, for the two considered configurations, kriging based The stability behavior is analyzed versus the two parameters  and  11 .This is carried out by determining the number of unstable modes.The latter is estimated by the kriging model, displayed and compared to the reference results in Figure 10.It can be observed that the stability and instability regions are accurately predicted.Otherwise, it can be noted that even if more complex coalescence patterns are considered, it is still possible to reproduce the stability behavior of the system and the successive modes couplings by using the kriging model with LHS.
In conclusion, the efficiency of the proposed kriging surrogate model with LHS is clearly demonstrated.So, for almost the same level of error, the size  of the training data with an LHS plan is drastically reduced in comparison with the previous linear grid (see Section 4.2.1).

4.2.3.
Study with p = {,  11 ,  22 }.In this last section, the kriging metamodel is considered to predict mode coupling instabilities within a 3-dimensional parameter space defined as p = {,  11 ,  22 }.Tuning parameters determined previously with the 2-dimensional study are kept in the following.The effect of the increased complexity in terms of the dimension of the space parameter can then be analyzed.So the general exponential correlation and second-order regression functions are used.The training data are generated from a 3D LHS plan for parameters  11 and  22 with dispersions equal to 10% around their mean values and  within the unitary interval.Its size has been fixed to  = 200.Envelopes of the system's eigenvalues induced by the dispersions of considered parameter are estimated and plotted in the -dimension.Results are presented in Figure 11 for both configurations.Figures 11(c) and 11(d) have to be compared to Figures 5(c) and 5(d), respectively.
For the first configuration, comparing Figures 11(a) and 11(c) to Figures 5(a) and 5(c), respectively, the variation of the stiffness  22 has an effect on eigenvalues with the two lower pulsations but not on the others.This sensitivity towards the  22 parameter is well modeled by the used kriging metamodel as well as the stability behavior.The coalescence patterns are suitably predicted so increasing the dimension of the parameter space has not negatively impacted the performance of the kriging metamodel.
For the second configuration, comparing Figures 11(b) and 11(d) to Figures 5(b) and 5(d), respectively, the envelopes of eigenvalues are more spread and an eigenvalue can have a positive real part from  = 0, and the system is always unstable.In this case, the kriging metamodel has also well predicted the spread character of the envelopes and the resulting unstable behavior.Increasing the dimension of the parameter space has involved the increase of the size of the LHS and thus the size of the training data to capture the induced  22 effect on the eigenvalues.
In order to probe the contribution of the use of kriging metamodels comparing to the classical parametric CEA, computational costs induced by both methods in the stability analysis are given in Tables 11 and 12.These results highlight the less computational cost of the proposed method with LHS and the nonnegligible contribution of using kriging metamodels.Moreover the numbers of CEA realized in order to perform the stability analysis with kriging metamodels are drastically reduced as shown in Table 13.The reduction level is more important when the kriging surrogate model is determined using the LHS plan.

Conclusion
The potential of kriging based metamodeling in the prediction of mode coupling instabilities has been analyzed in this study.This analysis has been carried out in order to propose a suitable alternative to the parameter-dependent CEA which is known to be too costly in most cases.In this perspective, a four-degree-of-freedom mechanical system subjected to friction-induced vibrations and mode coupling instabilities has been considered with two configurations.These yield two different stability behaviors characterized by mode coupling having coalescence patterns with different complexities.The aim was to assess performances of kriging metamodels with respect to the complexity of the stability behavior.Results have shown an interesting efficiency for kriging metamodels comparing to the parameter-dependent CEA method since, for the same level of accuracy, the number of calculuses has been shown to be drastically decreased.This efficiency is however shown to be dependent on the complexity of the stability behavior.Indeed, among all the kriging tuning parameters, the size of the experience plan is the most impacted one.The size is closely dependent on the number of coalescence points and thus the number of mode coupling instabilities needed to be predicted.Hence, results in this study give strong indicators for kriging aptitudes to be efficiently used for the prediction of squeal which is a work in progress.Otherwise, the design space of the used system has been defined to capture some complexity patterns related to mode coupling instabilities, the main aim being the assessing of kriging abilities to model and predict these instabilities.
Other phenomena with different complexity levels are not necessary captured with the considered parameter set.This issue will be addressed in a future study.

Figure 2 :
Figure 2: Evolution of real parts (a, b) and imaginary parts (c, d) of eigenvalues versus  and  11 for configurations 1 (a, c) and 2 (b, d).

Figure 3 :Figure 4 :
Figure 3: Two-dimensional parameter space: real (a, b) and imaginary (c, d) parts of eigenvalues versus  and evolution in the complex plan (e, f) for configurations 1 (a, c, e) and 2 (b, d, f).

Figure 5 :
Figure 5: Three-dimensional parameter space: real (a, b) and imaginary (c, d) parts of eigenvalues versus  and evolution in the complex plan (e, f) for configurations 1 (a, c, e) and 2 (b, d, f).

Figure 8 :
Figure 8: Evolution of the quadratic mean error and the standard deviation error versus the size of the experimental design on real (a, b) and imaginary (c, d) parts of each eigenvalue for configurations 1 (a, c) and 2 (b, d).

Table 2 :
Set of parameters under study for the mechanical model.

Table 3 :
Configuration 1: quadratic mean of absolute error on real parts of eigenvalues for different regression and correlation kriging functions.

Table 4 :
Configuration 1: maximum of absolute error on real parts of eigenvalues for different regression and correlation kriging functions.

Table 5 :
Configuration 1: quadratic mean of absolute error on imaginary parts of eigenvalues for different regression and correlation kriging functions.

Table 6 :
Configuration 1: maximum of absolute error on imaginary parts of eigenvalues for different regression and correlation kriging functions.

Table 7 :
Configuration 2: quadratic mean of the absolute error on real parts of eigenvalues for different regression and correlation kriging functions.

Table 8 :
Configuration 2: maximum of the absolute error on real parts of eigenvalues for different regression and correlation kriging functions.

Table 9 :
Configuration 2: quadratic mean of absolute error on imaginary parts of eigenvalues for different regression and correlation kriging functions.

Table 10 :
Configuration 2: maximum absolute error on imaginary parts of eigenvalues for different regression and correlation kriging functions.

Table 11 :
Comparison of computation time between the CEA and the metamodel for two parameters.

Table 12 :
Comparison of computation time between the CEA and the metamodel for three parameters.

Table 13 :
Comparison of the number of samples for each of the case studies p = {,  11 } and p = {,  11 ,  22 } (i.e., for 2 or 3 parameters) in order to compute the deterministic results via CEA or to construct a kriging metamodel.p = {,  11 } p = {,  11 ,  22 }