# Calculs de cycles limites de systèmes autonomes multi-instables par une méthode fréquentielle : Application au crissement des freins à disque

Nicolas Coudeyras, Jean-Jacques Sinou, Samuel Nacivet

► To cite this version:

Nicolas Coudeyras, Jean-Jacques Sinou, Samuel Nacivet. Calculs de cycles limites de systèmes autonomes multi-instables par une méthode fréquentielle : Application au crissement des freins à disque. 9e Colloque national en calcul des structures, CSMA, May 2009, Giens, France. hal-01413843

HAL Id: hal-01413843

https://hal.science/hal-01413843v1

Submitted on 11 Dec 2016

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>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Calculs de cycles limites de systèmes autonomes multi-instables par une méthode fréquentielle : Application au crissement des freins à disque

N. Coudeyras¹,², J-J. Sinou² and S. Nacivet¹

¹PSA Peugeot Citroën
18 rue des Fauvelles, 92256 La Garenne Colombes Cedex
{nicolas.coudeyras,samuel.nacivet}@mpsa.com

² Ecole Centrale Lyon
Laboratoire de Tribologie et Dynamique des Systèmes UMR-CNRS 5513
36 avenue Guy de Collongue, 69134 Ecully Cedex
jean-jacques.sinou@ec-lyon.fr

**Résumé** — Ce papier présente une méthode fréquentielle permettant de calculer l'état dynamique stationnaire d'un système autonome soumis à plusieurs modes instables et générant une réponse pseudo-périodique. La méthode est dérivée des méthodes de balance harmonique couramment utilisées pour des systèmes dynamiques soumis à des excitations périodiques. La méthode est appliquée sur un modèle réduit de frein automobile et les résultats sont comparés à une méthode temporelle.

**Mots clés** — Crissement de frein, multi-instabilité, non-linéaire, cycle limite.

## 1 Introduction

Le crissement des freins à disque est un problème majeur pour l'industrie automobile car il dégrade les prestations acoustiques et défavorise la qualité perçue du véhicule. Le mécanisme du crissement est issu d'un couplage de modes du système de frein dû au frottement entre les plaquettes et le disque. De nombreuses études ont été publiées sur le sujet et [1, 2] en font une revue globale. L'analyse de stabilité du système de frein qui consiste en un calcul de valeurs propres complexes est la méthode numérique classiquement employée dans la prédiction du crissement. Cette méthode a néanmoins tendance à surestimer le nombre de modes instables du système et une analyse transitoire permet d'analyser plus finement le modèle de frein car elle fournit les amplitudes de vibrations (cycles limites) et indique les modes instables réellement impliqués dans la dynamique [3]. Néanmoins cette démarche est prohibitive en temps CPU et on propose dans ce papier une démarche dérivée des méthodes de balance harmonique (HBM) qui permettent de calculer directement l'état dynamique stationnaire de la réponse et qui présentent des performances intéressantes [5]. Dans une première partie un modèle numérique du système de frein est présenté et sa stabilité analysée. Le développement de la méthode harmonique est présenté dans la deuxième partie et les résultats de calculs sont présentés et analysés puis comparés avec une méthode d'intégration numérique dans la dernière partie.
# 2 Analyse linéaire

## 2.1 Modèle numérique

Le modèle étudié considère une plaquette de frein et le disque en acier (Figure 1), ces éléments étant suffisants pour obtenir des instabilités. Une procédure de Craig & Bampton est employée sur les modèles éléments finis de ces deux sous-systèmes afin de ne conserver que quelques noeuds à l'interface de contact et une base de modes dynamiques. Ces modèles réduits sont représentatifs jusqu'à la fréquence de 10 kHz. L'équation de mouvement s'écrit classiquement :

$$ M\ddot{U} + D\dot{U} + KU + F_{nl}(U) = F_{ext} \quad (1) $$

avec respectivement M, D, K les matrices de masses, d'amortissements et de rigidités, Fnl (U) le vecteur des efforts non linéaires de contact et de frottement et Fext le vecteur des efforts extérieurs. Les vecteurs de l'accélération, de la vitesse et du déplacement s'écrivent respectivement U, U, U. Une méthode de pénalité permet de représenter le contact à l'interface disque/plaquette. Afin de prendre en compte le caractère non linéaire de déformation de la plaquette sous charge, nous choisissons des raideurs linéaires et non linéaires, c'est à dire dépendantes de l'interprétation entre les deux corps en contact. De plus, un contact unilatéral est considéré, autorisant localement le décollement de la plaquette et du disque. La forme mathématique du contact est la suivante :

$$ F_{\text{contact},i} = \begin{cases} k_l(u_1 - u_2) + k_{nl}(u_1 - u_2)^3 & \text{si } (u_1 - u_2) > 0 \\ 0 & \text{sinon} \end{cases} \quad (2) $$

avec u1 et u2 les déplacements respectifs de la plaquette et du disque au noeud de contact i, kl et knl les raideurs linéaires et non linéaires. Une loi de Coulomb est employée pour définir les efforts de frottement :

$$ F_{f,i} = \mu F_{\text{contact},i} \operatorname{sgn}(v_{r,i}) \quad (3) $$

$\mu$ est le coefficient de frottement, constant au cours du freinage et $v_{r,i}$ la vitesse relative entre les deux corps au noeud de contact i. On considère un état toujours glissant sans phénomène de stick-slip.

<figure><img src="image_2.png" /><figcaption>Figure 1 – Modèle en vue éclatée du disque de frein et d'une plaquette</figcaption></figure>

## 2.2 Équilibre statique et analyse de stabilité

La phase de freinage consiste à appliquer un effort sur la plaquette afin qu'elle rentre en contact avec le disque. Il s'agit d'une phase où l'inertie et l'amortissement des différents constituants n'ont pas d'influence sur l'équilibre des efforts du système. Ainsi l'équation d'équilibre se met sous la forme :

$$ KU_0 + F_{nl}(U_0) = F_{ext} \quad (4) $$
La position d'équilibre statique $\mathbf{U}_0$ est calculée et sa stabilité analysée grâce à une analyse aux valeurs propres complexes du système (1) linéarisé autour de sa position d'équilibre :

$$ \mathbf{M}\ddot{\mathbf{U}} + \mathbf{D}\dot{\mathbf{U}} + (\mathbf{K} + \mathbf{J}_{nl})\tilde{\mathbf{U}} = \mathbf{0} \quad (5) $$

avec $\mathbf{J}_{nl}$ la matrice Jacobienne des efforts non linéaires à la position d'équilibre statique $\mathbf{U}_0$. Les valeurs propres complexes $\lambda = a + i\omega$ sont extraites et le signe des parties réelles analysé. Si au moins une des parties réelles est positive, alors on est en présence d'un mode instable à la pulsation $\omega$ et le crissement peut apparaître. La Figure 2 présente l'évolution des parties réelles et des fréquences en fonction de $\mu$.

<figure><img src="image_3.png" /><figcaption>Figure 2 – Evolution des parties réelles (a) et des fréquences (b) (lignes pointillées : modes instables)</figcaption></figure>

On voit très clairement le phénomène de coalescence des fréquences dû au frottement (b). A $\mu = 0$ on a les modes de la structure couplée puis avec $\mu$ qui augmente, certains de ces modes se rapprochent et deviennent égaux (phénomène de verrouillage). Les parties réelles correspondantes se séparent et l'une devient positive (instable) tandis que l'autre reste négative (stable). Ici, deux modes M1 et M2 à la fréquence respective de 920Hz et 1514Hz sont instables pour $\mu > 0.28$.

# 3 Analyse non linéaire

Cette section présente la méthode de balance harmonique contrainte généralisée (GCHBM) permettant de calculer des réponses pseudo-périodiques de systèmes autonomes.

## 3.1 Fonctions quasi-périodiques

Lorsqu'un système dynamique est excité par $p$ fréquences incommensurables entre elles avec $p > 1$, la réponse engendrée est dite quasi-périodique. Si $p < +\infty$, la réponse peut se décomposer sous forme de série de Fourier :

$$ \mathbf{U}(t) \approx \sum_{k_1=-N_h}^{N_h} \dots \sum_{k_p=-N_h}^{N_h} \mathbf{a}_{\mathbf{k}_1 \dots \mathbf{k}_p} \cos(k_1 \omega_1 + \dots + k_p \omega_p) t + \mathbf{b}_{\mathbf{k}_1 \dots \mathbf{k}_p} \sin(k_1 \omega_1 + \dots + k_p \omega_p) t \quad (6) $$

$N_h$ correspond au nombre d'harmoniques retenues dans la solution, $\mathbf{a}_{\mathbf{k}_1 \dots \mathbf{k}_p}$ et $\mathbf{b}_{\mathbf{k}_1 \dots \mathbf{k}_p}$ sont les vecteurs des coefficients de Fourier de la série. Ainsi non seulement les fréquences fondamentales et leurs harmoniques sont prises en compte, mais aussi les couplages entre les $p$ fréquences. En posant $\tau = \omega t$ et en considérant $(,)$ comme le produit scalaire, l'équation (6) peut se mettre sous une forme adimensionnelle en temps :
$$ U(\tau) = a_0 + \sum_{\mathbf{k} \in \mathbb{Z}^p} a_{\mathbf{k}} \cos(\mathbf{k}, \tau) + \sum_{\mathbf{k} \in \mathbb{Z}^p} b_{\mathbf{k}} \sin(\mathbf{k}, \tau) \quad (7) $$

Dans (7), chaque fréquence évolue dans un espace temporel qui lui est propre et dont la période est égale à $2\pi$ [4]. Ainsi, on a autant de dimensions temporelles que de fréquences incommensurables dans la réponse. En injectant (7) dans (1) et en orthogonalisant le résidu par une procédure de Galerkin, l'équation du mouvement devient l'équation algébrique suivante :

$$ \Gamma \tilde{z} + \tilde{F}_{nl}(\tilde{z}) = \tilde{F}_{\text{ext}} \quad (8) $$

où $\Gamma$ représente la matrice de raideur dynamique :

$$ \Gamma = \begin{bmatrix} \mathbf{K} & \mathbf{0} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \gamma_1 & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \ddots & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \gamma_{N_c} \end{bmatrix} \quad (9) $$

avec $\gamma_i$ une matrice de raideur élémentaire, dépendante en fréquence :

$$ \gamma_i = \begin{bmatrix} -(\mathbf{k}, \omega)^2 \mathbf{M} + \mathbf{K} & (\mathbf{k}, \omega) \mathbf{D} \\ -(\mathbf{k}, \omega) \mathbf{D} & -(\mathbf{k}, \omega)^2 \mathbf{M} + \mathbf{K} \end{bmatrix} \quad \text{pour } i \in [1, N_c] \quad (10) $$

$\tilde{z}$, $\tilde{F}_{nl}$ et $\tilde{F}_{\text{ext}}$ représentent respectivement le vecteur des coefficients de Fourier des déplacements, des efforts non linéaires et des efforts externes. $N_c$ est le nombre total de composantes harmoniques. La résolution de (8) nécessite d'utiliser des méthodes itératives étant donné son caractère non linéaire. Nous utiliserons un algorithme dérivé des méthodes de Newton. Le calcul des coefficients de Fourier non linéaires étant prohibitif dans le domaine fréquentiel, nous employons une méthode AFT proposée par Cameron [7] permettant d'évaluer les efforts non linéaires dans le domaine temporel puis d'en déterminer les coefficients de Fourier par l'usage de transformées de Fourier directes et inverses, discrète ou rapide. Étant donné le caractère multi-dimensionnel de l'équation (7), on utilisera une transformée de Fourier rapide (FFT) de dimension $p = 2$ puisque nous sommes en présence de deux modes instables.

## 3.2 Equations et inconnues additionnelles

La méthode HBM telle qu'elle est décrite par l'équation (8) s'applique sur des systèmes soumis à des excitations périodiques mais ne permet pas de déterminer les niveaux vibratoires des systèmes autonomes. En effet, l'unicité de la solution n'est plus respectée lorsque ces systèmes sont instables et une réponse statique ou dynamique sont solutions du système. Ainsi il est nécessaire d'ajouter des équations supplémentaires à (8) pour ne calculer que la réponse dynamique. Cette contrainte porte sur la ou les parties réelles des modes instables en dynamique. En effet, la dynamique du système autonome est gouvernée par les modes instables et leur parties réelles donnent l'évolution des amplitudes. Une partie réelle nulle indique que le mode oscille de façon stationnaire à sa propre fréquence. Ainsi le système est à l'équilibre dynamique si les parties réelles des modes instables sont nulles. Le calcul des parties réelles du système en dynamique fait référence à la linéarisation équivalente proposée par Iwan [6] où la partie temporelle non linéaire $F_{nl}(U)$ de (1) est remplacée par une matrice Jacobienne équivalente indépendante du temps $J_{eq}$. Le système linéarisé engendré permet de procéder à un calcul de modes propres complexes où les parties réelles d'intérêts sont extraites et ajoutées à (8) comme équations supplémentaires à résoudre. De plus les fréquences de la réponse dynamique ne sont à priori pas connues et sont donc considérées comme des inconnues supplémentaires. En résumé, la méthode de balance harmonique contrainte généralisée (GCHBM) se met sous la forme suivante :
$$ \begin{cases} J_1(\tilde{z}, \Omega) = \| \Gamma(\tilde{z}, \Omega) + \tilde{F}_{nl}(\tilde{z}, \Omega) - \tilde{F}_{\text{ext}} \| < \varepsilon_1 \\ J_2(\tilde{z}, \Omega) = \| \{ \text{Re}(\lambda_j(\tilde{z}, \Omega)) \} \| < \varepsilon_2 \text{ avec } j \in [1, p] \end{cases} \quad (11) $$

On dispose d'un système bien déterminé avec autant d'inconnues que d'équations qui permet de calculer la réponse dynamique de systèmes autonomes sujets à plusieurs modes instables.

# 4 Résultats et performances

Cette section présente les résultats de calculs issus de la méthode de balance harmonique généralisée avec $N_h = 3$ puis la comparaison aux résultats d'une analyse transitoire.

## 4.1 Réponse non-linéaire comportant deux modes instables $\mu = 0.35$

La figure 3 montre le cycle limite d'un degré de liberté (ddl) du disque dans la direction normale à la surface de frottement, pour la GCHBM (a) et pour l'intégration temporelle (b).

<figure><img src="image_4.png" /><figcaption>Figure 3 – Cycles limites à $\mu = 0.35$ : (a) GCHBM, (b) Integration numérique. (c) Spectre de densité de puissance : - - GCHBM, — Integration numérique</figcaption></figure>

La comparaison des résultats entre les deux méthodes montre une bonne corrélation, les cycles limites ayant la même forme et les amplitudes extrêmes étant quasiment identiques avec une erreur inférieure à 2%. On voit très bien le caractère pseudo-périodique de la réponse puisque les cycles limites ne sont pas fermés et que la réponse occupe toute une surface. Afin de caractériser plus en détail le signal, la figure 3 (c) montre les spectres de puissance des réponses. On distingue les composantes fondamentales des modes M1 et M2 aux fréquences $f_1 = 897$Hz et $f_2 = 1497$Hz et leurs harmoniques correspondantes $2f_1, 3f_1, 2f_2$ et $3f_2$. On voit ici l'intérêt de considérer les fréquences fondamentales comme des inconnues dans la GCHBM puisqu'elles diffèrent légèrement de celles obtenues dans l'analyse linéaire. De plus on observe d'autres harmoniques qui sont issues des couplages entre les modes M1 et M2 tel que $2f_1 - f_2, f_2 - f_1, f_1 + f_2$ etc.

## 4.2 Réponse non-linéaire comportant un seul mode instable $\mu = 0.40$

La figure 4 montre les réponses issues de la GCHBM (a) et de l'intégration numérique (b) à $\mu = 0.40$. Bien que l'analyse de stabilité montre deux modes instables à $\mu = 0.40$, l'analyse non linéaire montre une réponse harmonique où seul un mode semble participer à la dynamique. La trajectoire des cycles limites est fermée et les spectres (c) ne montrent qu'une fréquence fondamentale $f_1$ et ses harmoniques $2f_1$ et $3f_1$ du mode M1. Le choix de 3 harmoniques est suffisant pour décrire le comportement non-linéaire du système. On voit très clairement l'absence du mode M2 car aucune de ses harmoniques pures ou couplées avec le mode M1 ne sont présentes. Le mode M1 semble occulter le mode M2 du fait de sa probable plus grande énergie. On voit bien l'intérêt d'effectuer une analyse non linéaire puisque cette dernière permet de détecter des comportements
a

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

b

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

c

<figure><img src="image_7.png" /><figcaption>Figure 4 – Cycles limites à μ = 0.40 : (a) GCHBM, (b) Intégration numérique. (c) Spectre de densité de puissance : - - GCHBM, — Intégration numérique</figcaption></figure>

qui n'apparaissent pas dans l'analyse linéaire classique. De plus, on remarque que les amplitudes extrêmes sont plus faibles pour μ = 0.40.

Le temps de calcul de la GCHBM est inférieur d'un facteur d'au moins 50 par rapport à une intégration temporelle classique et ne nécessite que quelques kilo-octets pour stocker les coefficients de Fourier du signal stationnaire. De plus, à partir d'un résultat de calcul, il est aisé d'obtenir la réponse d'un autre jeu de paramètres en réutilisant les résultats comme conditions initiales.

# 5 Conclusion

Ce papier présente une méthode originale qui permet de calculer les cycles limites de systèmes autonomes soumis à un ou plusieurs modes instables. Les réponses correspondantes peuvent être périodiques ou pseudo-périodiques selon le nombre de modes instables présents dans la dynamique du système. Le modèle considéré montre bien l'intérêt de cette méthode puisque non seulement elle permet d'affiner l'analyse en donnant les amplitudes de vibrations mais aussi les modes impliqués dans la dynamique du système. De plus elle présente des performances intéressantes et des analyses paramétriques peuvent être conduites plus rapidement.

# Références

[1] N.M. Kinkaid, O.M. O'Reilly, P. Papadopoulos. Automotive disc brake squeal. *Journal of Sound and Vibration*, 105-166, 2003.

[2] H. Ouyang, W. Nack, Y. Yuan, F. Chen. Numerical analysis of automotive disc brake squeal : a review. *Int. J. Vehicle Noise and Vibration*, 207-231, 2005.

[3] F. Massi, L. Baillet, O. Giannini, A. Sestieri. Linear and non-linear numerical approaches. *Mechanical Systems and Signal Processing*, 2374-2393, 2007.

[4] Y.-B. Kim, S.-K. Choi. A multiple harmonic balance method for the internal resonant vibration of a non-linear jeffcot rotor. *Archives of Computational Methods in Engineering - State of the art reviews*, 257-344, 2004.

[5] J.-J. Sinou, F. Thouverez, L. Jezequel. Methods to reduce non-linear mechanical systems for instability computation. *Archives of Computational Methods in Engineering - State of the art reviews*, 257-344, 2004.

[6] W.D. Iwan. A generalization of the concept of equivalent linearization. *International of Journal Non-Linear Mechanics*, 279-287, 1973.

[7] T.M. Cameron, J.H. Griffin. An Alternating Frequency/Time Domain Method for Calculating the Steady-State Response of Nonlinear Dynamic Systems. *Journal of Applied Mechanics*, 149-154, 1989.