Generalized Modal Amplitude Stability Analysis for the prediction of the nonlinear dynamic response of mechanical systems subjected to friction-induced vibrations

The numerical prediction of the dynamic behaviour of mechanical systems subjected to friction-induced vibrations is still a tedious problem. Different methodologies exist nowadays to study it. The first one is the complex eigenvalue analysis, which is widely used by the scientists and the industrials to predict the appearance of instabilities despite its disadvantages. Other methodologies, namely temporal integration and frequential approaches, have been developed to determine the transient and/or the steady-state response to assess the history of the dynamic response, and so to identify the unstable modes involved in the nonlinear dynamic response as well as the vibration levels. However, because of their complex implementation, their high numerical cost and sometimes the strong assumptions made on the form of the solutions, these methods are not widely and currently used in industry. To cope with the limitations of the CEA, namely the over- or under-predictability and the lack of information about modal participations in the nonlinear dynamic response, developing complementary tools is necessary. Thus, this paper is devoted to the extension and generalization of a nonlinear approach, called the modal amplitude stability analysis, to the multi-instability case. The method, called the Generalized Modal Amplitude Stability Analysis (GMASA), allows to identify the evolutions and contributions of unstable modes involved in the nonlinear self-sustaining vibration response and to estimate the limit cycles. The method is applied on a phenomenological system for which it is easy to provide an understanding of the unstable mode(s) contribution to the nonlinear dynamic response of the system and for which the calculations can be performed with reasonable computational times. Thus, the efficiency and validity of the GMASA approach are investigated by comparing the GMASA results with those of the reference results based on temporal approach.


Introduction
Instabilities for mechanical systems subjected to frictioninduced vibrations (squeal, for example) are a complex phenomenon studied by both industrials and academics since several decades [1][2][3]. At the beginning, the studies focused on the tribological aspect of the problem and a velocity friction dependency was considered as the squeal mechanism [4]. Then, Spurr [5] developed a sprag-slip model that highlights squeal with a constant friction coefficient. A generalization of this phenomenon was proposed by some researchers [6][7][8] by considering geometrical coupling as the origin of friction-induced vibrations and squeal noise. This mode coupling theory allows the appearance of the squeal phenomenon by assuming a constant friction coefficient. This approach is often accepted as one of the first causes of friction-induced vibration and squeal noise in the automotive industry.
In terms of numerical strategies, different methodologies exist to predict the appearance of instability. The first strategy is the well-known Complex Eigenvalue Analysis (CEA), which consists in the determination of the stability of the static sliding nonlinear equilibrium of the system. Instability is then characterized by the notion of unstable modes: an eigenvalue with a positive real part indicates the presence of an instability, and the associated eigenvector defines the unstable mode shape of the system around the equilibrium point. Vibrations are then generated by considering a small perturbation of the sliding equilibrium position that leads to an increase in the oscillation levels. This CEA method is relatively easy to implement, and results are easy to analyse since the unstable modes identified by the method are potentially the modes present in the nonlinear dynamic response. However, this method is often over-or under-predictive [9,10]. It means that all the unstable modes identified by the CEA are not necessarily present in the nonlinear dynamic response and some frequencies of the dynamic response might not be identified by the CEA. Moreover, one of the major drawbacks is that the nonlinear dynamic behaviour of the system (i.e. the self-sustained vibrations) cannot be predicted by this approach.
In order to cope with the limitations of the CEA, it is necessary to determine the dynamic response and the associated vibratory levels with other numerical strategies complementary to the CEA that allows for the obtaining of the transient and/or the steady-state regime. Two classes of strategies exist to do that, namely the temporal integration and nonlinear methods to predict self-sustaining vibrations. The temporal integration solves for each time step the solution of the nonlinear dynamic problem. Thus, one of the major advantages of this approach is that no assump-tions are made on the form of the nonlinear solution. Moreover, this method makes it possible to predict the transitory regimes as well as the steady-state solution if it exists [11][12][13][14][15]. Thus, it is possible to analyze the history of the system dynamic and to identify the different modal interactions that occurred during the transient part. However, the numerical cost is often prohibitive and performing a study for different initial conditions is not possible in an industrial context. Since performing temporal integration is not a reliable solution because of the too high numerical cost, numerical strategies have been developed to predict the steadystate regime. It is worth reminding that the frequencies of the self-sustained vibrations are unknown and so remain an unknown of the problem. A first approach consisted in an extension of the classical harmonic balance method called the Constrained Harmonic Balance Method (CHBM) [16]. The solution is assumed to be periodic or multi-periodic, and the self-sustained vibrations are approximated with a Fourier or a generalized Fourier series. One of the difficulties of this approach is due to the fact that the frequencies of the nonlinear solution are unknown. So a constraint has to be added to solve the nonlinear dynamic problem and to find the self-sustaining vibrations, as well as the associated frequencies. While the CHBM method gave interesting results on a simplified finite element model of an industrial brake system in [16], the initialization and the convergence of the method remain difficult to control. Moreover, the assumptions made on the form of the approximated self-sustained vibration are a major drawback. Other nonlinear approaches can also be found in the literature, such as the shooting method. This method is based on a temporal method that identifies a periodic solution of a problem by transforming it into an initial value problem [17]. For mechanical systems subjected to friction-induced vibrations, Charroyer et al. [18] proposed to iterate on specific initial conditions (based on an energetic criterion) and the fundamental period of the nonlinear system to predict the self-sustaining vibration. However, no extension to the multi-instability case exists today.
The nonlinear approaches mentioned above have undeniable advantages and allow to approximate the vibratory levels of the system for mono-or multiinstabilities. However, the a priori choice of the form of the approximate solution, the potential complexity of the numerical implementation and the high associated numerical cost make their use in an industrial context limited. For these reasons, the CEA remains the most widely used method in industry despite its limitations. Indeed, this last one does not allow to identify all the frequency contributions that are involved in the nonlinear dynamic response. Some recent developments try to give additional information to complete the CEA results in order to identify the modes involved in the squeal. These methods are aimed to be applicable on real industrial systems. One method developed by Brunetti et al. [19] defined the Modal Absorption Index (MAI) of each mode. It describes the capacity of each mode to absorb energy and to generate squeal. This method was applied on an academic system and gave good prediction results. A second strategy is the Modal Amplitude Stability Analysis (MASA) developed by Nacivet and Sinou [20] in the case of mono-instability. The method assumes that the dynamic response is driven by only one unstable mode; thus, the dynamic response is approximated by a linearization of the initial dynamic response around a nonlinear sliding equilibrium position for a given modal amplitude. The evolution of the dynamic behaviour (i.e. increase or decrease in the amplitude of the unstable mode) is studied.
The objective of the present paper is to extend the MASA method to the case where several instabilities are identified by the CEA and to cope with the coupling effects that may exist between unstable modes. This extended approach is named generalized Modal Amplitude Stability Analysis (GMASA). One of the main objectives is to demonstrate the validity of this new extension of the MASA methodology by considering an application on a phenomenological model for which it is possible to perform temporal integration to obtain rapidly reference results. The use of this simplified model also allows to provide a deep understanding of the contribution of each unstable mode to the nonlinear dynamic response of the system and so to make extensive comparisons between the GMASA and the reference solutions.
The paper is composed of two main parts. The first one consists in the presentation of the extension of the MASA method and its theoretical and numerical description. In the second part, the GMASA method is applied on an academic mechanical model subjected to friction-induced vibration. Thus, the prediction results of the GMASA method are compared to the reference results from the temporal integration and the validity of the method is demonstrated.

Presentation of the Generalised Modal Amplitude Stability Analysis
This section is devoted to the description of the modal amplitude stability analysis generalized to the multiinstability case. In a first time, the problem under study is presented and the CEA is briefly summarized since it is the starting point of the GMASA methodology. Then, the different steps of the method are described. Finally, in the last subsection, the global scheme of the method is given.

Formulation of the problem
Let us consider a mechanical system subjected to friction-induced vibrations; for example, the associated nonlinear dynamic equation can be written as: where U is the displacement vector and the dot denotes the time derivative. M, C and K are, respectively, the mass, damping and stiffness matrices. F nl is the nonlinear effort vector and represents the friction forces and contact conditions. F ext represents the external forces that are applied initially on the system, as the piston pressure for a brake system. It is worth noticing that the latter is constant and does not depend on the time. The CEA is the first step of the GMASA method, and it identifies the unstable modes to be taken into consideration initially in the GMASA. The global strategy is the following: a scanning on the modal amplitude of these modes is performed. Thus, for each amplitude level, an approximate displacements field is obtained. The divergent behaviour of this approximate solution is then studied by considering a linearization of the nonlinear efforts. During the process, new instabilities (i.e. the appearance of new unstable modes) can be identified by the GMASA. If this is the case, these new unstable modes are added in the group of unstable modes and are taken into consideration for the displacements field approximation.
A complete and detailed description of each step of the GMASA is presented and discussed in the following sections.

Stability analysis
The first step of the method consists in the identification of the different unstable modes with a CEA procedure. The methodology is reminded here.
The nonlinear static equilibrium U S of the system is determined by solving: In order to study the stability of this equilibrium, a small perturbation U is considered, hence: Since the perturbation is small, the nonlinear efforts can be linearized around the nonlinear static equilibrium U S with a first-order Taylor development. By denoting J nl the Jacobian matrix of the nonlinear efforts F nl at the point U S , it comes: By injecting the previous equation in the dynamic equation (1), it comes: Assuming a solution of the form U = e λt where λ represents the eigenvalues of the system and the matrix of the eigenvectors, system (5) becomes: with the following characteristic equation: Friction contributions bring nonsymmetric terms in the Jacobian matrix J nl , so the eigenvalues λ j = a j + iω j and the eigenvectors j are complex. ω j corresponds to the angular frequency of the mode j and a j to its real part. According to the Lyapunov theory, the asymptotic stability of the static equilibrium position of the system is obtained by considering the sign of the a j . If at least one of the real parts a j is strictly positive, then the nonlinear static equilibrium is unstable. It is important to remind that this analysis is locally valid since the matrix J nl is defined only in a neighbourhood of U S .
For the rest, it is assumed that N i unstable modes are identified with the CEA, characterized by their angular frequencies ω k for k ∈ [1, N i ]. The system is subjected to self-sustained vibrations that can be decomposed into N f incommensurable frequencies corresponding to the frequencies of the unstable modes. It is important to note that the number of unstable modes N i via the CEA is not necessarily identical to the number of unstable modes N f whose contributions are present in the vibration response. Indeed, the over-and underpredictive aspect of the CEA [9] means that some of the N i unstable modes are not part of the real dynamic response and that other modes may become unstable. The objective of the GMASA is to identify which modes are part of the real dynamic response.

Quasiperiodic solution
From Eq. (4), let us define as: represents the error made on Taylor's development at first order. The dynamic equation (5) becomes: In the following, the solution is supposed to be quasiperiodic and characterized by N f incommensurable frequencies corresponding to the frequencies of the N f unstable modes. It is worth noticing that N f can be different of the value N i as explained in the previous section. The oscillations are characterized by frequencies that are a linear combination of the fundamental incommensurable frequencies. Thus, if N h harmonics of the solution are retained, then the different angular frequencies that might be observed are of the form: k 1 ω 1 + · · · + k j ω j + · · · + k N f ω N f (10) where k j ∈ [−N h , +N h ] and j ∈ [1, N f ], and where only the positive terms are retained. Then, the dynamic response of the system can be decomposed in a generalized Fourier series: where the U C where τ = τ 1 , τ 2 , . . . , τ N f is the hyper-time variable [21]. It is 2π -periodic on each dimension τ j . Thus, it is possible to define k = k 1 , k 2 , . . . , k N f , the vector of the number of harmonics in each dimension τ j . By denoting as (.) the dot product, Eq. (11) can be written as: Instead of manipulating a single time domain t ∈ R + for the solution U(t), this transformation allows to deal with a multiple periodic time domain τ ∈ R N f + where each dimension τ j represents an incommensurable frequency of the solution. Thus, Eq. (13) is 2π -periodic in each dimension and U(τ ) is a N fdimensional time function that describes the multifrequential function U(t). This function is often denoted as torus function. For example, in the case of dimension 2, the function U(t) is biperiodic and evolves in R 3 on a torus formed by the two hyper-time parameter [τ 1 , τ 2 ], as represented in Fig. 1.
Therefore, the different equations written in the mono-frequential Fourier domain remain valid in the multi-frequential domain. The displacements U can be written as: The nonlinear efforts are also expected to be quasiperiodic and so can also be decomposed on a generalized Fourier series:

Projection on the first harmonic and linearization of
As in the methodology proposed in [20], a first-order truncation of the Fourier series is considered. In other words: where v k is a vector of 0 of size N f where the kth value is equal to 1 (i.e. only the first harmonic of each frequency is retained). Similarly, the dynamic equation projected on firstharmonic terms is: And the orthogonality of the Fourier series gives the N f following equations: From Eq. (8), the nonlinear efforts are equal to: By injecting Eq. (16) and projecting on the first harmonic, it comes: The term F nl (U S + U 1 ) = F nl (U S + N f k=1 U v k ) couples the different angular frequencies and so the different unstable modes. Cameron [22] proposed the Alternating Frequency/Time (AFT) method to compute the periodic nonlinear forces in the time domain and then to extract their exact Fourier coefficients. In the case of a hyper-time domain where the functions are 2π -periodic on each dimension, the AFT can be extended to a N f -dimensional frequency domain by using N f -dimensional FFT and 1 can be determined.
The originality of the GMASA method relies in the following linearization of each v k [20]: where K v k and C v k are the equivalent stiffness and damping matrices on the kth angular frequency, respectively. They can be computed analytically with the following method. On the one hand, for each contact element n, the normal relative displacement δ n k at the point n projected on the kth angular frequency is given by: On the other hand, the nonlinear efforts at the same contact element n in the dimension d, which can be tangential or normal, between two contact nodes are: And so, from the assumption done in Eq. (22), for each contact element n and for each direction d, these efforts must verify: The matrix K v k and C v k can be determined analytically from Eqs. (23,24,25); the resolution gives:

Subsystems definition
Hence, injecting Eq. (22) in (19), it comes: This leads to N f subsystems of the form:

and the dynamic matrices
A k are given by: Thus, each instability k is characterized by a dynamic subsystem of dynamic matrix A k . From the latter, it is possible to analyze the temporal evolution of the dynamic solution linearized around an equilibrium point for a given modal amplitude (i.e. increase or decrease in the vibratory amplitudes of the mechanical system). This analysis is performed, for each mode k, from the associated eigenvalue in the kth subsystem.
Each Y k is only related to ω k , indeed there is no contribution from modes j with j = k in Y k . Moreover, the projection on the first harmonic does not show any term in A j Y k . Hence, if the real part of the eigenvalue associated with the kth mode in the kth subsystem is positive, the kth mode brings a divergence in the global solution of the system defined by Eq. (18). If it is negative, then the associated vibrations decrease, and if it is zero, then the associated vibrations remain constant, which then corresponds to obtaining self-sustained vibrations.
If a mode j is unstable in a subsystem k, where j = k, then it is necessary to consider the subsystem j to know whether the divergence associated with the mode j is effectively translated by a divergence in the dynamic equation (18) because the contribution of a mode is only given by its own subsystem. In the case where the mode j is identified as stable with the CEA (i.e. a j = 0), it is necessary to introduce it in the GMASA methodology and to consider a GMASA modal basis of N f + 1 unstable modes (i.e. N f is equal to N i plus the number of instabilities detected by the method).

Determination of the modal amplitudes and of the displacements field
In the GMASA, the displacements U 1 are known (based on assumptions explained in the following) and depend on the modal amplitudes of the unstable modes. For each modal amplitude, the increase or decrease in the amplitudes of U 1 is studied. U 1 expressed in Eq. (16) is equivalent to: If U 1 is nonzero and known, then each K v k and C v k is known, and so each A k . So for each unstable mode k, its eigenvalue λ k is deduced from the matrix A k by performing an eigenvalue analysis on it.
The quasiperiodic response of a system defined by Eq. (5) is mainly driven by the N f identified unstable modes linearized around the equilibrium, so that the displacements can be approximated by Let us consider that N i unstable modes of mode shapes 0 k k∈ [1,N i ] and of eigenvalues λ 0 k = a 0 j + iω 0 k are identified by the CEA. The GMASA is initialized with these N i unstable modes and the associated frequencies (i.e. N f = N i ); then, for each scan step m, the modal amplitudes where each p m k corresponds to the amplitude of the mode k for the scanning index m, are evaluated. The number of unstable modes is updated each time a new instability appears as explained in Section 2.5, and the modal contribution of this new instability is added in the expression of the displacements. The GMASA is based on the assumptions that the amplitudes are assumed to be low enough to consider that the mode shapes 0 k and the angular frequencies ω 0 k remain almost constant. Thus, the angular frequencies in the generalized Fourier series are the ω 0 k , and for each mode, it comes: Hence, an approximation of the displacements U 1 is obtained by summing the contributions of the different modes: So for each modal amplitude p, the displacements U 1 are known and correspond to an approximate solution of problem (5). Since U 1 is known, is also known (see Eq. (21)) and the matrices K v k and C v k can be determined (see Eq. (26)). Finally, each matrix A k is completely known for a given modal amplitude p. The GMASA proposes to perform a scanning on this parameter to study the dynamic behaviour of the system.
For the determination of the modal amplitude p, an approach that proposes a physical description of the vibratory response is considered and consists in coming back to a temporal approach in order to establish a connection between the different modal amplitudes. In this way, if each mode is characterized by a modal amplitude p k (t m ) at an instant t m and an eigenvalue λ m k = a m k + iω m k , then the modal amplitude p k (t m+1 ) at the moment t m+1 = t m + δt is related to p k (t m ) and λ m k by: which is equivalent to: The scanning is then carried out in time instead of modal amplitude. The time step is chosen so that Eq. (32) remains valid, i.e. considering a constant increase a n k between the instants t m and t m+1 is a valid approximation. For a reminder, each λ m k is determined from the matrix A m k associated with the system of the mode k at the instant t m . In the case of an unique instability, the scanning can be realized directly with the modal amplitude instead of the time as proposed in [20]. It is worth noticing that the time defined here is not a real time but a fictional one that corresponds to the discretization of the method. Thanks to this formulation, displacements are known, and the different matrices defined previously can be determined for each level of amplitude.
For the sake of clarity, it is worth pointing out that contrary to the CHBM where the Fourier coefficients and the angular frequency are to be found and are the solution of an optimization problem, here they are imposed and the stability of the corresponding behaviour is studied.

General scheme of the GMASA procedure
In a nutshell, the general approach can be summarized as follows: • Initialization (m = 0): a CEA is performed on the initial system. The N i unstable modes are identified as well as their eigenvectors ( 0 k ) k∈[1,N i ] and their associated angular frequencies (ω 0 k ) k∈ [1,N i ] . This step corresponds to a zero modal amplitude factor p 0 = 0. An initialization with a perturbation p is considered, and N f = N i .  The analysis of the divergence of each unstable mode k is performed with the study of the evolution of its eigenvalue λ k = a k + iω k in the subsystem k. At the step m, if a m k is positive, then the vibrations associated with the mode k increase. If a m k is negative, then the associated vibrations decrease, and if a m k = 0, then the amplitudes of the oscillations are constant. The total behaviour is deduced by summing the individual contributions. It can be noted that if one or more subsystems have one or more eigenvalues with a zero real part (the others being negative), while all the other subsystems have eigenvalues with negative real parts, then the solution corresponds to self-sustained vibrations of the mechanical system subjected to friction-induced vibrations.
If during the process, in a subsystem k the real part of an eigenvalue λ j is positive (with k = j), then it is necessary to check the sign of the eigenvalue of the mode j in its own subsystem. If the mode j was not taken into consideration in the GMASA process, it means that the GMASA has identified a new potential instability and it is necessary to add it in the process presented (i.e. N f is incremented).
The main interest of this method is that it finally relies on performing only eigenvalue analysis and so can be applied on large finite element models with potential significant gain in computation time compared to other methods such as temporal integration.

Application
In this part, the GMASA method is applied on a phenomenological model subjected to friction-induced vibrations. In a first time, the model is introduced, and the CEA is performed to determine the stabil- ity behaviour of the mechanical system. Then, the GMASA method is applied and its results are compared to the results of the temporal integration.

Mechanical system under study
The system under study is a phenomenological model, which makes it possible to study stability and nonlinear vibrations for mechanical systems subjected to frictioninduced vibration. This four-degree-of-freedom model is represented in Fig. 2. It corresponds to an extension of the two-degree-of-freedom model proposed by Hulten [6,23]. In this study, an extension of this minimal system is considered in order to investigate the case of multi-instabilities. It has already been used in [24] to study the nonlinear behaviour of the system.
The model has two masses m 1 and m 2 held against moving bands, as displayed in Fig. 2. Each contact between masses and bands is modelled by a plate supported by springs and damping. Each spring k 11 , k 12 , k 21 and k 22 is a nonlinear stiffness defined so that the efforts are equal to k i j u + k i j,nl u 3 for the stiffness k i j where u are the displacements. The two masses m 1 and m 2 are coupled by a stiffness k a and a damping c a . The bands and the masses are supposed to be always in contact. The friction forces are determined with a classical Coulomb's law, and the friction coefficient at the different contacts is supposed to be constant and equal to μ; the stick-slip phenomenon is not taken into consideration. The velocity of the moving bands is supposed to be constant, and the relative velocity between the bands and the masses are positive so that the tangential friction forces F t do not change. According to the Coulomb's law, the tangential friction forces F t are related to the normal forces F n by the following formula: F t = μF n . In addition, in order to get a Jacobian matrix of the nonlinear efforts which is nonzero at the sliding equilibrium position, an effort P is initially applied on the system. It is worth noticing here that the dynamic of the system studied is smooth.
Thus, the dynamic equation of the system is: The vector of the nonlinear efforts is defined as follows: The Jacobian matrix of nonlinear efforts at the point The objective of the present study is to illustrate the validity of the GMASA method to predict the stability of systems subjected to friction-induced vibrations.  For this reason, it is applied here to a phenomenological model for which the dynamic response can be determined by time integration. Thus, two configurations of the system are retained, their characteristics are given in Table 1 and an initial pressure of P = [100, −100, −100, −100] T is applied on the system.

kg) (kg) (N/m) (N/m) (N/m) (N/m) (N/m) (Ns/m) (Ns/m) (Ns/m) (Ns/m) (Ns/m) (N/m) (N/m) (N/m) (N/m)
In the rest of the study, the following strategy will be applied. First, the CEA is used to determine the eigenvalues and so the stability behaviour of the system for different values of the friction coefficient μ. Five different cases with a different number of unstable modes (predicted by the CEA) are retained. Then, for these different cases, a temporal integration is performed to get the dynamic behaviour of the system. Finally, the GMASA method is used and its results are compared to the temporal integration results.

Stability analysis and preamble of the GMASA methodology
In a first time, the CEA is used to determine the eigenvalues and the stability behaviour of the system for different values of the friction coefficient μ. The evolution of the eigenvalues is given in Fig. 3 for the two configurations. The stability behaviours are different according to the considered configuration. Indeed, the first configuration has a unique interval of instability for μ > 0.48, whereas the second configuration has two intervals of instability for μ ∈ [0.17, 0.23] and μ > 0.55. For the first configuration, two mode coupling phenomena take place: the first one at μ = 0.48 for an unstable mode detected at 5.8 Hz and the second at μ = 0.54 at 10 Hz. Concerning the second configuration, a unique mode is unstable between μ = 0.17 and μ = 0.23 at 7.4 Hz. Then, two mode coupling phenomena are observed, the first one at μ = 0.55 at 4.8 Hz and the second at μ = 0.58 at 9.9 Hz. In the following sections, different cases are studied in order to illustrate the capacity of the GMASA method to predict the evolution of the vibration behaviour of the system around its nonlinear equilibrium, in other word the prediction of the increase or the decrease in the self-excited vibrations. Five different cases with different dynamic behaviours are selected. They are summarized in Table 2; the characteristics as well as the unstable frequencies with the associated real parts computed by the CEA are given. The first two cases (i.e. Cases 1 and 2) correspond to cases where a unique instability is detected by the CEA, and the last three cases (i.e. Cases 3, 4 and 5) correspond to cases where two instabilities are identified with the CEA.

Cases 1 and 2: one instability
In a first time, the two first cases are studied. They correspond to cases where a unique instability is detected with the CEA. The corresponding parameters are given in Table 2.

Temporal integration: reference results
First, a small perturbation of the equilibrium is considered and a time integration is performed on the system to determine its dynamic response for Cases 1 and 2. The steady-state regime is reached in a few seconds, and the corresponding limit cycles associated with the different degrees of freedom (dofs) are given in Fig. 4 on the right column in green and purple, respectively.
The FFT are also determined, and the different frequencies observed are given in Table 3. In both cases, the signal is composed of a main frequency and its harmonics. The main frequency is equal to 5.18 Hz for Case 1 and to 7.57 Hz for Case 2. In those cases, only three harmonics are observed so the dynamic response is driven by a unique frequency and mostly by the first harmonic. For these cases, the frequency of the unstable mode has slightly evolved between the CEA prediction and the time integration: from 5.29 to 5.18 Hz and from 7.39 to 7.57 Hz for the first and the second cases,  respectively. A unique unstable mode is identified in the dynamic response, and only a few harmonics are involved in the nonlinear dynamic response.

Application of the GMASA method
Then, the GMASA method is applied on Cases 1 and 2. Because there is a unique unstable mode predicted by the CEA, it is possible to perform directly a modal  amplitude scanning for the identified unstable mode. Hence, the method is initialized with the mode shape and the angular frequency of the unstable mode identified with the CEA and a scanning on the modal amplitude between 0 and 1 is performed. For each value of the modal amplitude p, the equivalent nonlinear efforts and the associated dynamic matrix are determined and the eigenvalues of the latter are determined. The evolution of these eigenvalues is given in Fig. 5 for both cases. The behaviour is similar for both configurations. The eigenvalue corresponding to the unstable mode is the purple one, and the evolution of the real part of this eigenvalue allows to conclude on the divergence of the approximate dynamic solution for each amplitude level. When this real part is positive, then the approximate solution is divergent; when it is negative, then the vibratory levels decrease. If the real part is zero, then the vibratory levels remain constant and the approximate solution corresponds to a limit cycle.
For the first configuration, the real part of the eigenvalue of the unstable mode decreases and becomes null for p = 0.33. Hence, when the modal amplitude is inferior to 0.33, the real part of the associated eigenvalue is positive, and the vibratory levels increase. And when the modal amplitude p is superior to 0.33, then the real part of the associated eigenvalue is negative, and so the approximate solution has decreasing vibratory levels. Then, when p = 0.33, the real part is equal to zero and so the vibratory levels remain constant. So, if a small initial perturbation of the mode is considered (i.e. an initial low value of p), then the vibratory levels will increase since the real part is positive. As long as the real part is positive, the growth of the limit cycles will be slower and slower since the associated real part decreases when the modal amplitude increases. Finally, when the modal amplitude reaches p = 0.33, then the real part is equal to zero, which implies that the approximated solution describes a steady-state regime and that the limit cycle of the self-excited vibrations of the nonlinear mechanical system is reached. This solution is composed of a unique contribution with a frequency equal to 5.19 Hz, which is almost equal to the first frequency peak observed on the FFT of the nonlinear solution obtained from temporal integration (i.e. 5.18 Hz). These results are summarized in Table 4.
Results are similar for the second case: when p increases, the real part decreases and becomes null for p = 0.79 (see Fig. 5 and Table 4). The GMASA identifies a steady-state regime driven by one frequency at 7.54 Hz, which is very close to the first frequency peak observed on the FFT from temporal integration equal to 7.57 Hz.
In these two cases, the GMASA identifies modal amplitudes for which the real part of the eigenvalue becomes zero, i.e. an approximate steady-state solution. The latter one is characterized by a limit cycle given by Eqs. (30) and (31). These limit cycles obtained with the GMASA method can be compared to the limit cycles obtained with temporal integration. The comparison between the limit cycles is given in Fig. 4 for the four dofs and for the five cases, on the left column for GMASA and right column for temporal integration. The limit cycles corresponding to the Cases 1 and 2 are in green and purple, respectively. Moreover, the maximal amplitudes of the displacements and of the   velocities are determined and summarized in Table 6. The associated amplitudes of the main frequency peak in the FFT are given in Table 5 for each dof. For the reader comprehension, f i,GMASA and f i,Temporal refer to the resonance peak for the ith frequency predicted by the GMASA and the temporal approach, respectively. From Fig. 4, it is clear that the evolutions of the limit cycles obtained with the GMASA are similar to the evolution of the limit cycles obtained with the temporal integration. For example, by considering Case 1 for the first two dofs, the limit cycles obtained with the GMASA and the temporal integration are visually represented just by a point compared to the other limit cycles due to the very low levels of vibration; this is confirmed by the amplitudes given in Table 6. The ratio of amplitudes for x 1 between Cases 1 and 2 is equal to 0.025 (resp. 0.039) with the GMASA (resp. temporal integration). It demonstrates well that amplitudes of x 1 are negligible for Case 1 compared to Case 2 and that the GMASA estimates well this tendency. Similar analysis can be bone for the other dofs. This is also confirmed by comparing the participation of the main peak in the FFT in Table 5. For Case 1, considering x 2 and y 2 , the limit cycles are larger than those of x 1 and y 1 ; this behaviour is observed on the limit cycles obtained by temporal integration as well as on those obtained with the GMASA. Moreover, in the case of y 2 , the amplitude of displacements is greater in Case 1 than in Case 2 (see Table 6 and Fig. 4) and this difference is well predicted by the GMASA. A last remark concerns the prediction of the dynamic equilibrium position. Indeed, for both Cases 1 and 2, the limit cycles are not centred around the (0,0) point, and for x 1 , the limit cycle is shifted on the right, whereas it is shifted on the left for the other dof (see Fig. 4). This is well predicted by the GMASA since the position of the limit cycles follows the same behaviour. Considering the comparison of the associated amplitudes of the main frequency peak in the FFT (see Table 5), the amplitudes given by the GMASA are good approximation of the amplitudes obtained from the temporal integration for both Cases 1 and 2. This demonstrates that in these cases, the GMASA identifies well the contribution of the first harmonic in the dynamic response of the system. It should be noted that the FFT obtained from the limit cycles of the GMASA are composed only of one peak at the frequency identified by the GMASA, whereas the FFT obtained from temporal integration have many peaks as described in Table 3.
As a conclusion, in the case of one instability, the GMASA method identifies the frequency of the selfsustained vibrations as well as the modal amplitude of the unstable mode in the steady-state regime. Moreover, the limit cycles of the steady-state response obtained with the GMASA give an approximation of the limit cycles obtained by temporal integration. If the vibratory levels are not exactly predicted, it is still possible to compare the vibratory levels between different configurations of the system. Moreover, the contribution of the first harmonic of the main frequency is well approximated by the GMASA. The good level of prediction in these cases may be due to the fact that the dynamic response is mainly driven by one frequency (i.e. the first harmonic of the frequency of the original unstable mode) and that a low number of harmonics are involved making the truncation at the first harmonic a good approximation of the real solution.

Cases 3, 4 and 5: multi-instabilities
The GMASA method shows a good ability to approximate the dynamic behaviour of the system when one instability is present. The main objective of the following study is to undertake the effectiveness and the ability of the method to treat multi-instability cases. For this purpose, the GMASA methodology is applied on Cases 3, 4 and 5 (see Table 2 for more details on each case), where two unstable modes are identified with the CEA. The validation process is identical to the previous part: a time integration is performed to get the reference solution (i.e. the self-excited vibration of the nonlinear system). The effectiveness of the GMASA methodology is performed by comparing the GMASA results to this reference.
Since a few parameters must be tuned by the user for the GMASA method, a preliminary analysis is proposed to discuss the choice of these tuning parameters. This will be more specifically illustrated for Case 3.

Preamble and GMASA parameters tuning
As previously explained, a preliminary study is performed on Case 3 to discuss the choice of the GMASA tuning parameters. As a reminder, two unstable modes are identified with the CEA for Case 3 (see Table 2) and the two associated frequencies (resp. the real parts) are equal to 5.47 Hz and 9.27 Hz (resp. 3.3 and 0.007). It is worth noticing that the second instability has a real part that is largely inferior to the other one.
By using temporal integration, the dynamic behaviour of the system is obtained. The steady-state regime is reached in a few seconds, and the associated limit cycles are given in Fig. 4 on the right column in yellow. The FFT are determined, and the frequencies observed are given in Table 3. Compared to Case 1, the vibratory levels have largely increased (see the limit cycle amplitude about 100 times larger for x 1 in Table 6), and the dynamic behaviour is more complex since the signal is composed of two incommensurable frequencies, equal to 5.31 Hz and 10.13 Hz, their harmonics and their linear combinations (see Table 3). Twelve peaks are observed in the FFT when only four were detected in Case 1. However, the dynamic response is mainly driven by the first harmonic of each incommensurable frequency since they correspond to the main peaks on the FFT. Moreover, it appears that the first mass m 1 (i.e. x 1 and y 1 ) oscillates principally at 10.13 Hz, and the second mass m 2 (i.e. x 2 and y 2 ) oscillates at 5.31 Hz. Therefore, each mass oscillates mainly at one frequency and the contributions related to the other frequency as well as the interaction between the two masses are relatively weak.
The GMASA method is then applied on the system. In a first time, it is applied for predefined values of the modal amplitudes p 1 and p 2 associated with each unstable mode. Thus, for each amplitude pair, the dynamic matrices A 1 and A 2 of each subsystem are computed and their eigenvalues are determined. The evolution of the eigenvalues of each subsystem is given in Fig. 6. It appears that the evolutions of the eigenvalues are different in each subsystem. Indeed, in the first subsystem, the real part of the purple eigenvalue remains almost constant and equal to −1.5, whereas in the second subsystem, it increases when p 1 > 2.3. In the same way, a variation in p 1 has no impact on the frequency of this mode in the first subsystem, whereas in the second, when p 1 > 2.3, the frequency increases and is impacted by a variation in p 2 . In both subsystems, the evolution of the eigenvalues is complex. For example, the real parts of the orange and yellow modes evolve abruptly on certain zones, whereas they remain almost constant on the rest of the space. So it seems complicated to identify areas where the unstable modes will evolve and so to conclude on the increase or decrease in the vibratory levels of the solution. Moreover, this kind of results requires a large number of calculations which cannot be done easily on a system with a large number of degrees of freedom. For these reasons, it is necessary to use relations between the different p i under the form of Eq. (32) in order to follow the evolution of the system.
The use of a recurrence formula such as Eq. (32) requires the tuning of two parameters, namely the time discretization, i.e. the choice of δt, and the initial conditions on the modal amplitudes, i.e. p 0 1 and p 0 2 . The influence of these two parameters will be studied in the following.
In a first time, the influence of the initialization of the modal amplitude on the results is studied. In this way, different initializations are considered, namely . Because these initializations correspond to a perturbation of the system, they are voluntarily low. The evolutions of the different parameters that characterize the evolution of the unstable modes are given in Fig. 7. For each mode, it corresponds to the evolution of its modal amplitude p i , of the real part of its eigenvalue and of its frequency. It is important to remind that the real parts and the frequencies are determined from the eigenvalues of the dynamic matrices associated with the mode at each iteration. Moreover, at t = 0, values of real parts and frequencies are those of the CEA, in other words when p = 0. Because the problem is not contin- uous in 0, large variations may be observed. According to Fig. 7, whatever the initialization is, whether it is identical or not for the two modes, the different parameters always converge towards the same value. Indeed, the two real parts converge always to 0, and the frequency of the first mode converges to 4.76 Hz, and to 9.13 Hz for the second one. The modal amplitude of the first mode always converges to p 1 = 0.47, and to p 2 = 0.97 for the second mode. These results are summarized in Table 4. Of course, the time needed to converge to the solution is highly influenced by the initial conditions. Indeed, the initializations (0.1, 0.1) and (0.5, 0.5) lead to a faster convergence of the method. In other words, less calculations are required to get the results. However, the (0.5, 0.5) initialization is slightly different from the other: the first mode has a decreasing modal amplitude, and the evolution of the real part of the second mode has a different evolution compared to the other cases. This demonstrates that an initialization with a perturbation too far from the origin implies a different "path" on the surface given in Fig. 6. For the reader interest and to be more comprehensive on this analysis, initializations with higher values of ( p 0 1 , p 0 2 ) have been tested and systematically lead to a divergence of the method. Therefore, the choice of the initialization of the modal amplitudes is important and a compromise has to be found between a low number of iterations and the convergence of the method. For the rest of the study, an initialization equal to ( p 0 1 , p 0 2 ) = (0.1, 0.1) is considered.  Fig. 8. First, it can be noticed that the choice of a δt too large generates a divergence of the solution (see δt = 10 s in Fig. 8). This result can easily be explained by the fact that, when this interval is too large, considering a constant growth during this time interval is a rough assumption that prevents from following the evolution of the modes. For the other time intervals, a modification of the discretization has no impact on the solution. Indeed, whatever the value of δt in [0.01, 0.1, 1] s, the different parameters converge to the same value. However, between the case δt = 1 s and the case δt = 0.1 s, 100 times less calculations are needed. In conclusion, the choice of the time interval is crucial. A time interval that is too large does not allow the evolution of the different modes to be followed correctly, and a time interval that is too small will require a high number of iterations. In the rest of the study, a time interval equals 0.1 s is chosen.

Case 3
Once those parameters are determined, it is possible to get the evolutions of the eigenvalues of each subsystem for Case 3. Results are given in Fig. 9. As a reminder, the diverging behaviour or not (i.e. the increase or the decrease in the vibratory levels) of each unstable mode in the global dynamic solution is given by its behaviour in its subsystem. If a subsystem shows a divergence on another mode, it is necessary to consider the subsystem related to this other mode to conclude on its participation in the global dynamic response. In this case, in each subsystem, the real part associated with each unstable mode decreases until it stabilizes to a zero value.
Finally, the results show that the modal amplitudes evolve until the vibration levels remain constant (i.e. the associated two real parts are equal to 0). As long as both real parts are not null and that at least one of them is positive, the dynamic response approximated by the GMASA method is divergent (i.e. a growth of the vibratory levels is found) and the modal amplitudes of the unstable modes evolve. After a while, the different parameters do not longer evolve: the two real parts are equal to zero, which means the two contributions corresponding to the evolution of the initial unstable modes are involved in the nonlinear self-sustained dynamic response. These results are consistent with the reference results where the two instabilities are present (see Fig. 4b, d, f, h in yellow). Thus, the GMASA identifies a steady-state response composed of two modes: one at 4.76 Hz with a modal amplitude p 1 = 0.46, and the second at 9.13 Hz with a modal amplitude p 2 = 0.97. The solution found by the GMASA defines quasiperiodic limit cycles given in Fig. 4 in yellow, and the peak amplitudes of the FFT are given in Table 5. It is worth noticing that for the GMASA, since only the first harmonic is retained, the number of peaks corresponds to the number of unstable modes. The GMASA identifies well the evolution of the limit cycles compared to the other cases (see Fig. 4, Tables 6 and 5): for x 1 , the displacement amplitudes are about 100 times larger compared to Case 1; for x 1 and y 2 , the limit cycles are larger compared to Case 2; and for y 1 and x 2 , the limit cycles are of the same order of magnitude for Cases 2 and 3, etc. The levels of vibrations are not accurately predicted but the tendencies of the evolution of the limit cycles are well predicted and comparison between the different systems is possible. Moreover, the estimation of the participation of the first harmonic of the incommensurable frequencies in the FFT is also well approximated. Indeed, the GMASA identifies well that the contribution of the first frequency is more important for the degrees of freedom x 2 and y 2 than for x 1 and y 1 (and vice versa); see the two orders of magnitude of difference in Table 5).
In conclusion, the GMASA makes it possible to follow the contributions of the two unstable modes for different modal amplitudes and leads to decide on the participation (or not) of each unstable mode initially identified by the CEA in the final steady-state regime of the system. In this case, the estimated frequencies of the two unstable modes that contribute to the steadystate regime approximated by the GMASA method are close to those of the reference solution. In addition, the CEA highlights two instabilities, one with a real part almost null and much lower than the other. The GMASA method succeeds in identifying that the two modes are present in the dynamic response with modal amplitudes of the same order of magnitude. This demonstrates the relevance of the proposed method for obtaining the contribution and evolution of the initial unstable modes in the final self-sustained solution. On the other hand, in the present case, the GMASA method has identified an approximate solution of the steadystate regime of the system. By comparison of the limit cycles to the reference for Cases 1 and 2, it appears that the GMASA is able to predict well the growth and the evolution of the limit cycle amplitudes between the different cases and can predict well the participation of the first harmonic of each frequency.

Cases 4 and 5
Finally, two last cases are investigated. As a reminder, physical parameters for Cases 4 and 5 are given in Table 1. In both cases, the CEA identified two unstable modes (see Table 2). The major difference with Case 3 is that the associated real parts of the two unstable modes are of the same order of magnitude. The limit cycles obtained from temporal integration are given in Fig. 4 in orange and blue for Cases 4 and 5, respectively, and the frequencies of the FFT are given in Table 3. Compared to Case 3, the dynamic is more complex for both cases: several harmonics and more linear combinations of the two incommensurable frequencies are involved in the dynamic of the system (see Table 3).
The GMASA method is then applied on the system. The evolutions of the modal amplitudes, of the real part of the eigenvalues and of the frequencies are displayed in Fig. 10. Results are summarized in Table 4. The real parts of the two unstable modes converge to 0 in both cases. For Case 4 (resp. Case 5), the frequencies converge to 4.98 Hz and 10.32 Hz (resp. 4.86 Hz and 9.92 Hz), and the modal amplitudes to p 1 = 0.64 and p 2 = 1.34, respectively (resp. p 1 = 0.58 and p 2 = 1.11). Even if the real parts of the unstable modes obtained with the CEA are of the same order of magnitude, the GMASA predicts that the modal participation of the mode with the highest frequency is almost two times higher than the other one in both cases.
As in previous cases, the GMASA identifies a set of modal amplitudes that characterizes a steady-state response, corresponding to the limit cycles given in Fig. 4. The amplitudes of the limit cycles are given in Table 6, and the amplitudes of the FFT peaks of the two incommensurable frequencies that describe the dynamic are given in Table 5. As previously, the limit cycles obtained by the GMASA follow the same evolution than those obtained by temporal integration: limit cycles of Cases 4 and 5 are larger than those of the other cases, the displacement amplitudes of x 2 are of the same order of magnitude for Cases 3 and 5, the limit cycles of x 2 are very spread (2 incommensurable frequencies involved), the limit cycle of y 2 is larger for Cases 4 and 3 than for Case 5, etc. If the prediction is not exact, the global tendency is well predicted, and it is possible to compare the vibratory levels of each configuration with the vibratory levels given by the GMASA. By considering the amplitudes of the FFT peaks of the two incommensurable frequencies detected in the FFT, the GMASA estimates well the participation of the first harmonics. Indeed, the dynamic response of x 1 and y 1 is more driven by the second frequency than by the first one: see for Case 4 and for x 1 and y 1 FFT peaks of the first harmonic are equal to 5.4 10 −4 m and 1.3 10 −3 m for the first frequency and equal to 4.2 10 −3 m and 1.1 10 −1 m for the second frequency for the temporal integration. This tendency is also well predicted by the GMASA (2.2 10 −4 m and 8.6 10 −4 m versus 1.9 10 −2 m and 2.9 10 −2 m, respectively). A similar analysis can be done for the contribution of the first frequency, which contributes more to the dynamic response of x 2 and y 2 . Similar analyses can be done for Case 5.
In both cases, the approximation of the frequencies is close to the frequencies of the reference results. As in the previous case, the GMASA method predicts the contribution of two initial unstable modes in the dynamic response of the system, which is consistent with the results obtained with the temporal integration. Moreover, even if the two real parts predicted by the CEA are equal, the GMASA identifies that the modal amplitudes of the two modes are different and that the modal amplitude of the second one is about the double of the modal amplitude of the first unstable mode. So even if the dynamic behaviour of the system is more complex, the GMASA method is able to predict well the dynamic behaviour of the system. The comparison of the dynamic solution approximated by the GMASA with the results obtained by temporal integration demonstrates that the GMASA is able to give an estimation of the vibratory levels and more particularly it gives the correct global evolution of the limit cycles with regard to the different configurations of the system. Moreover, the frequencies involved in the selfsustained vibrations are also well estimated.

Conclusion
The objective of the present study is to propose a nonlinear approach complementary of the CEA. This methodology, called the Generalised Model Amplitude Stability Analysis (GMASA), consists in the extension and the generalization of the modal amplitude stability analysis from the mono-instability to the multiinstability case. It consists in considering an approximation of the displacements as a function of the modal amplitudes of the unstable modes identified by the CEA. For each level of amplitude, the evolution of the transient nonlinear vibrational responses (i.e. increase or decrease in the amplitude of a predefined approx-imate vibratory solution) is estimated and studied. This evolution of the response behaviour estimated by GMASA characterizes the potential divergence or mitigation of the nonlinear dynamic response of the system. Hence, the GMASA approach identifies the unstable modes involved in the nonlinear dynamic response of the system and gives a good approximation of the associated frequencies that are present in the dynamic response of the self-sustaining system.
Moreover, this approach can also estimate the contribution and the associated modal amplitude of each initial unstable mode (predicted by the CEA) in the final steady-state vibrational response. Indeed, the GMASA is able to predict an approximation of the vibratory levels by estimating quite well the participation of the first harmonic of frequencies of the self-sustained vibrations. The GMASA methodology allows the comparison of the vibratory levels for different configurations of the mechanical system subjected to friction-induced vibration and can be used to identify design tendencies.
The application of the GMASA methodology on a finite element model of an automotive brake system subjected to multi-instabilities will be the subject of future developments.