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

G. Corradi, Jean-Jacques Sinou, S. Besset

► To cite this version:

G. Corradi, Jean-Jacques Sinou, S. Besset. Performances of the double modal synthesis for the prediction of the transient self-sustained vibration and squeal noise. Applied Acoustics, 2021, 175, pp.107807. 10.1016/j.apacoust.2020.107807. hal-03260508

HAL Id: hal-03260508

https://hal.science/hal-03260508

Submitted on 15 Jun 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
<figure><img src="image_2.png" /></figure>

ELSEVIER

<figure><img src="image_3.png" /></figure>

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

G. Corradi<sup>a</sup>, J.-J. Sinou<sup>a,b,*</sup>, S. Besset<sup>a</sup>

<figure><img src="image_4.png" /></figure>

## ARTICLE INFO

**Article history:**

Received 15 July 2019

Received in revised form 28 August 2020

Accepted 24 November 2020

Available online 17 December 2020

**Keywords:**

Brake

Squeal noise

Friction-induced vibration

Non-linear vibration

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

© 2020 Elsevier Ltd. All rights reserved.

## 1. 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–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–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

<footnote><sup>a</sup>Laboratoire de Tribologie et Dynamique des Systèmes, UMR CNRS 5513, École Centrale de Lyon, 36 avenue Guy de Collongue, 69134 Écully Cedex, France</footnote>

<footnote><sup>b</sup>Institut Universitaire de France, 75005 Paris, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail addresses: gregoire.corradi@ec-lyon.fr (G. Corradi), jean-jacques.sinou@ec-lyon.fr (J.-J. Sinou), sebastien.besset@ec-lyon.fr (S. Besset).</footnote>
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 and 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.

## 2. 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.

### 2.1. 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–16], and possible loss of contact

<table><caption>**Table 1**
Model parameters.</caption><thead><tr><th>Parameter</th><th>Disc</th><th>Pads</th></tr></thead><tbody><tr><td>Thickness (Z direction, m)</td><td>0.02</td><td>0.0128</td></tr><tr><td>Outer radius (m)</td><td>0.151</td><td>0.147</td></tr><tr><td>Inner radius (m)</td><td>0.034</td><td>0.091</td></tr><tr><td>Young modulus E (Pa)</td><td>1.25 10<sup>11</sup></td><td>2 10<sup>9</sup></td></tr><tr><td>Poisson's ratio</td><td>0.3</td><td>0.1</td></tr><tr><td>Density (kg.m<sup>-3</sup>)</td><td>7200</td><td>2500</td></tr><tr><td>Angular pad (degree)</td><td>-</td><td>50</td></tr><tr><td>Number of nodes in thickness</td><td>4</td><td>3</td></tr><tr><td>Number of nodes on the circumference</td><td>150</td><td>-</td></tr></tbody></table>

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 $\delta_i$ is the relative normal displacement between one pad and the disc:

$$F_{contact,disc}^{i} = \begin{cases} k_{L}\delta_{i} + k_{NL}\delta_{i}^{3} & \text{if } \delta_{i} \ge 0 \\ 0 & \text{otherwise} \end{cases} \quad (1)$$

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:

$$F_{contact,pad}^{i} = -F_{contact,disc}^{i} \quad (2)$$

The tangential reaction direction are given by considering the classical Coulomb's law with a constant friction coefficient $\mu$ and a steady sliding motion without any stick-slip event:

$$F_{friction,disc}^{i} = -F_{friction,pad}^{i} = \mu \|F_{contact,disc}^{i}\| \vec{t}_{i} \quad (3)$$

where $\vec{t}_i$, that defines the local orthoradial in-plane direction, can be defined using the polar coordinates: $\vec{t}_i = \vec{r}_i \times \vec{z}$. Hence the global nonlinear force $F_{NL}$ is:

$$F_{NL} = F_{contact} + F_{friction} \quad (4)$$

The set of equations describing the dynamic brake system can be written by:

$$M\ddot{X} + C\dot{X} + KX = F_{NL} + F_{ext} \quad (5)$$

where $X$ is the generalized displacements vector while the dot denotes derivative with respect to time. $F_{ext}$ corresponds to the

<figure><img src="image_5.png" /><figcaption>**Fig. 1.** Finite Element Model of the simplified brake system.</figcaption></figure>
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:

$$ \mathbf{M} = \begin{pmatrix} \mathbf{M}_{\text{disc}} & 0 & 0 \\ 0 & \mathbf{M}_{\text{pad1}} & 0 \\ 0 & 0 & \mathbf{M}_{\text{pad2}} \end{pmatrix} \quad \text{and} \quad \mathbf{K} = \begin{pmatrix} \mathbf{K}_{\text{disc}} & 0 & 0 \\ 0 & \mathbf{K}_{\text{pad1}} & 0 \\ 0 & 0 & \mathbf{K}_{\text{pad2}} \end{pmatrix} \qquad (6) $$

The damping matrix **C** is estimated by considering a Rayleigh's damping under the form of **C** = α**K** + β**M**. The coefficients α and β are chosen to have approximately 0.5% of damping for the two frequencies at f₁ = 2kHz and f₂ = 5kHz (i.e. α = 2ζ(ω₁ω₂/ω₁+ω₂) and β = 2ζ(ω₁+ω₂) with ζ = 0.5% and ωᵢ = 2πfᵢ for i = 1,2). As a reminder, the use of the Rayleigh's damping formulation tends to overdamped modes outside the interval of the two chosen frequencies f₁ and f₂, 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.

## 2.2. 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**<sub>eq</sub>, which is defined by **K**<sub>**X**<sub>eq</sub></sub> = **F**<sub>ext</sub> + **F**<sub>NL</sub>(**X**<sub>eq</sub>). Assuming small perturbations δ**X** around **X**<sub>eq</sub> such as **X** = **X**<sub>eq</sub> + δ**X** the following linearized equations of motion are obtained:

$$ \mathbf{M}\delta\ddot{\mathbf{X}} + \mathbf{C}\delta\dot{\mathbf{X}} + \underbrace{(\mathbf{K} - \mathbf{J}_{\text{FNL}}(\mathbf{X}_{\text{eq}}))}_{\mathbf{K}^{\text{tot}}} \delta\mathbf{X} = \mathbf{0} \qquad (7) $$

where **J**<sub>FNL</sub>(**X**<sub>eq</sub>) is the Jacobian of the nonlinear reactions **F**<sub>NL</sub> in **X**<sub>eq</sub>. 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**<sub>eq</sub>. The eigenvalue problem associated is given by:

$$ (\lambda_k^2 \mathbf{M} + \lambda_k \mathbf{C} + \mathbf{K}^{\text{tot}}) \mathbf{X}_k = \mathbf{0} \qquad (8) $$

where (λ<sub>k</sub>, **X**<sub>k</sub>) are respectively the eigenvalues and eigenvectors of the discretized problem, with complex eigenvalues λ<sub>k</sub> = a<sub>k</sub> + ib<sub>k</sub>. The imaginary part b<sub>k</sub> defines the pulsation of the mode while the real part a<sub>k</sub> is related to the stability of the system. If at least one eigenmode have a strictly positive real part (∃k, a<sub>k</sub> > 0), the system is locally unstable around the equilibrium point **X**<sub>eq</sub>.

# 3. 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.

## 3.1. 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 (Ψⁿ) and fixed interface modes (Φⁿ), 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<sub>max</sub> (in the present study f<sub>max</sub> = 20 kHz) and a coefficient γ as such as, for each substructure we keep the fixed interface modes which frequency is less than f<sub>c</sub> = γf<sub>max</sub>. The internal dof (u<sub>int</sub><sup>d</sup>, u<sub>int</sub><sup>p1</sup>, u<sub>int</sub><sup>p2</sup>) are expressed in function of the generalized dof (ξ<sup>d</sup>, ξ<sup>p1</sup>, ξ<sup>p2</sup>, corresponding to the fixed interface modes Φⁿ) and interface dof (u<sub>l</sub><sup>d</sup>, u<sub>l</sub><sup>p1</sup>, u<sub>l</sub><sup>p2</sup>):

$$ \mathbf{u}_{\text{int}}^{n} = \mathbf{\Phi}^{n} \xi^{n} + \mathbf{\Psi}^{n} \mathbf{u}_{l}^{n} \qquad (9) $$

This allows to define the following transfer matrix **T** between the complete model and the condensed Craig & Bampton model:

$$ \mathbf{X} = \begin{pmatrix} \mathbf{u}_l^d \\ \mathbf{u}_\text{int}^d \\ \mathbf{u}_l^{p1} \\ \mathbf{u}_\text{int}^{p1} \\ \mathbf{u}_l^{p2} \\ \mathbf{u}_\text{int}^{p2} \end{pmatrix} = \underbrace{\begin{pmatrix} \mathbf{I}^d & 0 & 0 & 0 & 0 & 0 \\ \mathbf{\Psi}^d & \mathbf{\Phi}^d & 0 & 0 & 0 & 0 \\ 0 & 0 & \mathbf{I}^{p1} & 0 & 0 & 0 \\ 0 & 0 & \mathbf{\Psi}^{p1} & \mathbf{\Phi}^{p1} & 0 & 0 \\ 0 & 0 & 0 & 0 & \mathbf{I}^{p2} & 0 \\ 0 & 0 & 0 & 0 & \mathbf{\Psi}^{p2} & \mathbf{\Phi}^{p2} \end{pmatrix}}_{\mathbf{T}} \begin{pmatrix} \mathbf{u}_l^d \\ \xi^d \\ \mathbf{u}_l^{p1} \\ \xi^{p1} \\ \mathbf{u}_l^{p2} \\ \xi^{p2} \end{pmatrix} = \mathbf{T}\mathbf{X}_{CB} \qquad (10) $$

Finally matrices of the equation of motion given in Eqs. (5) and (7) can be expressed in the Craig & Bampton basis as follows:

$$ \mathbf{M}_{CB} = \mathbf{T}^T \mathbf{M} \mathbf{T} \quad \text{and} \quad \mathbf{K}_{CB} = \mathbf{T}^T \mathbf{K} \mathbf{T} \quad \text{and} \quad \mathbf{C}_{CB} = \mathbf{T}^T \mathbf{C} \qquad (11) $$

$$ \mathbf{F}_{\text{NL,CB}} = \mathbf{T}^T \mathbf{F}_{\text{NL}} \quad \text{and} \quad \mathbf{F}_{\text{ext,CB}} = \mathbf{T}^T \mathbf{F}_{\text{ext}} \qquad (12) $$

$$ \mathbf{K}_{CB}^{\text{tot}} = \mathbf{T}^T (\mathbf{K} - \mathbf{J}_{\text{FNL}}(\mathbf{X}_{\text{eq}}))\mathbf{T} \qquad (13) $$

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.

## 3.2. 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 **T**<sub>2</sub>

$$ \begin{pmatrix} \mathbf{u}_l^d \\ \mathbf{u}_l^{p1} \\ \mathbf{u}_l^{p2} \end{pmatrix} = \underbrace{\begin{pmatrix} \mathbf{I}^d & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & \mathbf{I}^{p1} & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & \mathbf{I}^{p2} & 0 \end{pmatrix}}_{\mathbf{T}_2^T} \begin{pmatrix} \mathbf{u}_l^d \\ \xi^d \\ \mathbf{u}_l^{p1} \\ \xi^{p1} \\ \mathbf{u}_l^{p2} \\ \xi^{p2} \end{pmatrix} = \mathbf{T}_2^T \mathbf{X}_{CB} \qquad (14) $$

it is possible to define the following reduced mass and stiffness matrices M̃ and K̃ given by
$$ \tilde{\mathbf{M}} = \mathbf{T}_2^t \mathbf{M}_{CB} \mathbf{T}_2 \quad \text{and} \quad \tilde{\mathbf{K}} = \mathbf{T}_2^t \mathbf{K}_{CB}^{\text{tot}} \mathbf{T}_2 \qquad (15) $$

As $\mathbf{K}_{CB}^{\text{tot}}$ depends on $\mu$, a DMS reduction for each equilibrium point and so for each set of parameters chosen has to be performed unlike the Craig & Bampton reduction which is valid for any set of parameters. The $i^{th}$ right and left eigenmodes $\mathbf{X}_{Br}^i$ and $\mathbf{X}_{Bl}^i$ are then computed, sorted by order of increasing frequencies, and we define the matrices $\varphi_r$ and $\varphi_l$ which contain the first $N_{DMS}$ eigenmodes:

$$ \left\{ \begin{array}{l} (-\omega_l^2 \tilde{\mathbf{M}} + \tilde{\mathbf{K}}) \mathbf{X}_{Br}^i = 0 \quad \text{with} \quad \varphi_r = [\mathbf{X}_{Br}^i \cdots \mathbf{X}_{Br}^{N_{DMS}}] \\ (-\omega_l^2 \tilde{\mathbf{M}}^t + \tilde{\mathbf{K}}^t) \mathbf{X}_{Bl}^i = 0 \quad \text{with} \quad \varphi_l = [\mathbf{X}_{Bl}^i \cdots \mathbf{X}_{Bl}^{N_{DMS}}] \end{array} \right. \qquad (16) $$

The size of this step of reduction is set by the number of interfaces modes kept $N_{DMS}$, associated with the degrees of freedom $\eta_l$. Using the right and left matrices of transfer $\mathbf{T}_r$ and $\mathbf{T}_l$ defined by

$$ \mathbf{T}_r = \begin{pmatrix} \varphi_r & 0 \\ 0 & \mathbf{I}^{d+p_1+p_2} \end{pmatrix} \quad \text{and} \quad \mathbf{T}_l = \begin{pmatrix} \varphi_l & 0 \\ 0 & \mathbf{I}^{d+p_1+p_2} \end{pmatrix} \qquad (17) $$

and assuming a rearrangement of the dof

$$ \mathbf{X}_{\text{DMS}} = \begin{pmatrix} \eta_l \\ \xi^{d} \\ \xi^{p_1} \\ \xi^{p_2} \end{pmatrix} = \mathbf{T}_r^{-1} \mathbf{X}_{\text{CB}} = \mathbf{T}_r^{-1} \begin{pmatrix} \mathbf{u}_l^d \\ \mathbf{u}_l^{p_1} \\ \mathbf{u}_l^{p_2} \\ \xi^d \\ \xi^{p_1} \\ \xi^{p_2} \end{pmatrix} \qquad (18) $$

the reduced matrices for each contribution of the FEM brake system are given by

$$ \mathbf{M}_{\text{DMS}} = \mathbf{T}_r^t \mathbf{M}_{CB} \mathbf{T}_r \quad \text{and} \quad \mathbf{C}_{\text{DMS}} = \mathbf{T}_r^t \mathbf{C}_{CB} \mathbf{T}_r \qquad (19) $$

$$ \mathbf{F}_{\text{NL.DMS}} = \mathbf{T}_r^t \mathbf{F}_{\text{NL.CB}} \quad \text{and} \quad \mathbf{F}_{\text{ext.DMS}} = \mathbf{T}_r^t \mathbf{F}_{\text{ext.CB}} \qquad (20) $$

$$ \mathbf{K}_{\text{DMS}}^{\text{tot}} = \mathbf{T}_r^t \mathbf{K}_{CB}^{\text{tot}} \mathbf{T}_r \qquad (21) $$

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.

# 4. 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.

## 4.1. 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 $\mu$ for 501 samples equally distributed between [0; 1] is chosen. Then eigenvalues of the system for each value of $\mu$ are computed. Within the range of study for $\mu = [0; 1]$ there are two unstable modes: the first one is detected at $\mu_1 = 0.67$ and the second one at $\mu_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.

## 4.2. 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 $\gamma = 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 $\gamma = 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 $\delta r$ and error on the imaginary parts $\delta \omega$ are defined by:

$$ \delta r = \frac{1}{N^{\text{tot}}} \sum_{k=1}^{N^{\text{tot}}} \left| \frac{a_k^{\text{ref}} - a_k}{a_k^{\text{ref}}} \right| \qquad (22) $$

$$ \delta \omega = \frac{1}{N^{\text{tot}}} \sum_{k=1}^{N^{\text{tot}}} \left| \frac{b_k^{\text{ref}} - b_k}{b_k^{\text{ref}}} \right| \qquad (23) $$

where $N^{\text{tot}}$ denotes the number of eigenvalues in the frequency range of interest (i.e. [0; 20000] Hz for the present study). $a_k^{\text{ref}}$ ($b_k^{\text{ref}}$, 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 $\delta r^{\text{un}}$ and $\delta \omega^{\text{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

$$ \delta r^{\text{un}} = \frac{1}{N^{\text{un}}} \sum_{k=1}^{N^{\text{un}}} \left| \frac{a_k^{\text{ref}} - a_k}{a_k^{\text{ref}}} \right| \qquad (24) $$

$$ \delta \omega^{\text{un}} = \frac{1}{N^{\text{un}}} \sum_{k=1}^{N^{\text{un}}} \left| \frac{b_k^{\text{ref}} - b_k}{b_k^{\text{ref}}} \right| \qquad (25) $$

where $N^{\text{un}}$ defines the number of unstable modes. In the present study, $N^{\text{un}}$ is equal to 0, 1 or 2 depending on the value of the coefficient of friction. In this case, $a_k^{\text{ref}}$ ($b_k^{\text{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 $\mathbf{M}^{\text{Mac}}$ (see [18–20]) is used to compare the eigenmodes shapes of the reference $\mathbf{X}_i^{\text{ref}}$ and those issued from a condensed system $\mathbf{X}_j$:

$$ \mathbf{M}_{ij}^{\text{Mac}} = \frac{\left| \mathbf{X}_i^{\text{ref}} \mathbf{t} \mathbf{X}_j \right|^2}{\left| \mathbf{X}_i^{\text{ref}} \right| \left| \mathbf{X}_j \right|} \qquad (26) $$

This matrix $\mathbf{M}^{\text{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
<figure><img src="image_6.png" /><figcaption>Fig. 2. Evolution of (a) real parts of eigenvalues and (b) coalescence pattern for the unstable modes 1 and 2 ($f_1$ = 857 Hz and $f_2$ = 2223 Hz).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 3. Unstable modes for $\mu$ = 0.8 (a) 1st mode at 857 Hz and (b) 2nd mode at 2223 Hz.</figcaption></figure>

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.

## 4.3. Results

Considering the previous results given in Section 4.1, four different cases are undertaken in order to validate the efficiency the DMS reduction:

* $\mu$ = 0 : the case with no friction;

* $\mu$ = 0.7 : a case with one unstable mode predicted by the stability analysis;

* $\mu$ = 0.8 : a case with two unstable modes with real parts very close;

* $\mu$ = 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 $\delta r$ and $\delta \omega$ 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 ($\mu$ = 0, $\mu$ = 0.7, $\mu$ = 0.8 and $\mu$ = 1) a number of DMS modes above 150 is sufficient to allow a lower error on the real parts $\delta r$ (on the imaginary parts $\delta \omega$, 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. $\delta r^{un}$ and $\delta \omega^{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.

Then, Figs. 6 show the adapted MAC-criterion defined by the matrix $M^{Mac}$ in the case of $\mu$ = 0.8 for the reduction basis composed of 50 and 150 interface modes. This statistic indicator illustrates the degree of consistency between mode shapes of the reference model and those of the DMS reduction basis thereby indicating the limit of validity of the reduction basis depending on the size of the DMS basis. Along the diagonal, every mode of the DMS basis is compared to the identical mode coming from the reference model. Not surprisingly the number of values larger than 0.9 (that indicate consistent correspondence whereas small values indicate poor resemblance of the mode shapes) increases along the diagonal of the matrix $M^{Mac}$ with the size of the DMS basis. Off of the diagonal, the MAC values are generally very low. If the value is close to one this implies that the modes shapes have similarities and are practically identical even if each mode should be uniquely observed and have a different shape than the other modes.
<figure><img src="image_8.png" /><figcaption>Fig. 4. Error criterion (a) $\delta r$ and (b) $\delta\omega$ (for $\mu = 0, \mu = 0.7, \mu = 0.8$ and $\mu = 1$ from top to bottom).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 5. Error criterion (a) $\delta r^{un}$ and (b) $\delta\omega^{un}$ (for $\mu = 0.7, \mu = 0.8$ and $\mu = 1$ from top to bottom).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 6. Comparison of MAC matrices for $\mu = 0.8$ by considering (a) 50 interface modes and (b) 150 interface modes.</figcaption></figure>
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 $\mu = 0$) has been previously proposed and discussed in [10].

## 5. 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 $\theta$-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 self-excited 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 $\mu = 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 $\mu = 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).

### 5.1. Time series - transient dynamics and self-excited vibrations

First of all the time series and the transient nonlinear behavior of the system in the cases of $\mu = 0.8$ and $\mu = 1$ are illustrated. Figs. 8(a-b) and 9(a-b) show the transient time evolutions for both Point1 and Point2 in the X, Y and Z directions for $\mu = 0.8$ and $\mu = 1$.

<figure><img src="image_11.png" /><figcaption>**Fig. 7.** Position of Point1 (red) and Point2 (blue). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

by considering the reference model. Considering the case of $\mu = 0.8$ (see Figs. 8(a) and (b)), the transient vibrations of Point1 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)). Finally the last phase(for $t > 0.7$ s) is characterized by a global saturation of the all the nonlinear oscillations which therefore reflects the establishment of a stationary regime of self-sustained vibrations. For the case $\mu = 1$, a succession of three main phases is also visible: the first one for $t = [0; 0.1]$ s, the second one for $t = [0.1; 0.4]$ s and the last one for $t > 0.4$ s. The main difference compared to the previous case $\mu = 0.8$ lies in the final establishment of a stabilized regime with regular waves of oscillations, which characterizes a quasi-periodic behaviour of self-sustained vibrations.

In order to evaluate the performance of the reduction on accurately reproducing these transient behaviours, it was chosen to compare the reference results to four choices of DMS sizes with 50, 150, 500 and 1000 interface modes, respectively. Results are provided in Figs. 8(c-j) and 9(c-j) for $\mu = 0.8$ and $\mu = 1$, respectively. By looking at the results in an overall way, it is clearly observed that the evolution of the transient vibrations for the reference model and the different reduced systems are quite similar. The three successive phases are reproduced with the final occurrence of a stationary regime of self-sustained vibrations for $\mu = 0.8$ and a quasi-periodic regime for $\mu = 1$.

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 $\mu = 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 $\mu = 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.
<figure><img src="image_12.png" /><figcaption>**Fig. 8.** Transient signals for the case $\mu = 0.8$ (a) Point1 - Reference (b) Point2 - Reference (c) Point1 - $N_{DMS}=50$ (d) Point2 - $N_{DMS}=50$ (e) Point1 - $N_{DMS}=150$ (f) Point2 - $N_{DMS}=150$ (g) Point1 - $N_{DMS}=500$ (h) Point2 - $N_{DMS}=500$ (i) Point1 - $N_{DMS}=1000$ (j) Point2 - $N_{DMS}=1000$.</figcaption></figure>

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
<figure><img src="image_13.png" /><figcaption>**Fig. 9.** Transient signals for the case $\mu = 1$ (a) Point1 - Reference (b) Point2 - Reference (c) Point1 - $N_{DMS}=50$ (d) Point2 - $N_{DMS}=50$ (e) Point1 - $N_{DMS}=150$ (f) Point2 - $N_{DMS}=150$ (g) Point1 - $N_{DMS}=500$ (h) Point2 - $N_{DMS}=500$ (i) Point1 - $N_{DMS}=1000$ (j) Point2 - $N_{DMS}=1000$.</figcaption></figure>
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.

## 5.2. 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 $\mu = 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 \times$ 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 $\mu = 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 \times$ 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 $\mu = 0.8$. In the same way the limit cycles and FFT are given in Figs. 14 and 15 for $\mu = 1$.

For the case $\mu = 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 $\mu = 1$ results based on DMS reductions are also in agreement with the reference for the prediction of the quasi-periodic 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 quasi-periodic established regime with frequency contents very similar to those of the reference. The main frequency contributions of

<figure><img src="image_14.png" /><figcaption>**Fig. 10.** Spectrograms of the signals for the case $\mu = 0.8$ (a) Point1 - Reference (b) Point2 - Reference.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Fig. 11. Spectrograms of the signals for the case $\mu = 1$ (a) Point1 - Reference (b) Point2 - Reference.</figcaption></figure>

<table><caption>**Table 2**

Frequency of the two unstable modes for the reference and the reduced models.</caption><thead><tr><th></th><th></th><th>CEA</th><th>Reference</th><th>N<sub>DMS</sub> = 50</th><th>N<sub>DMS</sub> = 150</th><th>N<sub>DMS</sub> = 500</th><th>N<sub>DMS</sub> = 1000</th></tr></thead><tbody><tr><td rowspan="2">μ = 0.8</td><td>1st mode (Hz)</td><td>858</td><td>857</td><td>857</td><td>857</td><td>857</td><td>857</td></tr><tr><td>2nd mode (Hz)</td><td>2223</td><td>2205</td><td>2202</td><td>2203</td><td>2205</td><td>2205</td></tr><tr><td rowspan="2">μ = 1</td><td>1st mode (Hz)</td><td>857</td><td>857</td><td>856</td><td>857</td><td>857</td><td>857</td></tr><tr><td>2nd mode (Hz)</td><td>2224</td><td>2195</td><td>2195</td><td>2195</td><td>2195</td><td>2195</td></tr></tbody></table>

<figure><img src="image_16.png" /><figcaption>Fig. 12. Limit cycles on the Point1 and associated Fast Fourier Transform for $\mu = 0.8$, (black) Reference model, (red) DMS = 50, (green) DMS = 150 and (blue) DMS = 500. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

order 4 or less ($|m| + |n| \le 4$) spotted in the Fast Fourier Transforms of the limit cycles for the latter case $\mu = 1$ are depicted in 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.
<figure><img src="image_17.png" /><figcaption>**Fig. 13.** Limit cycles on the *Point*2 and associated Fast Fourier Transform for $\mu = 0.8$, (black) Reference model, (red) DMS = 50, (green) DMS = 150 and (blue) DMS = 500. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Fig. 14.** Limit cycles on the *Point*1 and associated Fast Fourier Transform for $\mu = 1$, (black) Reference model and (red) DMS = 150. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>**Fig. 15.** Limit cycles on the Point2 and associated Fast Fourier Transform for $\mu = 1$, (black) Reference model and (red) DMS = 150. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

### 5.3. 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 $\mathbf{Y}(t)$ in state-

space variables (i.e. $\mathbf{Y} = (\dot{\mathbf{X}} \ \mathbf{X})^t$), the evolution of the modal participation $\beta_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:

$$ \beta_i(t) = \frac{\overline{\mathbf{Y}}_i^{\text{ref.ins}} t \tilde{\mathbf{A}} \mathbf{Y}(t)}{\overline{\mathbf{Y}}_i^{\text{ref.ins}} t \tilde{\mathbf{A}} \mathbf{Y}_i^{\text{ref.ins}}} \quad \text{with} \quad \tilde{\mathbf{A}} = \begin{pmatrix} \mathbf{M} & \mathbf{C}/2 \\ \mathbf{C}/2 & \mathbf{K} \end{pmatrix} \qquad (27) $$

where $\mathbf{Y}_i^{\text{ref.ins}}$ define the modes shapes in state space. The matrix $\tilde{\mathbf{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 $\langle \beta_i(t) \rangle$ 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 $\mu = 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><caption>**Table 3**

Frequencies and corresponding combinations for the case $\mu = 1$ on the limit cycles.</caption><thead><tr><th>Frequency (Hz)</th><th>Combinations</th></tr></thead><tbody><tr><td>857</td><td>f<sub>1</sub></td></tr><tr><td>1366</td><td>-f<sub>1</sub> + f<sub>2</sub></td></tr><tr><td>1714</td><td>2f<sub>1</sub></td></tr><tr><td>2223</td><td>f<sub>2</sub></td></tr><tr><td>2571</td><td>3f<sub>1</sub></td></tr><tr><td>2732</td><td>-2f<sub>1</sub> + 2f<sub>2</sub></td></tr><tr><td>3080</td><td>f<sub>1</sub> + f<sub>2</sub></td></tr><tr><td>3428</td><td>4f<sub>1</sub></td></tr><tr><td>3589</td><td>-f<sub>1</sub> + 2f<sub>2</sub></td></tr><tr><td>3937</td><td>2f<sub>1</sub> + f<sub>2</sub></td></tr><tr><td>4446</td><td>2f<sub>2</sub></td></tr><tr><td>4794</td><td>3f<sub>1</sub> + f<sub>2</sub></td></tr><tr><td>6669</td><td>3f<sub>2</sub></td></tr></tbody></table>
<figure><img src="image_20.png" /><figcaption>Fig. 16. Evolution of the modal participations for the case $\mu = 0.8$ (dashed lines = mode 1, lines = mode 2; black = reference, red = DMS50, green = DMS150, blue = DMS500, magenta = DMS1000). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Fig. 17. Evolution of the modal participations for the case $\mu = 1$ (dashed lines = mode 1, lines = mode 2; black = reference, red = DMS50, green = DMS150, blue = DMS500, magenta = DMS1000). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

* $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 $\mu = 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 non-linear 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: $\mu = 0.8$ and $\mu = 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.

## 6. Acoustic radiation

The prediction of acoustic radiated noise from friction-induced vibration is a recent topic of interest [14,15,25-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.

### 6.1. 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 $\mathcal{S}$, which consists in the resolution of the Kirchhoff-Helmholtz surface integral, giving the field pressure $\mathbf{P}$ at a point of coordinate $\mathbf{r}$, such as:

$$ \epsilon(\mathbf{r})\mathbf{P}(\mathbf{r}) = \int_{\mathcal{S}} \left( \mathbf{P}(\mathbf{s}) \frac{\partial G(\mathbf{r}, \mathbf{s})}{\partial \mathbf{n}} - G(\mathbf{r}, \mathbf{s}) \frac{\partial \mathbf{P}(\mathbf{s})}{\partial \mathbf{n}} \right) d\mathcal{S} \quad (28) $$

with $\epsilon$ a parameter equal to 1 when strictly outside of the surface and taking a value in $[0, 1]$, when $\mathbf{r} \in \mathcal{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 $\mathbf{P}_s$ has to be determined. The pressure field $\mathbf{P}$ at any point outside of the system can then be computed. The surface of the system $\mathcal{S}$ is discretized in $\mathcal{S} = \cup_{i} \mathcal{S}_i$ 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 $\dot{\mathbf{X}}_n$ at each point of the mesh, and then selecting manually the $K$ main frequency contributions
<table><caption>**Table 4**
Values of the stabilized modal participations for both cases $\mu = 0.8$ and $\mu = 1$ and the different sizes of DMS reduction</caption><thead><tr><th></th><th></th><th>Reference</th><th>N<sub>DMS</sub> = 50</th><th>N<sub>DMS</sub> = 150</th><th>N<sub>DMS</sub> = 500</th><th>N<sub>DMS</sub> = 1000</th></tr></thead><tbody><tr><td rowspan="2">μ = 0.8</td><td>1st mode</td><td>19.9</td><td>21.1</td><td>20.3</td><td>20.1</td><td>19.9</td></tr><tr><td>2nd mode</td><td>2.1</td><td>2.3</td><td>2.2</td><td>2.1</td><td>2.1</td></tr><tr><td rowspan="2">μ = 1</td><td>1st mode</td><td>16.0</td><td>15.9</td><td>16.0</td><td>16.0</td><td>16.0</td></tr><tr><td>2nd mode</td><td>3.9</td><td>3.9</td><td>3.9</td><td>3.9</td><td>3.9</td></tr></tbody></table>

$\mathbf{f}_{BEM} = \{f_{BEM}^1, \dots, f_{BEM}^K\}$. Those frequencies may be the fundamental frequencies $f_1$ or $f_2$, their harmonics $nf_1$ or $mf_2$ (with $n > 1$ and $m > 1$ and integer), or linear combinations $mf_1 + nf_2$ (with $(m, n) \in \mathbb{Z}^2$). The velocity field $\dot{\mathbf{X}}_n$ can then be estimated as follows:

$$ 
\begin{aligned}
\dot{\mathbf{X}}_n(t) &\approx \sum_{i=1}^{K} \dot{\mathbf{X}}_n^i(t) \\
&= \sum_{i=1}^{K} \mathbf{a}_{BEM}^i \cos(\omega^i t) + \mathbf{b}_{BEM}^i \sin(\omega^i t) \quad \text{with } \omega^i \\
&= 2\pi f_{BEM}^i
\end{aligned} 
\quad (29) $$

where $\mathbf{a}_{BEM}^i$ and $\mathbf{b}_{BEM}^i$ define the Fourier vectors at each frequency $f_{BEM}^i$. Due to the fact that the numerical resources necessary to compute the external pressure responses for all the frequencies selected in $\mathbf{f}_{BEM}$ can rapidly become prohibitive, a selection on the optimal sample of these frequencies (denoted $\mathbf{f}_{BEM, opt}$) may be used to reduce the number of calculations to be performed. This method has been successfully tested in [14]. It is based on the convergence of the mean square velocity of the surface velocity field $\dot{\mathbf{X}}_n$ with the $K_{opt}$ different major components of $\mathbf{f}_{BEM}$. For more details the interested reader is referred to [14].

For each component $\dot{\mathbf{X}}_n^i$ (and so a specific frequency), the BEM method first computes the associated surface pressure $\mathbf{P}_s^i$ assuming that the velocities of the solid and the fluid are equal at the interface. Then, the associated pressure field $\mathbf{P}^i$ on any exterior set of

observation points can then be computed. Finally, squeal noise that corresponds to the acoustic field $\mathbf{P}_{ext}$ is defined by

$$ \mathbf{P}_{ext} = \sum_{i=1}^{K_{opt}} \mathbf{P}^i \quad (30) $$

where $K_{opt}$ is the optimal number of selected frequencies taken from the frequency contributions $\mathbf{f}_{BEM}$. Of course the pressure level in dB can be calculated such as

$$ \mathbf{L}_{BEM}^p = 10\log_{10}\left(\frac{\mathbf{P}_{ext}^*\mathbf{P}_{ext}}{P_{ref}^2}\right) \quad (31) $$

where $P_{ref}$ corresponds to the classical reference value used. (i.e. $P_{ref} = 2.10^{-5}$ Pa).

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.

<figure><img src="image_22.png" /><figcaption>**Fig. 18.** Acoustic patterns on the horizontal plane at a distance $Z_{square} = 0.75$ m for different Craig & Bampton reductions, scale in dB.</figcaption></figure>
## 6.2. Acoustic results

Figs. 18–21 illustrate the acoustic results on the horizontal plane square placed at a distance $Z_{square} = 0.75\,\text{m}$ and the sphere of radius $R_{sphere} = 0.5\,\text{m}$ for both different Craig & Bampton and DMS reductions. The convergence of the Craig & Bampton reduction is tested by considering five size reduction (i.e. $\gamma = [1; 1.25; 1.5; 2.5; 3]$). As a reminder, the performance of the DMS reduction is performed by considering the system with the Craig & Bampton reduction with $\gamma = 1.5$ as the reference. Four sizes

of DMS reductions are investigated (i.e. $N_{DMS} = [50; 150; 500; 100]$). Moreover, the squeal noise is calculated for $\mu = 0.8$ and $\mu = 1$ during the stationary regime as well as the transient vibratory regime for $\mu = 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$;

<figure><img src="image_23.png" /><figcaption>Fig. 19. Acoustic patterns on the horizontal plane at a distance $Z_{square} = 0.75\,\text{m}$ for different DMS reductions, scale in dB.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>Fig. 20. Acoustic patterns on the sphere of radius $R_{sphere} = 0.5\,\text{m}$ for different Craig & Bampton reductions, scale in dB.</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>Fig. 21. Acoustic patterns on the sphere of radius $R_{sphere}$ = 0.5 m for different DMS reductions, scale in dB.</figcaption></figure>

<table><caption>**Table 5**

Time intervals (in second as unit of time) of the transient signal used to compute the FFT for the acoustic calculations.</caption><thead><tr><th></th><th>CB γ = 1.5</th><th>N<sub>DMS</sub> = 50</th><th>N<sub>DMS</sub> = 150</th><th>N<sub>DMS</sub> = 500</th><th>N<sub>DMS</sub> = 1000</th></tr></thead><tbody><tr><td>Phase I, μ = 0.8</td><td>[0.37:0.39]</td><td>[0.16:0.18]</td><td>[0.21:0.23]</td><td>[0.26:0.28]</td><td>[0.35:0.37]</td></tr><tr><td>Phase II, μ = 0.8</td><td>[0.55:0.65]</td><td>[0.3:0.4]</td><td>[0.33:0.43]</td><td>[0.45:0.65]</td><td>[0.5:0.6]</td></tr><tr><td>Phase III, μ = 0.8</td><td>[1.3:1.4]</td><td>[1.3:1.4]</td><td>[1.3:1.4]</td><td>[1.3:1.4]</td><td>[1.3:1.4]</td></tr><tr><td>Phase III, μ = 1</td><td>[1.3:1.4]</td><td>[1.3:1.4]</td><td>[1.3:1.4]</td><td>[1.3:1.4]</td><td>[1.3:1.4]</td></tr></tbody></table>

* Phase II: the transition part of the signal between the first pseudo-established vibrational regime and the final vibrational regime;

* Phase III: the final vibrational regime.

For the reader comprehension, the various time intervals and the value of the friction coefficient are indicated on the left side of Figs. 18–21. The chosen reduction method is also indicated at the top of the figures. The first column corresponds to the reference solution.

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 μ = 0.8 and μ = 1), but also the transient vibrational regime (i.e. Phases I and II for μ = 0.8);

* from a general point of view, the acoustic patterns predicted by the different DMS reductions (for a specific phase and μ) 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).

## 7. 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-
diction of all the transient squeal phenomena or the estimation of squeal noise.

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.

## Acknowledgments

J.-J. Sinou acknowledges the support of the Institut Universitaire de France.

## References

[1] Ibrahim RA. Friction-induced vibration, chatter, squeal and chaos part 1: mechanics of contact and friction. Appl Mech Rev 1994;47:209-26.

[2] Ibrahim RA. Friction-induced vibration, chatter, squeal and chaos part 2: Dynamics and modelling. Appl Mech Rev 1994;47:227-53.

[3] Kindkaid NM, O'Reilly OM, Papadopoulos P. Automatic disc brake squeal. J Sound Vib 2003;267:105-66.

[4] Ouyang H, Nack W, Yuan Y, Chen F. Numerical analysis of automatic disc brake squeal: a review. Int J Veh Noise Vib 2005;1:207-31.

[5] Sinou J-J. Transient non-linear dynamic analysis of automotive disc brake squeal - on the need to consider both stability and non-linear analysis. Mech Res Commun 2010;37:96-105.

[6] Chevillot F, Sinou J-J, Hardouin N. Nonlinear transient vibrations and coexistences of multi-instabilities by friction in an aircraft braking system. J Sound Vib 2009;328:555-74.

[7] Vermot des Roches G. Frequency and time simulation of squeal instabilities. Application to the design of industrial automotive brakes PhD thesis. École centrale Paris; 2011.

[8] Brizard D, Chiello O, Sinou J-J, Lorang X. Performances of some reduced bases for the stability analysis of a disc/pads system in sliding contact. J Sound Vib 2011;330:703-20.

[9] Loyer A, Sinou J-J, Chiello O, Lorang X. Study of nonlinear behaviors and modal reductions for friction destabilized systems. application to an elastic layer. J Sound Vib 2012;331:1011-41.

[10] Monteil M, Besset S, Sinou J-J. A double modal synthesis approach for brake squeal prediction. Mech Syst Signal Process 2016;70-71:1073-84.

[11] Lai V-V, Chiello O, Brunel J-F, Dufrenoy P. Study of nonlinear behaviors and modal reductions for friction destabilized systems. application to an elastic layer. J Sound Vib 2019;331:197-215.

[12] Besset S, Sinou J-J. Modal reduction of brake squeal systems using complex interface modes. Mech Syst Signal Process 2017;85:896-911.

[13] Sinou J-J, Besset S. Simulation of transient nonlinear friction-induced vibrations using complex interfaces modes: application to the prediction of squeal events. Shock Vib 2017.

[14] Soobbarayen K, Besset S, Sinou J-J. A simplified approach for the calculation of acoustic emission in the case of friction-induced noise and vibration. Mech Syst Signal Process 2015;50-51:732-56.

[15] Soobbarayen K, Besset S, Sinou J-J. Noise and vibration for a self-excited mechanical system with friction. Appl Acoust 2013;74:1191-204.

[16] Soobbarayen K, Sinou J-J, Besset S. Numerical study of friction-induced instability and acoustic radiation - effect of ramp loading on the squeal propensity for a simplified brake model. J Sound Vib 2014;333:5475-93.

[17] Craig R, Bampton M. Coupling of substructures for dynamic analyses. Am Inst Aeronaut Astronaut 1968;6:1313-9.

[18] Allemang R. The modal assurance criterion - twenty years of use and abuse. Sound and Vibration, Aug:14-21; 2003.

[19] Allemang R, Phillips A. Techniques for evaluation of modal vector contamination. Sound Vib 2015(Jan):8-11.

[20] Brizard D, Chiello O, Sinou J-J, Lorang X. Performances of some reduced bases for the stability analysis of a disc/pads system in sliding contact. J Sound Vib 2011;330:703-20.

[21] Jean M. The non-smooth contact dynamics method. Comput Methods Appl Mech Eng 1999;177:235-57.

[22] Loyer A, Sinou J-J, Chiello O, Lorang X. Study of nonlinear behaviors and modal reductions for friction destabilized systems, applications to an elastic layer. J Sound Vib 2012;47:1011-41.

[23] Lorang X, Chiello O. Stability and transient analysis in the modeling or railway disc brake squeal. Notes Numer Fluid Mech Multidiscip Des 2008;99:447-53.

[24] Sinou J-J, Loyer A, Chiello O, Mogenier G, Lorang X, Cocheteux F, Bellaj S. A global strategy based on experiments and simulations for squeal prediction on industrial railway brakes. J Sound Vib 2013;332(20):5068-85.

[25] Oberst S, Lai JCS, Marburg S. Guidelines for numerical vibration and acoustic analysis of disc brake squeal using simple models of brake systems. J Sound Vib 2013;332:2284-99.

[26] Oberst S, Lai JCS. Squeal noise in simple numerical brake models. J Sound Vib 2015;352:129-41.

[27] Lee H, Singh R. Determination of sound radiation form a simplified disk-brake rotor by a semi-analytical method. Noise Control Eng J 2004;52:225-39.

[28] Akay A. Acoustics of friction. J Acoust Soc Am 2002;11:1525-48.

[29] Juhl PM. The Boudary Element Method for Sound Field Calculations PhD thesis. Technical University of Denmark; 1993.

[30] Bonnet M. Boundary integral equation methods for solids and fluids; 1999.

[31] Juhl PM, Henriquez VC. Open bem. http://www.openbem.dk; 2015. [Online].