Performances of the double modal synthesis for the prediction of the transient self-sustained vibration and squeal noise

Abstract This paper presents an efficient reduction method for predicting the nonlinear transient and steady state squeal events in mechanical systems subjected to friction-induced vibration and noise. This proposed reduction technique is based on the Double Modal Synthesis (DMS) method that involves the use of a classical Craig & Bampton modal reduction on each substructure considering the interface surfaces associated to a condensation at the frictional interface based on complex modes. In this paper, the performances of some reduced bases based on the DMS strategy are investigated in the case of a finite element model of a simplified disc/pads system. The originality of the present work is to propose a comprehensive study on the convergence of the DMS method in order to predict not only the stability or the limit cycles of a simplified brake system but also the transient nonlinear self-excited vibrations, as well as the squeal noise. A special attention is brought to the convergence of the DMS method and more precisely the number of interfaces modes required to provide satisfactory results in regard to various criterion used to characterize the squeal.


Introduction
Mechanical instabilities due to friction are still of great interest in both academic and industrial research. This is particularly the case for friction-induced vibrations and squeal noise in the automotive industry. Despite the great number of experimental and numerical studies on the subject of friction-induced vibration and noise, the squeal phenomena is still misunderstood. Reviews on the subject of friction-induced vibration and mechanism of the brake squeal phenomenon can be found in [1][2][3][4].
If we focus more specifically on simulating and predicting brake squeal, the usual and classical numerical approach is to use nowadays a Finite Element Model (FEM) of the mechanical system under study. Such numerical studies allows to predict first the propensity of brake squeal by performing a stability analysis using the Complex Eigenmodes Analysis (CEA). Additionally, a numerical integration of the differential equations of motion can be carried to obtain the transient and stationary evolutions of the nonlinear responses of the brake system, as well as the associated squeal noise. This second phase is often challenging because it requires prohibitive calculation times and the implementation of specific numerical techniques. However, it appears essential for a better understanding of the squeal phenomenon and a more accurate and robust design of brake systems in relation to the problem of squeal noise. Indeed it has been observed that generally the stability analysis may lead to an over-estimation or an underestimation of the unstable modes observed in the nonlinear time simulation [5,6] due to the fact that the stability analysis predicts only the onset of instabilities from the mode coupling phenomenon around an equilibrium position. A common way to reduce computation times is to implement model reduction techniques [7][8][9][10][11]. The challenge is to be able to define reduced models capable of reproducing the squeal propensity as well as the non-linear self-excited vibration and squeal noise.
According to the authors' knowledge, there have been few attempts to analyze the performance of the reduced bases with respect to the prediction of transient vibration squeal events as well as the associated acoustic radiation. Therefore one objective of this paper is to address the efficiency of the generalized Double Modal Synthesis method for the prediction of self-excited vibration and squeal noise. It has been demonstrated in previous studies [10,12,13] that this numerical technique can significantly reduce the size of the original finite element model of the brake system by combining a classical modal reduction and a condensation at https://doi.org/10.1016/j.apacoust.2020.107807 0003-682X/Ó 2020 Elsevier Ltd. All rights reserved. the frictional interface via a reduced complex modal basis. However, these previous conclusions are limited to the estimation of the squeal propensity [10,12] and the prediction of limit cycles [13].
This present paper is a complement to these previous studies to address the use of the generalized DMS method in conducting rigorous numerical simulations for both the prediction of transient squeal events ans the estimation of squeal noise. The paper is organized as follows: first the finite element model of a simplified brake system and its stability analysis are presented. Secondly the Craig & Bampton technique as well as the interface reduction based on complex eigenmodes are briefly discussed. Then numerical studies based not only on the stability analysis but also on the prediction of self-excited transient and steady state vibrations, as well as the squeal noise are investigated. Convergence of the Double Modal Synthesis is tested and discussed in detail for each criterion of brake squeal prediction.

Finite Element Model of the simplified brake system and stability analysis
In this section, the Finite Element Model (FEM) of the simplified brake system under study will be first presented. Then a rapid reminder about the stability analysis is proposed.

Finite Element Model
In this study a simplified model of brake system is used. It consists of a disc and two pads. These three mechanical components are considered to be the main contributors to brake squeal. The finite element model of this brake system is depicted in Fig. 1. The structured mesh uses hexahedrons with linear interpolation. The two pads slide face to face on either side of the annular disc. Both of the interfaces consist of 220 contact nodes which are distributed equally on the friction surfaces, forming matching meshes between the disc and each pad. This choice of refinement was previously discussed in [14,15] with a similar model using only one pad. The back plate of the lower pad (pad 2) is clamped as well as the inner surface of the disc. The four corners of the back plate of the upper pad (pad 1) are connected to a master node which is only allowed to move in the out of plane direction, its rotations are clamped as well. The hydraulic pressure is directly applied on the back plate of the upper pad. Details on the finite element model, geometry and materials are given in Table 1.
The reaction forces at the two interfaces are the sum of the normal and tangential reactions. The contact forces on the two interfaces follow a cubic law on the normal direction, as previously suggested by several studies [14][15][16], and possible loss of contact between the matching nodes on the disc and the nodes on the corresponding pad. So the normal reaction is chosen under the following form for the i th couple of matching nodes on any interface, where d i is the relative normal displacement between one pad and the disc: where k L and k NL are the linear and nonlinear stiffnesses at the friction interfaces (k L ¼ 10 3 N:m À1 and k NL ¼ 10 12 N:m À3 ). For the interested reader, this basic formulation choice allows to approximate the first and the third orders of pad compression curves obtained from experimental tests [5]. Such a definition for the contact model belongs to the class of penalty methods and ensures a regularized contact. The normal reaction forces on each pad is defined as follow: The tangential reaction direction are given by considering the classical Coulomb's law with a constant friction coefficient l and a steady sliding motion without any stick-slip event: where t i ! , that defines the local orthoradial in-plane direction, can be defined using the polar coordinates: Hence the global nonlinear force F NL is: The set of equations describing the dynamic brake system can be written by: where X is the generalized displacements vector while the dot denotes derivative with respect to time. F ext corresponds to the  hydraulic pressure on the back plate of the pad 1. M; C and K are the mass, damping and stiffness matrices, respectively. By concatenating diagonally the matrices of the disc, pad 1 and pad 2, expressions of the matrices M and K are given by: The damping matrix C is estimated by considering a Rayleigh's damping under the form of C ¼ aK þ bM. The coefficients a and b are chosen to have approximately 0:5% of damping for the two fre- . As a reminder, the use of the Rayleigh's damping formulation tends to overdamped modes outside the interval of the two chosen frequencies f 1 and f 2 , but remains valid over a frequency interval greater than in the case of using a hysteresis damping. However, it has the advantage of being easily implemented and it allows to have a realistic equivalent damping that remains valid over a frequency interval greater than in the case of using a hysteresis damping.

Stability analysis
The usual approach to determine system's stability is the Complex Eigenmodes Analysis (CEA). This study is performed locally around a sliding equilibrium position X eq , which is defined by Assuming small perturbations dX around X eq such as X ¼ X eq þ dX the following linearized equations of motion are obtained: where J F NL X eq À Á is the Jacobian of the nonlinear reactions F NL in X eq . It can be noted that an analytic expression of the Jacobian can be easily calculated via equations of the reaction forces at the interfaces. The complex eigenvalue analysis then provides the stability of the equilibrium point X eq . The eigenvalue problem associated is given by: where (k k ; X k ) are respectively the eigenvalues and eigenvectors of the discretized problem, with complex eigenvalues k k ¼ a k þ ib k . The imaginary part b k defines the pulsation of the mode while the real part a k is related to the stability of the system. If at least one eigenmode have a strictly positive real part (9k; a k > 0), the system is locally unstable around the equilibrium point X eq .

Reduction methods
The initial FEM system has 50000 degrees of freedom (dof), which can lead to prohibitive calculation costs. Therefore reduction techniques are applied on the FEM system to condense the number of dof. In the following the two successive steps of reduction are briefly presented: first the Craig & Bampton technique used to condense each substructure; and secondly the generalized Double Modal Synthesis (DMS), condensing the remaining dof at the two interfaces between the disc and the two pads, using a complex eigenmodes basis.

Craig & Bampton reduction
The Craig & Bampton technique ( [17]) consists in building a reduced model for each substructure (i.e. the disc, pad 1 and pad 2), using constraint modes (W n ) and fixed interface modes (U n ), with n taking the values of d; p1 or p2 respectively for the disc, pad 1 or pad 2. The validity of this approach and its convergence has already been discussed in [10].
The criterion proposed to choose the number of fixed interface modes is based on the maximum frequency of study f max (in the present study f max ¼ 20 kHz) and a coefficient c as such as, for each substructure we keep the fixed interface modes which frequency is less than f c ¼ cf max . The internal dof (u d int ; u p1 int ; u p2 int ) are expressed in function of the generalized dof (n d ; n p1 ; n p2 , corresponding to the fixed interface modes U n ) and interface dof (u d I ; u p1 I ; u p2 I ): This allows to define the following transfer matrix T between the complete model and the condensed Craig & Bampton model: Finally matrices of the equation of motion given in Eqs. (5) and (7) can be expressed in the Craig & Bampton basis as follows: In conclusion this first reduction allows the condensation of all internal dof of each substructure by keeping the physical dof that are located on the two frictional interfaces. So the aim of the second reduction based on the generalized Double Modal Synthesis will be to condense all these interfaces degrees of freedom (dof) at each of those nodes.

Interface reduction based on complex eigenmodes
The following condensation is based on interface reduction. Previous studies showed the efficiency of this method [12,13]. The aim is to compute the complex eigenmodes of the assembled structures only at the interfaces, leading to left and right eigenmodes of the interface system dof. Thus we need to extract the remaining dof (i.e. at the interfaces nodes). By introducing the transfer matrix such as it is possible to define the following reduced mass and stiffness matrices M $ and K $ given by As K tot CB depends on l, a DMS reduction for each equilibrium point and so for each set of parameters chosen has to be perform unlike the Craig & Bampton reduction which is valid for any set of parameters. The i th right and left eigenmodes X i Br and X i Bl are then computed, sorted by order of increasing frequencies, and we define the matrices u r and u l which contain the first N DMS eigenmodes: The size of this step of reduction is set by the number of interfaces modes kept N DMS , associated with the degrees of freedom g I . Using the right and left matrices of transfer T r and T l defined by and assuming a rearrangement of the dof the reduced matrices for each contribution of the FEM brake system are given by Finally the stability of the condensed system can be studied with the same equation as Eq. (8) with the adapted matrices issued from the Craig & Bampton reduction and the Double Modal Synthesis.

Stability analysis
In this section results for the stability of the reference FEM model are firstly presented. Then the effectiveness of the generalized Double Modal Synthesis for predicting the stability of the original system is investigated an discussed.

Stability of the reference model
The study of stability of the reference system is carried out by computing its complex eigenvalues, for a set of a varying parameter. In the present study, the evolution of the friction coefficient l for 501 samples equally distributed between 0; 1 ½ is chosen. Then eigenvalues of the system for each value of l are computed. Within the range of study for l ¼ 0; 1 ½ there are two unstable modes: the first one is detected at l 1 ¼ 0:67 and the second one at l 2 ¼ 0:76.
The evolution of all the real parts of eigenvalues with the coefficient of friction are illustrated in Fig. 2(a). The coalescence pattern for the two unstable modes are also illustrated in Fig. 2(b). Moreover, Figs. 3 illustrate the modes shape. It is observed that both the disc and the pads are involved in the unstable mode shapes.

Definition of the criteria
In this section we propose to undertake more precisely a convergence study for the Double Modal Synthesis. A preliminary study only based on the convergence associated with the Craig & Bampton reduction has been performed. It was found that a Craig & Bampton reduction with c ¼ 1:5 can be considered as an efficient configuration in terms of compromise between computation time, size of the reduced system and the accuracy of results. This classical analysis is not discussed in the present paper for the sake of brevity. The Double Modal Synthesis is therefore applied to this model reduced with the Craig & Bampton technique with c ¼ 1:5 and this latter serves as a reference to evaluate the convergence and effectiveness of the second step of reduction.
In order to quantify the quality of the convergence for DMS reduction on the CEA analysis, two criteria that are based on the relative error criterion between imaginary and real parts for all frequencies in the frequency range of interest are firstly proposed. Error on the real parts dr and error on the imaginary parts dx are defined by: where N tot denotes the number of eigenvalues in the frequency , respectively) defines the real part (imaginary part, respectively) of the k th mode for the ''reference" and a k (b k , respectively) is the approximate real part (imaginary part, respectively) of the k th mode calculated via the reduction technique. Two additional criteria dr un and dx un that are only based on the relative error criterion between imaginary and real parts of the unstable modes in the frequency range of interest are also undertaken. This two criteria are given by where N un defines the number of unstable modes. In the present study, N un is equal to 0, 1 or 2 depending on the value of the coefficient of friction. In this case, a ref , respectively) defines the real part (imaginary part, respectively) of the k th unstable mode for the "'reference"' and a k (b k , respectively) is the approximate real part (imaginary part, respectively) of the k th unstable mode calculated via the reduction technique. Finally the following adapted MAC-criterion defined by the matrix M Mac (see [18][19][20]) is used to compare the eigenmodes shapes of the reference X ref i and those issued from a condensed system X j : This matrix M Mac is composed of coefficients between 0 and 1, two modes shapes close to each other tend to have their coefficient near 1, whereas a lower coefficient means that the two modes have different shapes for this projection. When comparing the reference to a condensed model, we can assume that the convergence is good when the diagonal of this MAC-criterion matrix is filled with coefficients close to 1.

Results
Considering the previous results given in Section 4.1, four different cases are undertaken in order to validate the efficiency the DMS reduction: l ¼ 0 : the case with no friction; l ¼ 0:7 : a case with one unstable mode predicted by the stability analysis; l ¼ 0:8 : a case with two unstable modes with real parts very close; l ¼ 1 : a case with two unstable modes with one whose real part is much greater than the other.
Figs. 4 give the evolution of the two criteria dr and dx versus the number of modes retained for the DMS reduction. It is illustrated that the proposed reduction technique based on the generalized DMS converges very well. It appears that for the four cases (l ¼ 0; l ¼ 0:7; l ¼ 0:8 and l ¼ 1) a number of DMS modes above 150 is sufficient to allow a lower error on the real parts dr (on the imaginary parts dx, respectively) than 1% (0:3%, respectively).
Then Figs. 5 give the results by considering only the relative error criterion between imaginary and real parts of the unstable modes (i.e. dr un and dx un ). Again it is clear that convergence is achieved for a very small number of DMS modes. Moreover, it can be observed that the DMS convergence is faster than in the previous case (i.e. the DMS reduction versus all the modes in the frequency range of interest), which is consistent with the fact that the unstable modes are low frequency modes with respect to the frequency range selected for the study.     To be noted that comparison based on computing resources between the Craig & Bampton reduction and the Double Modal Synthesis with real eigenmodes (i.e. considering interface modes computed with l ¼ 0) has been previously proposed and discussed in [10].

Nonlinear self-sustained vibrations
As previously explained in [5], a time simulation has to be performed in addition to the stability analysis in order to predict more effectively friction-induced vibration and squeal noise. This is due to the fact that the stability focuses only on a notion of local behavior around an operating point and so it may lead to an underestimation or an over-estimation of the unstable modes observed in the nonlinear time simulation. Thus it seems essential to investigate the effectiveness of the generalized Double Modal Synthesis in order to predict the nonlinear behavior of the mechanical system under study subjected to friction-induced vibration.
In order to achieve such an objective, the prediction of the vibratory behavior of the brake system is performed by applying a small perturbation around the sliding equilibrium point and by using the time integrator based on the h-method (see the following contributions for more details [21,22]). It can be noted that the choice of these initial conditions (i.e. the perturbation around the equilibrium point) can have an impact on the prediction of the nonlinear self-sustained dynamic behavior of the mechanical system under study [22].
In the following, the efficiency of the reduction technique is first analyzed by comparing results based on complex interfaces modes with those from the reference model based on the finite element model without any interface condensation. Several quantities and types of results will be analyzed to discuss the relevance of the DMS reduction: the time series (transient dynamics and selfexcited vibrations), the spectral content by using short time Fourier transform, the limit cycles and the nonlinear signature of the signal in the frequency domain. Additionally a specific attention and a more detailed analysis of the contribution of the complex modes will be undertaken.
In this section, two cases will be more particularly treated and presented: the first one considers a coefficient of friction l ¼ 0:8. This case will lead to a periodic non-linear behaviour (even if the stability analysis has potentially detected two unstable modes); the second one considers a coefficient of friction l ¼ 1. This case will lead to a quasi-periodic non-linear behaviour.
Two specific points have been selected to illustrate all the results concerning the nonlinear analysis of friction-induced vibration. The first one, denoted Point1, corresponds to a node on the disc surface (outside of the interfaces between pads and disc). The second one, denoted Point2, corresponds to a node on the frictional interface of the upper pad. The geometrical position of Point1 and Point2 on the brake system under study is illustrated in Fig. 7. For each node, the results in the three directions will be presented: the normal direction to the surface of the disc (i.e. Z direction) and the two associated tangential directions (i.e. X and Y directions). and Point2 can be divided into a succession of three main phases: the first one (for t ¼ 0; 0:3 ½ s) corresponds to a classical exponential increase of the solution from the perturbed equilibrium. Then a transient nonlinear behavior of the self-excited system appears (for t ¼ 0:3; 0:7 ½ s) with small increases or decreases of the nonlinear amplitudes depending on whether looking at normal or tangential displacements. During this second phase, short time periods of stabilization of vibration levels can also be observed (see for example Point1 in Y direction on Fig. 8(a) or Point2 in Z direction on Fig. 8(b)). By analyzing and comparing more precisely the results of the different DMS reductions, it appears that small time offsets of the appearance of three successive phases are visible. For example, in the case of l ¼ 0:8, a significant time advance is detectable in the appearance of the three phases for the DMS reductions of small size (i.e. 50 or 150 interfaces modes, as illustrated in Figs. 8(c-f)). By increasing the size of the reduced system, this time delay disappears (see Figs. 8(g-j)). The temporal evolutions of the nonlinear vibrations for the reference model and the reduced systems are identical. Considering the results for l ¼ 1 similar conclusions can be stated. Moreover it is also clearly observed that using a reduction basis of small size leads to an acceleration of the temporal transient behavior. For example, the time interval of the second phase for Point1 in the Y direction (or for Point2 in the Z direction) is drastically reduced if a DMS reduction with 50 interface modes is used, whereas this time interval is exactly reproduced by applying a DMS reduction with 1000 interface modes. In conclusion, the DMS reduction allows to reproduce the transient non-linear behavior. Despite time offsets for the appearance of the transient successive phases, using a small size of DMS reduc-tion is efficient in order to obtain an accurate estimate of the transient behavior of the mechanical system under study. In all cases, the sizes of the DMS reduction bases are advantageous compared  to the size of the reference system, which shows the interest of using such reduction technique for the prediction of transient self-sustained vibration of brake system.

Spectral content -transient and stationary self-excited vibrations
In this section, the nonlinear signature during the transient and stationary regimes of self-excited vibrations will be analyzed.
First of all short time Fourier transform is used to visualize faster time-varying signals and the spectral content during the squeal event. Results are presented in Figs. 10 and 11 for the reference system. Spectrograms obtained for the different DMS reduction are not illustrated in this section for brevity purposes as the conclusions are similar to those stated above: apart from the fact that using DMS reduction of small size induces a time shift is on the evolution of the frequency signature over time (as previously discussed for the time evolution), the evolution of the spectral content is the same for the reference model and the reduced system. Therefore the nonlinear signature of the signal is briefly discussed. For the case l ¼ 0:8 the nonlinear signature is composed by the second unstable mode (i.e. f 2 ¼ 2205 Hz) and its second harmonic component (i.e. 2f 2 ¼ 4410 Hz) in the pseudo-stabilized phase (i.e. the second phase). Then, during the shift to the final stabilized regime (i.e. the third phase), a transition from the second unstable mode f 2 to the first unstable mode f 1 appears. The final signature is generally composed by the fundamental frequency of the unstable mode (i.e. f 1 ¼ 857 Hz) and the contribution of its second harmonic component. In the normal direction at the disc/pads interface, the contributions of the 3Â harmonics are also very visible (see Fig. 10(b) in Z direction). During the transition between second and third phases (for t ¼ 0:4; 0:7 ½ s), a brief presence of linear combinations of the two unstable modes can also be observed, like a contribution at 1347 Hz that corresponds to Àf 1 þ f 2 (see Fig. 10(b) in Z direction for example). For the case l ¼ 1 the nonlinear signature is also composed by the second unstable mode (i.e. f 2 ¼ 2194 Hz) and its second harmonic component (i.e. 2f 2 ¼ 4410 Hz) during the second phase. For Point1 it is also observed that the contribution of the 2Â harmonic components is greater than contribution of its fundamental frequency in the tangential Y direction (see Fig. 11(a)). For Point2 the major contribution in this second phase corresponds to the fundamental frequency f 2 regardless of the direction considered. During the final phase the nonlinear signature of the stationary self-excited vibration is more complex: not only the fundamental frequencies of the first and second unstable modes are visible (with f 1 ¼ 856 Hz and f 2 ¼ 2194 Hz) but also harmonic components (2f 1 and 2f 2 ) as well as harmonic combinations such as Àf 1 þ f 2 .
Finally Table 2 gives the frequency of the first and second unstable modes for the stability analysis (CEA), and the fundamental frequency of the first and second modes from the spectral content for the reference and the reduced models. It clearly appears that the fundamental frequencies predicted by the reduced models are very close to those of the reference even if the DMS basis with only 50 interface modes. It is also observed that these fundamental frequencies are slightly different from the one predicted by the stability analysis.
Finally the limit cycles and the associated Fast Fourier Transform (FFT) on the stationary self-excited vibrations are depicted in Fig. 12 (Fig. 13, respectively) at Point1 (Point2 respectively) for l ¼ 0:8. In the same way the limit cycles and FFT are given in Figs. 14 and 15 for l ¼ 1.
For the case l ¼ 0:8 results on vibratory levels and frequency content via the reduction technique are in perfect agreement with the reference, even if the number of interface modes retained are very small. Indeed, it appears that the smallest reduced base tested (i.e. DMS with only 50 interface modes) gives very similar and satisfactory results compared to the reference model. The reduced model is able to reproduce the same established regime than the original system which illustrates again the efficiency and robustness of the proposed reduction method. All the associated FFT clearly indicate that the stationary regime is periodic with the contribution of the fundamental frequency around 857 Hz as well as its harmonics.
For the case l ¼ 1 results based on DMS reductions are also in agreement with the reference for the prediction of the quasiperiodic limit cycles. In order to facilitate the reading of the graphical results, the comparison is only proposed for the DMS reduction with 150 interface modes. As previously explained the FFT content of the nonlinear quasi-periodic response is complex with the presence of not only the fundamental frequencies of the first and second unstable modes but also harmonic components as well as harmonic combinations. Even if the nonlinear signature is highly complex, the reduced model is able to reproduce the same quasiperiodic established regime with frequency contents very similar to those of the reference. The main frequency contributions of  Table 3.
In this case of quasi-periodic behaviour a lot more frequencies are present in the spectrum than for the periodic case, with linear combinations of the fundamental frequencies of the two unstable modes.

Modal contributions in transient evolutions
As previously proposed in [23,22,24] the nonlinear solution can be expanded on the basis of the complex modes in order to illustrate the energy contribution of each unstable mode during the transient response and to detect the contribution of one or more complex modes in the transient nonlinear response. This complex modal projection can also be used to validate the efficiency of the proposed reduced technique. In the following, this complex modal projection will be performed for different sizes of the Double Modal Synthesis reduction in order to illustrate the effects of advance or delay in the evolution of modal participations because of the use of reduced model. For a transient response Y t ð Þ in state-space variables (i.e. Y ¼ _ X X À Á t ), the evolution of the modal participation b i t ð Þ of the i th unstable mode, which can be seen as projections along the reference unstable modes of the time signal, can be expressed in the original physical basis: where Y ref ;ins i define the modes shapes in state space. The matrix A $ represents the computation of the total energy of the motion along a particular direction, here along the deformations of the two unstable modes. The original modal participations have to be processed with a low-pass filter to obtain an average evolution < b i t ð Þ > from Eq. 27. Figs. 16 and 17 give the evolution and the contribution of the two unstable complex modes during the squeal event. For example it can be observed that the three main zones previously detected for l ¼ 0:8 are easily detected again (see the reference black line in Fig. 16): t ¼ 0; 0:3 ½ s: the initial increase of the solution around the sliding equilibrium point is mainly governed by the growth of the participation of the second unstable mode; t ¼ 0:3; 0:7 ½ s: the second phase is first characterized by a stabilization of the contribution of the second unstable mode and an increase of the participation of the first unstable mode(see t ¼ 0:3; 0:6 ½ s). In a second time (for t ¼ 0:6; 0:7 ½ s) the contribution of the second unstable mode decreases whereas the participation of the first unstable mode remains constant;  Table 3 Frequencies and corresponding combinations for the case l ¼ 1 on the limit cycles. Àf 1 þ f 2 1714 2f 1 2223 3f 2 t > 0:9 s: the final phase is characterized by a stabilization of the contribution of the first and the second unstable modes. It leads to the appearance of the nonlinear stationary solution that is mainly governed by the participation of the first mode.
For l ¼ 1 the same analysis can be performed for the evolution of the modal contributions of the two unstable modes: t ¼ 0; 0:1 ½ s: the initial growth of the nonlinear solution is due to the growing evolution of the second unstable mode; t ¼ 0:1; 0:4 ½ s: a stabilization of the second unstable mode is observed while the contribution of the second mode increases. For t > 0:3 s the contribution of the second unstable mode decreases rapidly. t > 0:4 s: the quasi-periodic solution is composed by a non-yet stabilized oscillation of the contribution of the first and second modes.
In order to evaluate the performance of the reduction on accurately reproducing the evolution of the modal contribution, results for the DMS reduction with 50, 150, 500 and 1000 interface modes are provided in Figs. 16 and 17. As previously explained in Section 5.1, the four reduction bases tested allow to reach the same final nonlinear solution with exactly the same contributions for the first and the second modes. However it can be seen that a time shift on changes in growth or decay of each contribution is clearly observable and noticeable in the case of a reduced base of small size. When the size of the reduced base is increased, this shift diminishes and all the evolution of the contributions of the two modes become coherent with those of the original system. It demonstrates the convergence of the Double Modal Synthesis reduction and its capacity to reproduce exactly the same squeal event and its complexity during the transient nonlinear behavior both in terms of level and in terms of vibrational signature.
The values of the converged modal participations are depicted in Table 4. Even if the set of results are very close, we can observe the convergence of those values with the increase of the reduction size of the DMS. Those levels are also different between the two cases of study: l ¼ 0:8 and l ¼ 1, the periodic case shows a greater ratio between the levels of first unstable mode with the second, whereas the quasi-periodic case shows levels closer to each other.
In conclusion, in order to have a good compromise between precision and computational costs, the optimum number of interfaces modes has to be selected in connection with the main desired objective of the numerical simulation, that is the estimation of the final regime or the prediction of all the transient squeal events.

Acoustic radiation
The prediction of acoustic radiated noise from friction-induced vibration is a recent topic of interest [14,15,[25][26][27]. A comprehensive review on acoustic friction phenomena can be found in [28]. The purpose of this section is to investigate the effectiveness of the DMS reduction method on predicting both stationary and transient squeal noises.

Modeling of the acoustic problem using the Boundary Element Method
In this section the Boundary Element Method (BEM) used for the prediction of acoustic noise is discussed. This method allows to compute an approximation of the acoustic radiated by a surface S, which consists in the resolution of the Kirchhoff-Helmholtz surface integral, giving the field pressure P at a point of coordinate r, such as: with a parameter equal to 1 when strictly outside of the surface and taking a value in 0; 1 ½ , when r 2 S. G defines the classical Green function, which also depends on the frequency. Therefore, the Kirchhoff-Helmholtz integral has to be solved for each considered frequency. More details about the BEM can be found in [29,30].
The calculation of the integral is a two-step process: first the pressure at the surface of the system P s has to be determined. The pressure field P at any point outside of the system can then be computed. The surface of the system S is discretized in S ¼ [ l S l using the wrapping skin of the FE Model presented in Section 2.1, as previously suggested in [25]. The normal velocity field needed for BEM approach is provided by temporal calculation of self-excited vibration previously discussed in Section 5. These BEM calculations are conducted through the OpenBEM software [31] for a number of frequencies. The choice of these frequencies is performed by calculating a Fast Fourier Transform (FFT) of the normal velocity field _ X n at each point of the mesh, and then selecting manually the K main frequency contributions Þ2Z 2 ). The velocity field _ X n can then be estimated as follows: where a i BEM and b i BEM define the Fourier vectors at each frequency f i BEM . Due to the fact that the numerical resources necessary to compute the external pressure responses for all the frequencies selected in f BEM can rapidly become prohibitive, a selection on the optimal sample of these frequencies (denoted f BEM;opt ) may be used to reduce the number of calculations to be performed.This method have been successfully tested in [14]. It is based on the convergence of the mean square velocity of the surface velocity field _ X n with the K opt different major components of f BEM . For more details the interested reader is referred to [14].
For each component _ X i n (and so a specific frequency), the BEM method first computes the associated surface pressure P i s assuming that the velocities of the solid and the fluid are equal at the interface. Then, the associated pressure field P i on any exterior set of observation points can then be computed. Finally, squeal noise that corresponds to the acoustic field P ext is defined by where K opt is the optimal number of selected frequencies taken from the frequency contributions f BEM . Of course the pressure level in dB can be calculated such as where P ref corresponds to the classical reference value used. (i.e. For the present study, the calculation of squeal noise is proposed by considering two main surfaces of observation: the first one is a square plane, parallel to the disc and placed at a distance Z square from the upper surface of the disc and of a length of L square . The second one is a sphere of radius R sphere centered around the brake system. Considering these surfaces of observation, the results of the acoustic calculation are discrete patterns, approaching the acoustic field on a set of points. From a practical point of view, the acoustic field calculations needs an evaluation of the surface velocity spectrum. This velocity field comes from the previous calculation of the nonlinear self-excited vibrations. The FFT is then performed on the surface velocity field along a short period of time (between 0:05 s and 0:1 s) for which the signal is considered to have a slow transient evolution. Table 4 Values of the stabilized modal participations for both cases l ¼ 0: 8  Moreover, the squeal noise is calculated for l ¼ 0:8 and l ¼ 1 during the stationary regime as well as the transient vibratory regime for l ¼ 0:8 and three specific time intervals. The various time intervals selected to calculate the acoustic fields are given in Table 5 and can be described as follows: Phase I: the pseudo-established vibrational regime, driven by the second unstable mode f 2 ;  By analyzing all presented results, some general observations may be brought: the acoustic patterns estimated for each time interval are different (see and compare the reference for each time interval). This characterizes the evolution of the squeal noise during the squeal event; the ability of the Craig & Bampton reduction to reproduce the radiated noise is clearly observable for the four time intervals. Even for small reduction basis, the radiated fields are very close to the reference for both the horizontal plane square and the sphere. For example, the complexity of the color pattern for Phase I of the sphere (first line of 20) is very accurately replicated; the DMS technique is efficient to approximate the acoustic patterns not only for the final vibrational regime (i.e. Phase III for l ¼ 0:8 and l ¼ 1), but also the transient vibrational regime (i.e. Phases I and II for l ¼ 0:8); from a general point of view, the acoustic patterns predicted by the different DMS reductions (for a specific phase and l) are sensitively close. Increasing the DMS size results in the convergence of the calculated acoustic field compared to the reference. Moreover it can be noted that if the surface velocity field converges, the acoustic calculation tends to do so; the effects of the DMS size can be observed on the temporal shifts of the acoustic patterns (see Table 5). Of course, this observation is in total agreement with the previous results on the nonlinear vibrations (see Section 5.1).

Conclusion
In this paper, an extension of the Double Modal Synthesis technique [12] has been proposed and used to predict the nonlinear transient self-substantial vibrations as well as the squeal noise.
It was highlighted that the use of an efficient DMS for the estimation of stability analysis and the steady state limit cycles may be insufficient to correctly approximate the transient squeal events due to phenomena of delay or advance on the evolution of the modal contributions during transient squeal events. Therefore the size of the DMS basis needs to be increased in order to achieve proper convergence to predict transient squeal behaviors. Also it can be said that the definition of the most optimal reduced basis size in terms of calculation time and accuracy of results therefore requires special attention and depends on the criterion used to characterize the squeal, i.e. the detection of only squeal occurrences via stability analysis, the calculation of limit cycles, the pre- Table 5 Time intervals (in second as unit of time) of the transient signal used to compute the FFT for the acoustic calculations. Considering results for a simple disc/pads model, it appears that the DMS technique could be very interesting in a future context of robust design to eliminate or reduce brake squeal, since it provides fast convergence and requires only a recalculation of the transfer matrices for each change of physical parameters like the friction coefficient. However one of the major issues lies in the ability to determine an optimized cost/precision ratio which is not obvious at first glance since it depends directly on the outputs used to characterize squeal.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.