# Présentation d'un nouvel élément fini de contact

Etienne Arnoult, Isabelle Guilloteau, Bernard Peseux, Jérôme Bonini

► To cite this version:

Etienne Arnoult, Isabelle Guilloteau, Bernard Peseux, Jérôme Bonini. Présentation d'un nouvel élément fini de contact. Mechanics & Industry, 2001, 1 (2), pp.33-42. hal-01387692

HAL Id: hal-01387692

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

Submitted on 26 Oct 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
# Présentation d'un nouvel élément fini de contact

E. Arnoult<sup>a*</sup>, I. Guilloteau<sup>a</sup>, B. Peseux<sup>a</sup>, J. Bonini<sup>b</sup>

<sup>a</sup> Laboratoire Mécanique et Matériaux, École Centrale Nantes, BP 92101, 44321 Nantes cedex 3, France

<sup>b</sup> Département Méthodes, SNECMA Villaroche, 77550 Moissy-Cramayel, France

*Abstract* — Cet article présente une élément de contact spécifique, destiné à être implémenté dans un code de calcul par éléments finis implicite. La détection du contact est réalisée par mise en œuvre d'une méthode analytique basée sur l'utilisation des fonctions B-spline. Le traitement du contact est fait par la méthode des pénalités par élément de contact avec une raideur spéciale. © 2001 Éditions scientifiques et médicales Elsevier SAS

*Abstract* — A new contact finite element. This paper will present a new way to deal with contact within a finite element code. The detection of contact is realised with an analytical method using B-spline functions, and the analysis of the contact is done with a special finite element using a modified penalty method. © 2001 Éditions scientifiques et médicales Elsevier SAS.

## 1. INTRODUCTION

L'étude des turboréacteurs passe par une phase de modélisation numérique indispensable pour essayer de limiter le plus possible le nombre d'essais — coûteux et destructifs — réalisés sur banc réel. La façon dont le contact est géré pendant ces phases de simulation a beaucoup d'importance car il conditionne fortement les résultats : une mauvaise évaluation de l'effort de contact et sa mauvaise reprojection sur la structure peuvent conduire à des comportements simulés très éloignés de ce qui peut être observé expérimentalement.

Les développements qui sont présentés dans cette communication sont donc liés à ces deux points : il s'agit d'une part d'essayer de localiser avec précision les lieux de contact entre les structures (étude géométrique objet de la section 2), et d'autre part d'évaluer l'effort de contact qui en résulte (étude structurale, objet de la section 3). Les méthodes présentées sont alors regroupées dans un élément fini qui a été implémenté dans le code SAMCEF : les résultats obtenus sur un cas simple sont présentés dans la section 4.

## 2. RECHERCHE ANALYTIQUE DU CONTACT

### 2.1. Contexte

L'étude dynamique des contacts se produisant à l'intérieur des turboréacteurs a fait l'objet de nombreuses recherches. Le traitement numérique de ce problème a donné lieu au développement d'algorithmes spécifiques, souvent destinés à être intégrés à des codes éléments finis, dans le but d'essayer d'améliorer les simulations de mouvements complexes entre plusieurs pièces pouvant interagir. Cependant, les problèmes de facétisation, inhérents à la discrétisation spatiale des structures, peuvent limiter les performances de ces algorithmes. Les deux formes de difficultés que peut provoquer la facétisation sont d'une part un problème de discontinuité du vecteur normal entre les éléments (figure 1(a)), et d'autre part un problème d'approximation des courbures des géométries des structures (figure 1(b)).

Les deux palliatifs les plus utilisés pour contourner ces difficultés consistent, d'une part à attribuer à chaque nœud du maillage un vecteur normal « moyen », calculé en faisant la somme des vecteurs normaux des éléments auxquels le nœud appartient, et d'autre part à ajouter des
<figure><img src="image_2.png" /><figcaption>Figure 1. Problèmes liés à la discrétisation spatiale : (a) discontinuité du vecteur normal ; (b) mauvaise approximation de la courbure.</figcaption></figure>

nœuds au maillage pour obtenir une meilleur approximation de la courbure.

Ces solutions ne sont pourtant pas entièrement satisfaisantes. Dans le premier cas, sans même parler des approximations nécessaires pour obtenir un vecteur normal pour une facette qui compte plus de quatre nœuds, le problème de discontinuité du vecteur normal entre les éléments n'est pas résolu. Dans des configurations telles que la simulation dynamique d'un pallier, cette discontinuité peut s'avérer très pénalisante. Dans le second cas, l'ajout de nœuds au maillage a pour double conséquence l'augmentation de la taille du modèle, et la diminution de la taille relative des éléments (ce qui peut entraîner une diminution du pas de temps du calcul), ce qui se traduit, à chaque fois, par une augmentation du temps de calcul nécessaire à la simulation du phénomène.

La méthode proposée ici a été, à l'origine, motivée par l'étude des problèmes de contact aubes/carter dans les turboréacteurs pour lesquels une mauvaise détection du contact peut conduire à des simulations inexploitables («blocage» numérique des pièces entre elles, convergence vers des comportements non conformes à des observations expérimentales).

## 2.2. Démarche

La méthode proposée vient s'insérer dans la recherche de contact au moment de la sélection des nœuds qui, n'ayant pas vérifié un certain critère, vont se voir attribuer un effort dit effort de contact, dont le calcul précis fera l'objet de la section 3.

Le processus général comprend tout d'abord très classiquement une phase de sélection de domaines pouvant éventuellement entrer en contact sur les deux structures. Des entités géométriques sont ensuite construites à partir des nœuds des domaines retenus. Dans le cadre très

particulier de l'étude des contacts aubes/carter, les domaines en question seront d'une part une partie de sommet d'aube, qui sera ensuite approché par une courbe gauche, et d'autre part une partie du carter, qui sera approché par une surface.

La recherche de contact se fait alors avec ces deux entités géométriques, à l'aide d'une méthode récursive. Une fois les points de contact identifiés, il reste à reprojeter les résultats sur le maillage initial, c'est-à-dire qu'il faut pouvoir fournir à l'algorithme qui calcule l'effort de contact les grandeurs nécessaires pour chaque élément impliqué.

Dans la suite, nous supposerons que la première étape du processus, c'est-à-dire la sélection des nœuds candidats au contact, a été déjà réalisée. Il s'agit donc de construire, à partir de ces groupes de nœuds, des entités géométriques présentant un caractère de continuité suffisant pour s'affranchir des problèmes de facétisation. La théorie des splines, développée pour beaucoup par de Boor à la fin des années 70 [1, 2], présente une souplesse d'utilisation qui justifie son choix.

### 2.2.1. Modélisation des courbes gauches

Les courbes gauches seront représentées de façon paramétrique : chaque point P de la courbe est repéré par ses trois coordonnées x, y et z qui dépendent d'un même paramètre t, dont l'intervalle de définition doit être défini avec soin [3]. Pour toute valeur de t appartenant à cet intervalle, l'expression vectorielle de la courbe approchant un groupe de npt points $P_i$ est alors :

$$c(t) = \sum_{i=0}^{N} Q_i B_{ni}(t) \quad (1)$$

Où les $B_{ni}$ sont des fonctions de base de l'espace des splines appelées B-splines. Ces fonctions sont en fait
des polynômes de degré n définis par morceaux sur l'intervalle paramétré. Les N + 1 point $Q_i$ sont appelés points de contrôles de la courbe. Trois méthodes sont fréquemment utilisées pour déterminer ces points de contrôle :

* utiliser directement les points $P_i$ du jeu de données (méthode directe);

* calculer les $Q_i$ de façon à ce que la courbe $c(t)$ passe exactement par les points $P_i$ (interpolation);

* calculer les $Q_i$ de façon que la somme du carré des distances de chaque point $P_i$ à la courbe soit minimale (lissage au sens des moindres carrés).

### 2.2.2. Modélisation des surfaces

De la même façon que les courbes, les surfaces sont représentées avec une méthode paramétrique, et les points de contrôle peuvent également être obtenus à partir de l'une des trois méthodes exposées plus haut. Les données du problème sont alors constituées d'un groupe de points $nptu \cdot nptv$ notés $P_{ij}$. Deux paramétrages, $u$ et $v$ sont alors nécessaires pour décrire la surface dont l'expression vectorielle est la suivante :

$$s(u, v) = \sum_{i=0}^{Nu} \sum_{j=0}^{Nv} Q_{ij} B_{nu,i}(u) B_{nv,j}(v) \quad (2)$$

Les points $Q_{ij}$ sont les $(Nu+1) \times (Nv+1)$ points de contrôle. Comme on peut le voir sur cette formulation, une surface spline est équivalente à une modélisation croisée de plusieurs courbes gauches : celles définies par les points $Q_{ij}$, $j$ étant constant, et celles définies par les points $Q_{ij}$, $i$ étant constant. Il est alors possible d'utiliser des B-splines de degrés différents dans les deux directions (ces degrés ont été notés $nu$ et $nv$ dans l'expression (2)).

La prise en compte de l'épaisseur de la surface (cas d'éléments de coque par exemple) se fait en translatant les nœuds d'une quantité égale à la demi-épaisseur dans une direction normale moyenne des éléments. Tous les calculs de recherche d'intersection sont donc menés à partir de cette surface translatée. Noter que cette opération nécessite la définition d'un intérieur et d'un extérieur de la surface (démarche classique pour les méthodes de recherche de contact).

## 2.3. Principe de la méthode

Disposant de la définition analytique d'une courbe et d'une surface, il reste à en trouver les éventuels points

<figure><img src="image_3.png" /><figcaption>Figure 2. Situation initiale.</figcaption></figure>

d'intersection. Les solutions généralement utilisées dans les codes éléments finis (par exemple LS-DYNA [4], PLEXUS [5]...) commencent par classer les nœuds en deux groupes, « maîtres » et « esclaves ». Une première boucle est alors réalisée sur tout ou partie du groupe des esclaves : pour chacun d'eux, un test les reliant aux éléments maîtres permet de déterminer si un certain critère est vérifié (critère de présence dans un domaine de l'espace où le contact est susceptible de se produire par exemple). Si des couples maître/esclave ont pu être mis à jour lors de cette première étape, un second test est alors réalisé pour confirmer ou invalider le critère de contact (test sur la positivité de certains volumes). La phase d'appariement entre maîtres et esclaves peut se faire de façon globale (approche dite hiérarchique) ou locale (approche par voisinage).

La méthode présentée ici conserve la première étape de sélection : la courbe est considérée comme l'entité esclave. La seconde phase est réalisée via l'utilisation des formulations splines. Compte-tenu de la non-linéarité de cette formulation, une solution purement analytique n'est pas envisageable dans le cas général. Un processus numérique reposant sur une méthode itérative a donc été mise en place.

Le résultat de la recherche d'intersection n'est défini que pour une précision $\varepsilon$ donnée : un « point de contact » sera en fait représenté par un couple de points, le premier appartenant à la courbe, le second à la surface, tels que la distance les séparant soit inférieure à la précision $\varepsilon$. Pour arriver à ce résultat, les intervalles paramétrés des splines, à savoir $t$ pour la courbe, et $u$ et $v$ pour la surface, sont discrétisés. La taille des incréments $dt$, $du$ et $dv$ obtenus dépendent de la valeur $\varepsilon$ [6]. Les bornes initiales des intervalles paramétrés peuvent être obtenues après un premier tri grossier à l'aide d'une approche hiérarchique par exemple.

Dans la suite $P_t$ désignera le point de la courbe associé au paramètre $t$, et $P_{uv}$ le point de la surface associé au couple $(u, v)$. La courbe sert de support à la recherche (voir figure 2).
<figure><img src="image_4.png" /><figcaption>Figure 3. Principe de la méthode.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4. Zone d'intersection.</figcaption></figure>

Pour chaque valeur du paramètre *t*, la surface est balayée, et pour chaque couple (*u*, *v*), la distance $P_t P_{uv}$ est calculée (voir figure 3).

Si cette distance devient inférieure à la précision $\varepsilon$, les paramètres sont conservés en mémoire (paramètres dits d'initiation). Les paramètres sont ensuite incrémentés. Dès que la distance redevient supérieure à $\varepsilon$, les paramètres sont à nouveau gardés (paramètres de terminaison). Un certain nombre de zones potentielles d'intersection entre les deux entités géométriques sont ainsi définies (voir figure 4).

Après avoir terminé un premier balayage, un autre balayage est effectué avec une précision inférieure à $\varepsilon$. L'opération est répétée jusqu'à l'obtention de la précision voulue.

Les premières simulations réalisées avec cette méthode [7] ont montré que cet algorithme n'était pas si coûteux en temps de calcul que sa structure pourrait le laisser croire. Les premiers tests ont en particulier souligné que la possibilité offerte de travailler avec des maillages « allégés » en éléments permettait de compenser l'augmentation du temps de calcul induit par la méthode.

## 3. EVALUATION DE L'EFFORT DE CONTACT

### 3.1. Contact ponctuel

Le contact unilatéral entre deux solides est généralement caractérisé par deux paramètres : la distance *g* (gap) entre les deux solides candidats au contact et la pression de contact $\lambda_n$. Trois conditions découlent de la loi de contact unilatéral :

* la condition de non pénétration $g \ge 0$,

* la condition de pression $\lambda_n \le 0$, qui caractérise l'absence de traction sur la zone de contact,

* une condition de compatibilité : $g\lambda_n = 0$, qui traduit le fait que les deux solides sont en contact ($g = 0$) ou séparés ($\lambda_n = 0$).

Dans le cas d'une discrétisation par éléments finis, avec l'hypothèse de l'élastodynamique linéaire et en l'absence d'amortissement, le système matriciel à résoudre est :

$$M\ddot{U} + KU = F^{\text{ext}} - F^{\text{contact}} \quad (3)$$

avec *M* matrice masse, *K* matrice raideur, *U* vecteur des déplacements nodaux, $\ddot{U}$ vecteur des accélérations nodales, $F^{\text{ext}}$ vecteur des forces extérieures généralisées, $F^{\text{contact}}$ vecteur des forces de contact généralisées.

La condition de non pénétration s'écrit $G_i \ge 0$ où le vecteur gap linéarisé est donné par :

$$G = G_0 + QU \quad (4)$$

avec $G_i$ composante du vecteur des distances, *Q* matrice contrainte, qui dépend de la discrétisation des surfaces en contact.

Les modèles mathématiques principaux qui permettent de calculer les forces de contact sont basés sur la méthode des multiplicateurs de Lagrange et sur la méthode des pénalités.

La méthode des multiplicateurs de Lagrange si elle permet de respecter exactement les conditions de non pénétration, requiert l'introduction d'inconnues supplémentaires, les multiplicateurs de Lagrange $\Lambda$, qui représentent les forces de contact et telles que :

$$F^{\text{contact}} = Q^{\text{T}}\Lambda \quad (5)$$

La méthode des pénalités autorise des pénétrations entre les solides. Les forces de contact, également notées $\Lambda$, sont proportionnelles à la distance de pénétration *G* et sont égales aux réactions dues à des ressorts fic-
tifs introduits entre les entités en contact (nœuds ou segments) :

$$F^{\text{contact}} = Q^{\text{T}} \Lambda \qquad (6)$$

$$\Lambda = K_C G \qquad (7)$$

$K_C$ est la matrice raideur de contact. La principale difficulté de la méthode des pénalités est due au choix de la raideur de contact. En effet, une faible valeur de raideur de contact conduit à de grandes valeurs de pénétration, tandis qu'une valeur élevée de raideur de contact rigidifie le système et engendre des oscillations et des problèmes de convergence. Une raideur de contact couramment utilisée dans le cas de solides en contact, est celle proposée par Hallquist et al. [8] :

$$k_C = f_S \frac{K A^2}{V} \qquad (8)$$

avec $K$, $A$ et $V$ qui sont respectivement le module de compressibilité, l'aire de la surface de contact et le volume de l'élément maître. $f_S$ est un coefficient arbitraire qui doit être adapté au problème de contact traité.

## 3.2. Raideur de contact cinématique

Un des objectifs de ce travail, est de développer un élément de contact et de l'implémenter dans le code éléments finis MECANO (module du logiciel SAMCEF) qui utilise le schéma d'intégration temporelle classique de Newmark :

$$U_{n+1} = U_n + \Delta t \dot{U}_n + \left(\frac{1}{2} - \beta\right) \Delta t^2 \ddot{U}_n + \beta \Delta t^2 \ddot{U}_{n+1} \qquad (9)$$

$$\dot{U}_{n+1} = \dot{U}_n + (1 - \gamma) \Delta t \ddot{U}_n + \gamma \Delta t \ddot{U}_{n+1} \qquad (10)$$

La méthode des multiplicateurs de Lagrange est délicate à mettre en œuvre dans un code implicite. Elle conduit à des difficultés numériques (oscillations parasites et non convergence du processus itératif de Newton-Raphson) liées au caractère non linéaire du phénomène de contact. De plus les résultats sont très sensibles aux valeurs des paramètres de Newmark [9, 10].

La méthode des pénalités a été retenue car elle plus simple de mise en œuvre. Mais le problème délicat concerne le choix de la raideur de contact. Une étude comparative de la méthode des multiplicateurs associée à un schéma implicite de Newmark [9] et à un schéma explicite des différences finies centrales [11], a conduit

à proposer une nouvelle raideur de contact. Dans le cas des petites perturbations, le schéma implicite de Newmark est utilisé pour une méthode de prédiction (sans contact)-correction. L'état d'équilibre issu de la phase de prédiction est :

$$\left(\frac{M}{\beta \Delta t^2} + K\right) \Delta U^* = F_{n+1}^{\text{ext}} - F_n^{\text{int}} + \frac{M}{\beta \Delta t^2} \left(\Delta t \dot{U}_n + \left(\frac{1}{2} - \beta\right) \Delta t^2 \ddot{U}_n\right) \qquad (11)$$

avec $\Delta U^*$ incrément de déplacement prédit au temps $t_{n+1}$, $F_n^{\text{int}}$ vecteur des forces interne au temps $t_n$.

La valeur de la fonction distance est $G_{n+1}^* = G_0 + Q_{n+1}^* U_{n+1}^*$, avec la matrice $Q_{n+1}^*$ qui représente le prédiction des contraintes au temps $t_{n+1}$. On suppose que toutes les composantes de $G_{n+1}^*$ sont négatives, donc il y a pénétration, et que $Q_{n+1}^* = Q_{n+1}^*$, c'est à dire que les entités en contact au temps $t_{n+1}$ sont celles qui ont été détectées pendant la phase de prédiction. Ensuite durant la phase de correction, les vitesses et déplacements sont calculés pour satisfaire :

$$G_{n+1} = G_{n+1}^* + Q_{n+1}^* \Delta U^* = 0 \qquad (12)$$

où $U_{n+1} = U_n + \Delta U^* + \Delta U^c$, l'exposant c correspond aux valeurs corrigées.

Le nouvel état d'équilibre s'exprime par :

$$M \ddot{U}_{n+1} + F_{n+1}^{\text{int}} = F_{n+1}^{\text{ext}} - Q_{n+1}^{\text{T}} \Lambda_{n+1} \qquad (13)$$

Et les forces de contact sont données par :

$$\Lambda_{n+1} = \left[ Q_{n+1} \left( \frac{M}{\beta \Delta t^2} + K \right)^{-1} Q_{n+1}^{\text{T}} \right]^{-1} G_{n+1}^{*} \qquad (14)$$

La correction de déplacement est égale à :

$$\Delta U^c = - \left( \frac{M}{\beta \Delta t^2} + K \right)^{-1} Q_{n+1}^{\text{T}} \Lambda_{n+1}$$

La méthode des « multiplicateurs de Lagrange anti-cipés » proposée par Carpenter et al. [11] conduit à une expression du vecteur des forces de contact similaire à celle obtenue par la méthode des multiplicateurs associée à un schéma explicite des différences finies centrales :

$$\Lambda_n = \left[ \Delta t^2 Q_{n+1} M^{-1} Q_{n+1}^{\text{T}} \right]^{-1} G_{n+1}^* \qquad (15)$$

Pour une correction du contact au même temps $t_{n+1}$, les vecteurs des forces de contact en implicite ou en implicite (relations (14) and (15)) peuvent être considérés
comme le produit d'une raideur de contact équivalente et de la distance prédite $G_{n+1}^*$. En considérant que la méthode des multiplicateurs de Lagrange à un comportement similaire à celui de la méthode des pénalités, on peut conclure que les problèmes de convergence de la méthode des multiplicateurs sont dus à une valeur trop grande de la raideur de contact équivalente.

Aussi une nouvelle raideur de contact est proposée. Elle est déduite, par troncature, de l'expression de la force de contact (14) :

$$K_C = \frac{1}{\beta \Delta t^2} [QM^{-1}Q^T]^{-1} \quad (16)$$

La « matrice d'influence du contact » $QM^{-1}Q^T$ est le résultat d'une condensation de toutes les masses candidates au contact. Cette nouvelle raideur de contact est appelée « raideur de contact cinématique ».

Cette raideur de contact cinématique, contrairement à celle de Hallquist [8], ne dépend pas des caractéristiques des matériaux. De plus le paramètre de Newmark $\beta$, qui joue un rôle important dans la précision du schéma d'intégration joue un rôle supplémentaire, celui de coefficient de dureté de choc.

L'exemple simple de l'impact de deux masses ponctuelles, initialement traité Chaudhary et Bathe [12], permet de déduire un coefficient de dureté du choc associé à une raideur de contact cinématique égal à [13] :

$$4\beta - \gamma - \frac{1}{2} \quad (17)$$

Ce coefficient dépend des deux paramètres de Newmark. S'il est positif, alors le contact est relâché trop tard et s'il est négatif, le contact est relâché trop tôt. Un point particulier intéressant correspond au schéma de la méthode des trapèzes, $\beta = 1/4$, $\gamma = 1/2$, qui correspond au schéma d'intégration le plus précis et qui donne un coefficient de restitution de choc égal à zéro. Une analyse semblable à été menée avec la méthode des multiplicateurs de Lagrange. Il s'ensuit une autre expression du coefficient de dureté du choc [13] : $2\beta - \gamma - 1/2$. Les valeurs précoisées par Chaudhary et Bathe [12] ($\beta = \gamma = 1/2$) donne un coefficient de dureté du choc nul, donc le contact est correctement relâché.

## 4. IMPLEMENTATIONS ET RESULTATS

Le nouvel élément de contact a été implémenté dans le code éléments finis SAMCEF via les routines FORTRAN de l'élément utilisateur du module MECANO.

Cet élément comprend donc une phase de détection analytique du contact, comme décrit dans le paragraphe 2, puis une phase de calcul des efforts utilisant la méthode des pénalités associée à la raideur dynamique de contact dont le calcul est présenté dans le paragraphe 3.

### 4.1. Test sur la raideur dynamique de contact

La résolution du problème de l'impact monodimensionnel de deux barreaux élastiques différents [15] permet d'analyser les performances de la raideur de contact cinématique. Les caractéristiques (sans dimension) de l'impact sont données *tableau I*.

Les résultats obtenus par la méthode des multiplicateurs sont comparés aux résultats analytiques. Ces comparaisons sont faites pour différentes valeurs des paramètres de Newmark :

* la méthode des trapèzes, $\beta = 0,25$ ; $\gamma = 0,50$, correspondant au modèle à accélération constante et vérifiant la condition cinématique du choc $4\beta - \gamma - 1/2 = 0$;

* un couple spécial, $\beta = 0,505$ ; $\gamma = 0,51$, vérifiant le coefficient de dureté de choc des multiplicateurs $2\beta - \gamma - 1/2 = 0$.

Les figures 6 et 7 présentent l'évolution du déplacement des sections de contact des deux barreaux, les figures 8 et 9 présentent les variations de vitesses des sections de contact et les figures 10 et 11 montrent l'évolution de la force de contact.

<table><caption>TABLEAU I
Caractéristiques des barreaux</caption><thead><tr><th></th><th>Barreau 1</th><th>Barreau 2</th></tr></thead><tbody><tr><td>Module d'Young</td><td>E<sub>1</sub> = 0,49</td><td>E<sub>2</sub> = 1,0</td></tr><tr><td>Densité</td><td>ρ = 1,0</td><td>ρ = 1,0</td></tr><tr><td>Vitesse de propagation des ondes</td><td>c<sub>1</sub> = 0,7</td><td>c<sub>2</sub> = 1,0</td></tr><tr><td>Aire de la section droite</td><td>S = 1,0</td><td>S = 1,0</td></tr><tr><td>Longueur</td><td>L = 10</td><td>L = 10</td></tr><tr><td>Nombre d'éléments finis</td><td>n<sub>1</sub> = 100</td><td>n<sub>2</sub> = 70</td></tr><tr><td>Vitesse initiale</td><td>V<sub>0</sub><sup>1</sup> = 0,1</td><td>V<sub>0</sub><sup>2</sup> = 0,0</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>Figure 5. Impact de deux barres.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Evolution du déplacement des sections de contact. Méthode des trapèzes ($\beta = 1/4$; $\gamma = 1/2$).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Evolution du déplacement des sections de contact. Paramètres de Newmark vérifiant $(2\beta - \gamma - 1/2) = 0$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Evolution de la vitesse des sections de contact. Méthode des trapèzes ($\beta = 1/4$; $\gamma = 1/2$).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Evolution de la vitesse des sections de contact. Paramètres de Newmark vérifiant $(2\beta - \gamma - 1/2) = 0$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 10.** Evolution de la vitesse des sections de contact. Méthode des trapèzes ($\beta = 1/4$; $\gamma = 1/2$).</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 11.** Evolution de la vitesse des sections de contact. Paramètres de Newmark vérifiant $(2\beta - \gamma - 1/2) = 0$.</figcaption></figure>
<table><caption>TABLEAU II
Performance de l'impact de deux barreaux</caption><thead><tr><th></th><th>Raideur<br>cinématique</th><th>Multiplicateur<br>de Lagrange</th></tr></thead><tbody><tr><td colspan="3">β = 0,25 ; γ = 0,50</td></tr><tr><td>Nombre d'itérations moyen</td><td>1,05</td><td>1,20</td></tr><tr><td>Max. d'iter. pour un pas</td><td>2</td><td>2</td></tr><tr><td>Nombre total d'itérations</td><td>316</td><td>358</td></tr><tr><td colspan="3">β = 0,505 ; γ = 0,51</td></tr><tr><td>Nombre d'itérations moyen</td><td>1,01</td><td>1,17</td></tr><tr><td>Max. d'iter. pour un pas</td><td>2</td><td>2</td></tr><tr><td>Nombre total d'itérations</td><td>302</td><td>352</td></tr></tbody></table>

<figure><img src="image_13.png" /><figcaption>Figure 12. Maillage général.</figcaption></figure>

L'utilisation de la nouvelle raideur de contact, pour la méthode des pénalités dans un code implicite, présente différents avantages par rapport à la méthode des multiplicateurs classique :

* le nombre d'itérations de la méthode de Newton-Raphson est plus faible (réduction de plus de 12%, voir tableau II);

* moins de problèmes de convergence;

* compatibilité avec toutes les valeurs classiques des paramètres de Newmark. En particulier la méthode des trapèzes fournit de bons résultats avec une diminution des oscillations. Les valeurs vérifiant $2β - γ - 1/2 = 0$, donnent de meilleurs résultats avec la méthode des multiplicateurs qu'avec la méthode des trapèzes [13], et fournissent d'excellents résultats avec le nouvel élément de contact. En effet, ces valeurs introduisent un amortissement numérique bénéfique pour réduire les oscillations sur les forces de contact et vitesses des sections de contact.

## 4.2. Recherche du contact

La recherche analytique de contact a été implémentée dans le nouvel élément de contact. Le test présenté va permettre de comparer les résultats donnés par différentes méthodes : il consiste en la rotation d'une aube avec balourd à l'intérieur d'un carter dont la base est encastrée. L'aube est modélisée avec un ressort généralisé de type Bushing, et le carter avec des éléments coques de Mindlin. L'élément de contact est utilisé pour prendre en compte les contacts entre le sommet de l'aube et le carter. L'ensemble du maillage est représenté sur la figure 12.

La première méthode utilisée pour la détection du contact sera désignée sous le terme de méthode « classique » : il s'agit de travailler directement sur les géométries linéaires des éléments, à savoir des segments pour

l'aube, et des parties de plans pour le carter. Les autres méthodes utilisées sont les formulations splines dans le cas de l'interpolation et du lissage.

Les grandeurs observées sur le cas test sont les valeurs de l'effort normal de contact en sommet d'aube, ainsi que le temps de calcul CPU nécessaire à la simulation de cinq rotations complètes de l'aube dans le carter. Plusieurs configurations ont été envisagées, suivant d'une part le choix de la méthode utilisée pour la recherche de contact (classique, interpolation ou lissage), et d'autre part selon la finesse du maillage, et plus particulièrement, selon le nombre d'éléments dans la circonférence du carter (72, puis 36, puis 18).

### 4.2.1. Force de contact

Les premiers résultats présentés concernent l'évolution de l'amplitude de la force de contact en sommet d'aube : cette amplitude est normalisée par rapport à l'amplitude maximale obtenue dans la configuration 72 éléments/méthode classique. Comme cela est visible sur la figure 13 (72 éléments dans la circonférence du carter), les deux méthodes splines donnent des résultats très voisins (moins de 3 % d'écart entre les deux méthodes) sur le premier pic de réponse.

La méthode classique, qui utilise une approximation linéaire des éléments, sous-estime systématiquement le rayon de la structure (la valeur exacte du rayon n'étant atteinte qu'aux nœuds du maillage). Il en résulte des pénétrations plus importantes qu'avec des formulations splines, et donc des amplitudes de l'effort normal également plus importantes (environ 30 % au dessus des résultats donnés par les méthodes splines). Il en résulte également un décalage entre l'apparition du contact entre les deux types de méthodes (le contact avec les méthodes splines se produit plus tardivement).
<figure><img src="image_14.png" /><figcaption>Figure 13. Discrétisation avec 72 éléments.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 14. Discrétisation avec 18 éléments.</figcaption></figure>

Ces deux différences entre les méthodes ne font que s'accentuer lorsque le nombre d'éléments dans la circonférence du carter diminue. Sur la figure 14, qui représente la configuration avec seulement 18 éléments, la méthode classique détecte une forte pénétration initiale qui induit un effort de contact deux fois plus grand que dans la configuration 72 éléments.

En représentant sur un même graphique les résultats obtenus avec une méthode spline d'interpolation dans les trois configurations de maillage (72, puis 36, puis 18 éléments), il est intéressant de remarquer une très faible sensibilité de la méthode par rapport à l'amplitude de la force de contact (voir figure 15).

### 4.2.2. Temps CPU

Cette seconde série de résultats concerne le temps de calcul CPU requis pour simuler cinq tours complets de rotation de l'aube. Tous les temps CPU présentés sur la figure 16 sont normalisés par rapport au temps CPU

<figure><img src="image_16.png" /><figcaption>Figure 15. Méthode par interpolation.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Figure 16. Comparaison du temps CPU</figcaption></figure>

nécessaire à la réalisation du calcul avec la méthode classique dans la configuration 72 éléments.

Il apparaît que les formulations splines sont plus coûteuses que la méthode classique (environ 25% supplémentaires dans la configuration 72 éléments), ce qui est attendu compte-tenu de la forme récursive de la recherche. Il est néanmoins important de noter qu'une interpolation avec 36 éléments est moins coûteuse qu'une méthode classique avec 72 éléments (30% de différence environ), pour un résultat plus stable.

## 5. CONCLUSION

Deux méthodes dédiées au traitement du contact ont été développées dans un code éléments finis implicite. La première repose sur une approche géométrique pour déterminer avec précision les points de contact, et ainsi fournir des grandeurs affranchies des imprécisions induites par la facétisation. La seconde propose une nouvelle approche du calcul des efforts de contact, à par-
tir d'une méthode des pénalités qui utilise une nouvelle raideur dynamique de contact conduisant à des résultats meilleurs qu'avec des multiplicateurs de Lagrange.

## RÉFÉRENCES

[1] De Boor C., *A Practical Guide to Splines*, Applied Mathematical Sciences, Vol. 27, Springer Verlag, 1978.

[2] Numberger G., *Approximation by Spline Functions*, Springer Verlag, 1989.

[3] Daniel M., *Modélisation de courbes et surfaces par des B-splines* — Application à la conception et à la visualisation de formes, Thèse de doctorat, École Centrale de Nantes, 1989.

[4] Hallquist J.O., LS-Dyna3D, Theoretical Manual, Dynalis, 1995.

[5] CEA, Plexus, Theoretical Manual, CEA, 1997.

[6] Arnoult E., *Modélisation numérique et approche expérimentale du contact en dynamique* — Application au contact aubes/carter de turboréacteur, Thèse de doctorat, Université de Nantes, 2000.

[7] Arnoult E., Peseux B., Bonini J., *Recherche analytique du contact appliquée à un code éléments finis en dynamique*, in : Actes du 14ème Congrès Français de Mécanique AUM'99, Toulouse, CD-ROM 218, 1999, pp. 1-6.

[8] Hallquist J.O., Goudreau G.L., Benson D.J., Sliding interfaces with contact-impact in large-scale Lagrangian computations, Comput. Methods Appl. Mech. Engrg. 51 (1985) 107-137.

[9] Guilloteau I., Peseux B., Bonini J., *Simplified dynamic contact model — Application to the rotor/stator interaction*, in : Proceedings of the European Conference on Computational Mechanics ECCM'99, CD-ROM 592, 1999, pp. 1-13.

[10] Guilloteau I., Peseux B., Bonini J., *Modélisation du contact en implicite — Application à l'interaction rotor/stator*, Revue Européenne des Éléments Finis 9 (1/2/3) (2000) 331-348.

[11] Carpenter N.J., Taylor R.L., Katona M.G., Lagrange constraints for transient finite element surface contact, Int. J. Numer. Methods Engrg. 32 (1991) 103-128.

[12] Chaudhary A.B., Bathe K.J., *A solution method for static and dynamic analysis of three-dimensional contact problems with friction*, Comp. Struct. 24 (1986) 855-873.

[13] Guilloteau I., *Modélisation du contact en implicite — Application à l'interaction rotor/stator*, Thèse de doctorat, Université de Nantes, 1999.

[14] Géradin M., Rixen D., *Théorie des vibrations — Application à la dynamique des structures*, 2ème édition, Masson, Paris, 1996.

[15] Hughes T.J.R., Taylor R.L., Sackman J.L., Curnier A., Kanok-nukulchai W., *A finite element method for a class of contact-impact problems*, Comput. Methods Appl. Mech. Engrg. 8 (1976) 249-276.

[16] SAMCEF, MECANO/ROTORT Manual, Version 7.1, SAM-TECH S.A., 1997.