Research Article  Open Access
Kriging Surrogate Models for Predicting the Complex Eigenvalues of Mechanical Systems Subjected to FrictionInduced Vibration
Abstract
This study focuses on the kriging based metamodeling for the prediction of parameterdependent mode coupling instabilities. The high cost of the currently used parameterdependent 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 fourdegreeoffreedom mechanical system submitted to frictioninduced 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 parameterdependent CEA that the proposed kriging surrogate model can provide efficient and reliable predictions of mode coupling instabilities with different complex patterns.
1. Introduction
Studies of mechanical systems subjected to frictioninduced 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–5]. There are typically two different analyses and categories of mechanisms available for defining the origin of frictioninduced 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 frictioninduced vibrations is rather related to kinematic constraints or spragslip 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 frictioninduced vibrations and dynamic instabilities on mechanical systems: the Complex Eigenvalue Analysis (CEA) to detect unstable frequencies [7, 8] and time analysis to determine selfexcited 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 selfexcited 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 frictioninduced 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 frictioninduced vibrations under uncertainties [15–18]. This approach has been proposed for propagating uncertainties described by probability density functions in systems submitted to frictioninduced 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–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 secondorder 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 frictioninduced 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 frictioninduced 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 fourdegreeoffreedom 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.
2. Mechanical System and Stability Analysis
2.1. Description of the Phenomenological Model
Figure 1 shows the minimal fourdegreeoffreedom model to be used in the present study. This phenomenological model has its origins in the previous twodegreeoffreedom 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–17, 34]. In the present study, an extension of this minimal Hultén model is proposed in order to investigate the case of multiinstabilities. The model consists of two masses and 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 and 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 fourdegreeoffreedom model can be expressed aswhere is the displacement vector defined by . and are the associated acceleration and velocity vectors, respectively. , , and are, respectively, the mass, damping, and structural stiffness matrices of the mechanical system. The matrix corresponds to the frictional contributions between the four plates and the three bands. Expressions of the mass matrix , the damping matrix , the structural stiffness matrix , and the stiffness matrix due to frictional forces are given by
2.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 frictioninduced 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 selfexcited 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 frictioninduced 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 frictioninduced instabilities using the well known CEA is based on the numerical analysis of the system’s eigenvalues. Assuming a solution of the form (where corresponds to the complex eigenvalues of the system and is the associated eigenvector), system (1) becomesSo 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 , 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.
3. 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 frictioninduced 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.
3.1. 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 in the dimensional space parameter and the associated eigenvalues . Each output is obtained from the solution of the parameterdependent characteristic equation:where , and are the system’s matrices evaluated at the th sample of parameter .
Based on the kriging theory [26, 27, 37], a parameterdependent eigenvalue can be expressed aswhere 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 . The second term is a realization of a Gaussian process that is assumed to have a zeromean value and a covariance function given bywhere is the process variance and 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 and . 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 is the determinant of the correlation matrix whose entries are given by with and . Then, by considering the training matrix defined from the evaluation of regression functions at the generated samples and one of the entries defined by (with and ), the regression parameter vector is defined byas the least square solution of the following regression problem: where comprises the eigenvalue samples obtained from the solutions of the parameterdependent characteristic equation (5) associated with the generated samples with . The associated variance is given by Hence, according to [26], the best linear unbiased predictor of an eigenvalue can be obtained as follows: where is the vector containing values of the correlation between each of the input sample points with and the parameter variable .
In the present study, the surrogate kriging model is constructed by using the toolbox Matlab DACE developed by Sondergaard et al. [27].
3.2. 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. and 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 realvalued and complexvalued vector and is insensitive to the phase and to the norm of and . 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.
4. 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 parameterdependent 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 and . 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 and . 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 threedimensional 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 and 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 and the stiffness coefficients and around their nominal values with dispersions of .(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.
4.1. 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 for 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]). In this study, the number of samples is fixed to if (i.e., with a discretization of 100 values for both parameters and ) and if (i.e., with a discretization of 100 values for each of the parameters , and ).
First of all, evolutions of the system’s eigenvalues (both real and imaginary parts) against the two parameters and 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 whereas the second instability appears for . The associated frequencies of the unstable modes are Hz and 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 brings forward the apparition of the coalescence point. Figure 3(e) shows the evolutions of the real and imaginary parts versus and 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.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
(e)
(f)
Then, Figures 2(b) and 2(d) illustrate the coalescence patterns for the second case. Evolutions of the real and imaginary parts versus and 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 with the frequency of the unstable mode equal to Hz; the second and third mode couplings are detected for and , respectively. The frequencies of the associated unstable modes are equal to Hz and 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 (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.
(a)
(b)
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 and the stiffness and of 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 , , and . 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.
(a)
(b)
(c)
(d)
(e)
(f)
4.2. 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.
4.2.1. Study with : Training Data from a Linear Grid
A 2D linear sample grid is generated to cover the 2D parameter space defined by the friction coefficient and the linear stiffness (i.e., ). The generated samples are used to build the training matrix and the correlation matrix 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 as well as the quadratic mean of error on the real parts denoted by are calculated as the square root of the mean of the squares as follows:where defines the number of samples. and correspond to the pulsation of the reference model and the kriging surrogate model, respectively. and correspond to the real part of one eigenvalue for the reference model and the kriging surrogate model, respectively.

