# Vibration non-linéaire des structures avec frottement: application aux modèles possédant de nombreux ddls non-linéaires

Loic Salles

► To cite this version:

Loic Salles. Vibration non-linéaire des structures avec frottement: application aux modèles possédant de nombreux ddls non-linéaires. 11e colloque national en calcul des structures, CSMA, May 2013, Giens, France. hal-01722092

HAL Id: hal-01722092

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

Submitted on 2 Mar 2018

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 CC0 1.0 - Universal - International License
# Vibration non-linéaire des structures avec frottement: application aux modèles possédant de nombreux ddls non-linéaires

Loïc SALLES<fnref n="1" />

## Résumé —

Ce papier présente une méthode pour le calcul de la réponse forcée de structure possédant des non-linéarités de contact. La stratégie repose sur l'utilisation d'une méthode d'équilibrage harmonique avec procédure temps fréquence. Le problème complet dynamique est résolu de façon séquentiel avec une première étape statique puis la résolution d'un problème non-linéaire dynamique réduit autour de la position d'équilibre statique. La réponse forcée est construite à l'aide d'une procédure de continuation avec prédiction-correction. L'opération la plus couteuse est la résolution des systèmes linéaires à l'étape de correction. L'utilisation d'un solveur itératif GMRES est proposée et son efficacité étudiée.

**Mots clés** — Dynamique non-linéaires, HBM, frottement, synthèse modale, continuation, résolution de systèmes linéaires dense

## 1 Introduction

Le frottement tient un rôle important dans les systèmes mécaniques. Il peut être la conséquence de la nature même de la conception du système (roulement a bille) ou d'un fonctionnement anormal (contact entre rotor et stator par exemple). Le frottement peut être utile pour assurer un rôle mécanique lié à la conception. Dans cette étude nous nous intéressons en particulier aux systèmes assemblés. L'assemblage entre les éléments est souvent assuré par des liaisons comportant un ou plusieurs contacts avec frottement : liaisons boulonnées, attaches aube-disque, assemblage des bagues extérieur d'un roulement avec le carter...

<figure><img src="image_2.png" /><figcaption>Fig. 1 – Non-linéarités de contact dans un moteur d'avion</figcaption></figure>

Le contact avec frottement joue aussi un rôle important dans les vibrations des systèmes assemblés. Il modifie la raideur et surtout l'amortissement au niveau de l'interface entre les éléments mécaniques entrainant une diminution des amplitudes de vibration.

<footnote n="1">¹ Imperial College London, Vibration University Technology Centre, 1.salles@imperial.ac.uk</footnote>
Durant la phase de conception il est donc important de pouvoir modéliser les vibrations des systèmes avec contact et frottement. L'élément de contact avec frottement peut être intégré directement dans le modèle éléments finis (approche dite explicite) ou simplifié par l'utilisation d'un élément équivalent représentant tout le contact (élément implicite). Ces derniers comportent peu de degrés de libertés et permettent un calcul rapide de la réponse dynamique. Ils ne permettent pas toutefois de rendre compte de façon réaliste du comportement dans la zone de contact. Le contact avec frottement entraine une modification importante des gradients de contrainte dans les zones de contact. Pour pouvoir modéliser le comportement dans les zones de contact il faut utiliser un maillage raffiné dans ces zones. La résolution de problèmes éléments finis avec contact et frottement possédant de nombreux degrés de libertés est réalisable facilement aujourd'hui pour un problème statique tandis que le même problème de même taille en vibration reste un challenge.

## 2 équation de mouvement dynamique

Un système comprenant deux solides en contact est considéré. Le système dynamique est soumis à des forces volumiques centrifuges, des forces surfaciques aérodynamiques stationnaires et des excitations périodiques de pulsation $\omega$. L'équation de mouvement s'écrit :

$$[M]^1 \underline{\underline{U}}^1 + [C]^1 \underline{\underline{U}}^1 + [K]^1_T \underline{\underline{U}}^1 = \underline{F}_{\Omega}^1 + \underline{F}_T^1 + \underline{F}_{ex}^1(\omega,t) + \underline{F}_{NL}^1 ([B]\underline{\underline{U}}, [B]\underline{\underline{U}}) \quad (1)$$

$$[M]^2 \underline{\underline{U}}^2 + [C]^2 \underline{\underline{U}}^2 + [K]^2_T \underline{\underline{U}}^2 = \underline{F}_{\Omega}^2 + \underline{F}_T^2 + \underline{F}_{ex}^2(\omega,t) + \underline{F}_{NL}^2 ([B]\underline{\underline{U}}, [B]\underline{\underline{U}}) \quad (2)$$

où $[M]^i, [C]^i$ et $[K]^i_T$ sont respectivement les matrices de masse, d'amortissement et de raideur tangentielle. $[K]^i_T$ dépend de $\underline{\underline{U}}^i$ et se décompose en $[K]^i_T = [K]^i - [K]^i_{\Omega} + [K]^i_g (\underline{\underline{U}}^i)$ avec $[K]^i, [K]^i_{\Omega}$ et $[K]^i_g$ respectivement les matrices de raideur élastique, d'assouplissement centrifuge et de non-linéarités géométriques. $\underline{F}_{\Omega}^i, \underline{F}_T^i$ et $\underline{F}_{NL}^i$ sont les forces centrifuges, les forces surfaciques (forces aérodynamiques stationnaires) et les forces non-linéaires (contact avec frottement). Les forces de contact ne sont présent qu'à l'interface de contact entre les deux solides. Elles dépendent des déplacements relatifs entre des deux solides et en dynamique aussi des vitesses relatifs à l'interface. Les déplacements relatifs sont calculés par l'expression suivante :

$$[B]\underline{\underline{U}} = \underline{\underline{U}}_{int}^1 - \underline{\underline{U}}_{int}^2 \quad (3)$$

où $[B]$ est la matrice booléenne donnant les déplacements relatifs à l'interface à partir des déplacements absolus. La troisième loi de Newton donne la relation suivante entre les forces non-linéaires :

$$\underline{F}_{NL}^2 + \underline{F}_{NL}^1 = \underline{0} \quad (4)$$

$\underline{F}_{ex}^2(\omega,t)$ et $\underline{F}_{ex}^2(\omega,t)$ sont les forces d'excitation périodiques de fréquence $\omega$. Dans le cas d'une application turbomachine la fréquence d'excitation $\omega$ est le plus souvent reliée à la fréquence de rotation $\Omega$.

Le système (2) peut être résolu par intégration numérique pour chaque fréquence. Ces approches sont très chères numériquement. Comme la réponse est supposée périodique, des approches basées sur la méthode de Galerkin sont préférées. La solution est projetée sur une série de Fourier tronquée à l'ordre $N_h$. Les systèmes (2) devient :

$$[Z]\tilde{X} = \tilde{F}_{ex} - \tilde{F}_g - \tilde{F}_{nl}(\underline{[B]}\tilde{X}) \quad (5)$$

où $[Z]$ est la matrice de raideur dynamique définie par :

$$[Z] = \text{blkdiag} \left( [K], \begin{bmatrix} [K] - \omega^2[M] & \omega[C] \\ -\omega[C] & \text{matK} - \omega^2[M] \end{bmatrix}, \dots, \begin{bmatrix} [K] - N_h^2\omega^2[M] & N_h\omega[C] \\ -N_h\omega[C] & \text{matK} - N_h^2\omega^2[M] \end{bmatrix} \right) \quad (6)$$

et

$$[K] = \begin{bmatrix} [K]^1 & [0] \\ [0] & [K]^2 \end{bmatrix}, \quad [M] = \begin{bmatrix} [M]^1 & [0] \\ [0] & [M]^2 \end{bmatrix}, \quad [C] = \begin{bmatrix} [C]^1 & [0] \\ [0] & [C]^2 \end{bmatrix}$$
$\tilde{F}_{ex}$ et $\tilde{F}_{g}$ sont les vecteurs des coefficients de la serie de Fourier des forces d'excitations et les termes constants des forces centrifuges et surfaciques.

La résolution directe du système (5) est difficilement envisageable pour des applications industrielles. En effet le calcul d'une réponse forcée nécessite de résoudre $N_{freq}$ (nombre de fréquences calculées) fois un système non-linéaire de taille $(2N_h + 1) * N_{ddl}$. Il faut donc réduire la taille du système à résoudre. Des méthodes de synthèse modale sont utilisées pour réduire la taille du système (5). Ces méthodes sont présentées dans la section suivante. Les méthodes de synthèse modale nécessitent de calculer les modes de la structures. Or pour des structures chargées statiquement, les modes dépendent de ce chargement. Il faut donc séparer le calcul en trois étapes [3] : étape statique non-linéaire, calcul des modes autour de la position d'équilibre statique avec création d'un modèle réduit et enfin résolution du problème de dynamique non-linéaire sur le modèle réduit.

## 3 Résolution séquentiel des problèmes statique et dynamique

Les déplacements sont séparés en partie statique et dynamiques $\text{vect}U = \underline{U}_s + \underline{u}_d$. Le problème non-linéaire statique s'écrit :

$$[\mathbf{K}]_T^1 (\underline{U}_s^1) \underline{U}_s^1 = \underline{F}_{\Omega}^1 + \underline{F}_{\Gamma}^1 + \underline{F}_{ex}^1 (\omega, t) + \underline{F}_{NL}^1 ([\mathbf{B}] \underline{U}_s) \quad (7)$$

$$[\mathbf{K}]_T^2 (\underline{U}_s^2) \underline{U}_s^2 = \underline{F}_{\Omega}^2 + \underline{F}_{\Gamma}^2 + \underline{F}_{ex}^2 (\omega, t) + \underline{F}_{NL}^2 ([\mathbf{B}] \underline{U}_s) \quad (8)$$

### 3.1 Résolution du problème statique

Le système (8) est résolu de façon itérative en utilisant un solveur de Newton et une rampe pour appliquer le chargement lié à la rotation de la structure. À la convergence la position d'équilibre est obtenue. Un point important concerne le traitement des forces non-linaires liées au contact au niveau des interface entre les deux solides. Trois stratégies sont envisageables :

- Il est possible de considérer les déplacements relatifs à l'interface de contact liés au chargement statique petits et de les négliger. Les deux faces de l'interface sont alors collées durant tout le chargement statique. Le contact avec frottement est alors pris en compte que durant la phase dynamique sur le modèle réduit.

- Une deuxième approche plus rigoureuse que la première consiste à traiter le contact avec frottement à l'étape statique avec la même formulation que celle utilisée à l'étape dynamique. La formulation du contact repose sur une hypothèse de petit déplacement et un maillage coïncident au niveau de l'interface. Dans le cadre du code utilisé en dynamique le contact est modélisé par une méthode de pénalité.

- La dernière stratégie consiste à considérer un problème de contact en grand déplacement avec l'utilisation d'éléments surfaciques de contact. Les éléments hybrides de contact avec un solveur de Newton généralisée peuvent être utilises pour résoudre ce problème en statique. Aujourd'hui le code de dynamique non-linéaire en fréquentiel ne peut traiter que des éléments de contact nœuds a nœuds. Cela implique que le maillage doit être modifiée après l'étape statique pour faire coïncider les nœuds en contact des deux solides. Le développement d'éléments surfaciques de contact dans le code dynamique est envisagé dans le futur. Aujourd'hui cette approche pour le calcul statique est utilisée avec des éléments de contact nœud a nœud mais on perd l'égalité entre la projection de l'état d'équilibre statique sur la base réduite et l'état statique initial utilisé dans le problème dynamique.

Dans le cadre de cette étude la première approche est utilisée. On note $\lambda^s$ les forces de contact statiques, le problème (8) s'écrit :

$$\begin{bmatrix} [\mathbf{K}_T] & [\mathbf{B}]^t \\ [\mathbf{B}] & [0] \end{bmatrix} \begin{Bmatrix} \underline{U}_s \\ \underline{0} \end{Bmatrix} = \begin{Bmatrix} \underline{F}_{\Omega} + \underline{F}_{\Gamma} \\ \underline{0} \end{Bmatrix} \quad (9)$$

Une fois la position d'équilibre statique obtenue une analyse modale du système en considérant les interfaces de contact collées peut être effectuer. Le système (9) est résolu dans le programme élément fini code-aster en utilisant l'opérateur `STA_NON_LINE_`. On résout ensuite le problème au valeur propre suivant :
$$[K_T] \underline{X} - \omega^2 [M] \underline{X} = \underline{0} \qquad (10)$$

Cette analyse modale permet de connaitre la position (fréquence) du mode linéaire. La construction de la réponse forcée se fera en fonction des fréquences des modes linéaires, en effet les effets non-linéaires liés au contact ne font pas varier de façon trop importante les fréquences de résonance.

## 3.2 Réduction du système dynamique

Différentes approches sont possibles pour réduire la taille du système dynamique. Trois approches sont présentées, tout d'abord le modèle dynamique peut être projeté sur une base modale réduite. Les déplacements sont représentés par :

$$\underline{u}_d = [\underline{\Phi}_L] \underline{\eta} \qquad (11)$$

$[\underline{\Phi}_L]$ est la matrice des modes à interfaces libres.

Cette approche est la moins performante puisqu'elle nécessite un grand nombre de mode qui est souvent bien supérieur au nombre de degrés de liberté de contact pour converger. Pour diminuer l'erreur entre le modèle réduit et le modèle initial on ajoute les effets statiques en calculant la matrice de flexibilité aux nœuds ou des efforts sont appliquées. Dans notre cas, cela correspond à des nœuds où le contact intervient. Il n'est pas indispensable de prendre en compte tous les nœuds en contact. Bien entendu la précision de la base réduite augmente quand le nombre de nœuds retenus pour calculer la matrice de flexibilité tend vers le nombre de nœuds en contact.

$$\underline{u}_d = [\underline{\Lambda}] \underline{q} = [\underline{\Phi}_L] \underline{\eta}_L + [\underline{G}_r] \underline{\lambda}_r \qquad (12)$$

$\underline{q}$ est le vecteur des déplacements généralisés. $[\underline{G}_r]$ est la matrice de flexibilité obtenu en résolvant un problème linéaire avec membre droit de taille (RHS) $N_r$ avec une force unitaire sur chaque ddl retenu. $[\underline{G}_r]$ peut être définie par l'expression suivante :

$$[\underline{G}]_r = [\underline{K}]^{-1} [\underline{I}]_{[N_{ddl}*N_r]} \qquad (13)$$

Si la structure est libre il faut filtrer les modes de corps rigide puisque la matrice $[\underline{K}]$ est singulière. Cette méthode est souvent appelée méthode de Mac-Neal.

Une autre approche consiste à utiliser des modes à interface fixe comme base de Ritz enrichi avec les déplacements aux interfaces. Cette méthode est connue sous le nom de méthode de Craig et Bampton[1]. Comme pour la méthode de Mac-Neal il n'est pas indispensable de retenir tous les nœuds en contact. L'équation pour les déplacements s'écrit :

$$\underline{u}_d == [\underline{\Lambda}] \underline{q} = \left\{ \begin{matrix} \underline{u}_o \\ \underline{u}_r \end{matrix} \right\} = \left[ \begin{matrix} [\underline{\Phi}_F] & [\underline{\Psi}] \\ [\underline{0}] & [\underline{I}] \end{matrix} \right] \left\{ \begin{matrix} \underline{\eta}_F \\ \underline{u}_r \end{matrix} \right\} \qquad (14)$$

$\underline{q}$ est le vecteur des déplacements généralisés. $\underline{\Phi}_F$ est la matrice des modes à interface fixe et $\underline{\Psi} = -K_{oo}^{-1}K_{or}$ est la matrice des modes statiques. La matrice des modes statiques est calculée en résolvant un problème linéaire avec multiple membres droits.

Que ce soit pour la méthode McNeal ou de Craig-Bampton la résolution du problème linéaire pour obtenir les modes statiques est le goulot d'étranglement d'un point de vue mémoire vive utilisée. En effet quand le nombre de nœuds retenus devient très important la mémoire nécessaire pour construire les modes statiques peut dépasser les limites du système informatique utilisé alors que la construction de la réponse forcée sur le système réduit ne demandera comparativement que peu de mémoire. Les modes statiques sont calculés en utilisant l'opérateur MODE_STATIQUE dans code-aster.

La projection sur la base réduite du système (2) perturbé autour de la position statique d'équilibre donne le système dynamique suivant :

$$[\underline{M}_{red}^1] \underline{\ddot{q}}_1 + [\underline{C}_{red}^1] \underline{\dot{q}}_1 + [\underline{K}_{red}^1] \underline{q}_1 = [\underline{\Lambda}^1]^{\prime} \underline{F}_{ex}^1(\omega, t) + \underline{\Delta F}_{NL}(\underline{\lambda}_s, \underline{u}_c, \underline{\dot{u}}_c) \qquad (15)$$

$$[\underline{M}_{red}^2] \underline{\ddot{q}}_2 + [\underline{C}_{red}^2] \underline{\dot{q}}_2 + [\underline{K}_{red}^2] \underline{q}_2 = [\underline{\Lambda}^2]^{\prime} \underline{F}_{ex}^2(\omega, t) - \underline{\Delta F}_{NL}(\underline{\lambda}_s, \underline{u}_c, \underline{\dot{u}}_c) \qquad (16)$$
$\underline{\Delta F_{NL}}$ est le vecteur des coefficients de Fourier des forces dynamiques de contact. Leur calcul est présenté dans la section suivante.

Les déplacements relatifs aux interfaces de contact sont obtenus en utilisant l'expression suivante :

$$ \underline{u_c} = [\underline{B}^1] [\underline{\Lambda}^1] \underline{q}_1 - [\underline{B}^2] [\underline{\Lambda}^2] \underline{q}_2 \quad (17) $$

Finalement le problème (16) est projeté sur la base trigonométrique pour obtenir un système fréquentiel.

$$ [\underline{Z}_{red}^1] \underline{q}_1 = [\underline{\tilde{\Lambda}}^1] \underline{\tilde{F}}_{ex}^1(\omega, t) + \underline{\Delta \tilde{F}}_{NL}(\underline{\lambda}_s, \underline{u}_c, \underline{\dot{u}}_c) \quad (18) $$

$$ [\underline{Z}_{red}^2] \underline{q}_2 = [\underline{\tilde{\Lambda}}^2]^t \underline{\tilde{F}}_{ex}^2(\omega, t) - \underline{\Delta \tilde{F}}_{NL}(\underline{\lambda}_s, \underline{u}_c, \underline{\dot{u}}_c) \quad (19) $$

Les coefficients de Fourier des forces de contact ne peuvent pas être calculées directement en fréquentiel il faut repasser dans le domaine temporel. Le code utilise une approche discrète basée sur des transformées de Fourier discrètes inverse pour obtenir les déplacements et les forces dans le domaine temporel sur une période de vibration.

## 4 Prise en compte du contact avec frottement

Le contact avec frottement est traité en utilisant une formulation de contact régularisée élastique. Les déplacements $u_c = u_N \underline{n} + \underline{u}_T$ sont obtenus à partir des vecteurs harmoniques $\underline{U}_c$.

Selon l'axe normal on introduit entre les deux solides en contact la raideur normale $k_n$ tandis que dans le plan tangentiel on introduit la raideur $k_f$. Ces raideurs sont ajoutées au départ pour des considérations numériques. Toutefois d'un point de vue physique elles peuvent représenter le comportement élastoplastiques des aspérités de l'interface de contact [2] ou du lubrifiant solide. Une procédure de prédiction-correction est utilisée, la force normale prédite vaut :

$$ \lambda_{N}^{\text{pre},k} = \lambda_{s}^{N} + k_{n}u_{N} \quad (20) $$

La force tangentielle prédite s'écrit :

$$ \underline{\lambda}_{T}^{\text{pre},k} = \underline{\lambda}_{s}^{T} + k_{t} \left( \underline{u}_{T}^{k} - \underline{z}^{k-1} \right) \quad (21) $$

où $\underline{z}^{k-1}$ est le vecteur des déplacements relatifs tangentiel au point de contact au pas de temps $k-1$.

Les efforts de contact prédits sont corrigés de telle sorte que les contraintes soient vérifiées. Trois états sont possibles : séparation,adhérence et glissement.

1. Séparation : $\lambda_{pre}^{k,N} \le 0$, la force prédite est une force de traction, ce qui signifie que le contact est perdu au niveau de l'élément de contact. La force de contact doit donc être nulle, ce qui nécessite

2. Adhérence : $\lambda_{pre}^{k,N} > 0$ et $\|\lambda_{pre}^{k,T}\| < \mu|\lambda_{pre}^{k,N}|$, dans ce cas le déplacement normal est nul et le déplacement relatif tangentiel est constant. Ceci implique

$$ \{ \underline{z}^k = \underline{z}^{k-1} \} \quad (23) $$

3. Glissement : $\lambda_{pre}^{k,N} > 0$ et $\|\lambda_{pre}^{k,T}\| \ge \mu|\lambda_{pre}^{k,N}|$, à nouveau, le déplacement relatif normal est nul. La force normale est donc donnée par

$$ \lambda^{k,N} = \lambda_{pre}^{k,N} \quad (24) $$

La correction de la force de contact tangentielle s'appuie sur la connaissance de l'effort normal par (24) et sur l'hypothèse que la vitesse relative tangentielle est définie de manière implicite :

$$ \underline{v}^{k,T} = \frac{\underline{z}^{k,T} - \underline{z}^{k-1,T}}{t_k - t_{k-1}} \quad (25) $$
D'après les lois de frottement de Coulomb, la force de contact tangentielle doit avoir une norme égale à $\mu|\lambda_{pre}^{n,N}|$ et être colinéaire à $\underline{v}^{k,T}$ et de même sens ($\underline{\lambda}^k$ est défini comme l'opposé des forces de contact et agit donc dans le sens de la vitesse relative). Ces conditions sont résumées par

$$ \underline{\lambda}^{k,T} = \mu |\lambda_{pre}^{k,N}| \frac{\underline{v}^{k,T}}{\|\underline{v}^{k,T}\|} \qquad (26) $$

D'après (21) on peut écrire

$$ \underline{\lambda}_{pre}^{k,T} - \underline{\lambda}^{k,T} = k_t (\underline{\xi}^k - \underline{\xi}^{k-1}). \qquad (27) $$

En combinant (26) et (27) on trouve

$$ \underline{\lambda}_{pre}^{k,T} = \left( \frac{\mu |\lambda_{pre}^{k,N}|}{\|\underline{v}^{k,T}\|} + k_t (t_k - t_{k-1}) \right) \underline{v}^{k,T}, \qquad (28) $$

ce qui montre que $\underline{\lambda}_{pre}^{k,T}$ et $\underline{v}^{k,T}$ sont colinéaires et de même sens et donc

$$ \underline{\lambda}^{k,T} = \mu |\lambda_{pre}^{k,N}| \frac{\underline{\lambda}_{pre}^{k,T}}{\|\underline{\lambda}_{pre}^{k,T}\|} \qquad (29) $$

En reportant (29) dans (27), on trouve la relation qui permet de calculer $\underline{\xi}^k$ :

$$ \underline{\xi}^k = \underline{\xi}^{k-1} + \frac{1}{k_t} \underline{\lambda}_{pre}^{k,T} \left( 1 - \mu \frac{|\underline{\lambda}_{pre}^{k,N}|}{\|\underline{\lambda}_{pre}^{k,T}\|} \right). \qquad (30) $$

Finalement les forces de contact statique sont souhaitées aux forces de contact pour obtenir les forces dynamiques de contact.

$$ \underline{\Delta f}_{nl}^{k} = f_{nl}^{k} - \underline{\lambda}_{s} \qquad (31) $$

Une fois les forces de contact calculées en temporel l'étape finale consiste à transformer les forces de contact dans le domaine fréquentiel à l'aide d'une transformée de Fourier discrète (DFT).

# 5 Procédure de continuation

Le système (19) est réécrit en utilisant la fonction de transfert [2] $H(\omega)$ dont chaque composante harmonique $n$ est définie par :

$$ [H]_n = \sum_{j=1}^{N_m} \frac{\Phi_j^R \Phi_j^{R,T}}{\omega_j^2 - n^2\omega^2 + i2n\xi_j\omega_j\omega} \qquad (32) $$

$\Phi_j^R$ est la base modale entière du modèle réduit ($K_{red}$ et $M_{red}$), $\omega_j$ et $\xi_j$ sont la pulsation propre et le coefficient d'amortissement modal du mode $j$.

Le système (19) après transformation se ramène à l'annulation du résidu défini par :

$$ \underline{f} = \tilde{X}_n - [\underline{H}]_{nl} \underline{F}_{ex}^l - [\underline{H}]_{nn} \left( \underline{F}_{ex}^n + \underline{\lambda}(\tilde{X}_n) \right) \qquad (33) $$

où $l$ et $n$ réfèrent respectivement aux degrés de liberté linéaires et non-linéaires du système.

Le calcul de la réponse forcée est réalisé en annulant la fonction résidu (33) pour chaque pulsation $\omega$ nécessaire dans la procédure de continuation. Chaque étape de la procédure de continuation peut être divisée en quatre sous-étapes comme représentées sur la figure 2 :

1. À partir de la solution convergée à la pulsation $\omega^i$ les gradients du résidu en fonction de $\tilde{X}$ et $\omega$ sont calculés.

2. Un prédicteur de la solution pour la pulsation $\omega^{i+1}$ est construit en résolvant le système :

$$ \frac{D \underline{f}}{D \underline{\tilde{X}}} \Delta \underline{\tilde{X}}^{pre} = - \frac{D \underline{f}}{D \underline{\omega}} \Delta \underline{\omega}^i $$
<figure><img src="image_3.png" /><figcaption>Fig. 2 – étapes de la procédure de continuation pour construire la FRF</figcaption></figure>

3. Le résidu pour la fréquence $\omega^{i+1} = \omega^i + \Delta\omega^i$ est calculé avec le prédicteur $\tilde{X}^{pre} = \tilde{X}^i + \Delta\tilde{X}^{pre}$, si la convergence est suffisante on retourne a l'étape 1, sinon on lance une procédure de correction basée sur un solveur de Newton.

4. Le correcteur comprend la recherche d'une direction de Newton :

$$ \frac{Df}{D\tilde{X}} \Delta \underline{p}^k = -f(\tilde{X}^{i,k}) $$

où $k$ est l'itération du correcteur. Une fois la direction de Newton trouvée une recherche linéaire est effectuée pour minimiser le résidu $f$.

Que ce soit à l'étape de prédiction ou durant les itération du solveur de Newton un système linéaire dense doit être résolu. La résolution par un solveur directe basé sur une décomposition LU des matrices jacobienne est présente dans la version actuelle du code. Pour des systèmes avec de nombreux degrés de libertés cette approche prend beaucoup de temps CPU. Pour cette raison un nouveau solveur itératif basé sur la méthode GMRES [4] a été implémenté dans le code. GMRES a été retenu puisque la matrice jacobienne est non-symétrique, conséquence des non-linéarités de contact. Une comparaison entre les deux approches est donnée dans la section suivante.

# 6 Exemples numériques

La méthode présentée dans ce papier est appliquée sur un modèle de flasque avec boulons (Fig.3). Le contact est discrétisé par 600 éléments de contact. Un nœud est retenu pour excitation et un autre pour l'observation des amplitudes de vibration. La réponse fréquentielle autour du deuxième mode est présentée sur la Fig.3. Les résultats ont été obtenu avec des élément 3D (f3pen) dont la formulation a été présentée dans ce papier et un élément hybride 2D/3D [2] basé sur la séparation du problème tangentiel de contact suivant la direction $T_x$ et $T_y$. L'élément hybride 2D/3D repose sur un calcul analytique des coefficients de Fourier qui nécessite de trouver les temps de transition entre les différents états de contact pendant la période de vibration. Les tests numériques ont montrés que la formulation 3D basée sur un calcul discret en temps des coefficient de Fourier converge aussi vite que la formulation hybride 2D analytique. La formulation 3D permet une meilleure prédiction des niveaux vibratoires et des amortissement aux interfaces de contact.

L'analyse du temps nécessaire a chaque étape de la procédure de continuation a montré que l'étape la plus gourmande en temps de calcul était la factorisation de la matrice jacobienne (cf Fig.4 (a)). Les tests numériques effectués sur ce modèle avec le solveur directe et le solveur itératif GMRES ont montré un gain important en performance pour la version séquentiel du code mais cette avantage diminue avec la version parallèle du code. L'implémentation utilisée de GMRES ne semble pas profiter de la parallélisation (cf Fig.4 (b)).
<figure><img src="image_4.png" /><figcaption>Fig. 3 – Modèle EF d’un flasque vissé et réponse forcée autour du mode 2 avec deux types d’élément de contact</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4 – Analyse du temps utilisé par chaque étape de la procédure de continuation (a) et performance des solveurs linéaires</figcaption></figure>

# 7 Conclusions et perspectives

Une méthode d'obtention de la réponse forcée en présence de contact avec frottement a été présentée. Cette méthode repose sur une résolution en deux étapes : résolution d'un problème statique et résolution du problème dynamique sur un modèle réduit. Le problème dynamique est résolu autour de la position statique. Cette approche est adapte a des structures comportant de nombreux ddls. Le contact avec frottement est modélisé en dynamique de façon a conserver l'équivalence avec la formulation contact utilisée en statique. Cette formulation est plus rigoureuse que la plupart des formulations utilisées par d'autre auteur. La méthode a été appliquée sur un modèle large qui a permis de tester l'implémentation d'un nouveau solveur linéaire itératif basé sur l'algorithme GMRES. Des futurs travaux sont à mener pour améliorer ce solveur et utiliser les capacités de recyclage des sous-espaces de Krylov pour réduire le nombre d'itération du solveur linéaires a chaque fréquence.

# Références

[1] Craig, R.R., *Coupling of substructures for dynamic analysis*, AIAA journal, 6, 1313-1319, 1968, doi: 10.2514/3.4741

[2] Petrov, E.P., *A High-Accuracy Model Reduction for Analysis of Nonlinear Vibrations in Structures With Contact Interfaces*, Journal of Engineering for Gas Turbines and Power, 133, 10, 2011, doi: 10.1115/1.4002810

[3] Temis, Y.M and Salles, L., *Efficient CSM/CFD Method for Aeroelasticity Stability analysis of turbomachine blades*, proceeding of 14th ISROMAC, 2012

[4] Saad, Y. *Iterative methods for sparse linear systems*, Society for Industrial and Applied Mathematics, 2003