# Gestion du contact dans les simulations multi-domaines de crashs automobiles

Benjamin Bourel, Alain Combescure, Laurent Di Valentin

► To cite this version:

Benjamin Bourel, Alain Combescure, Laurent Di Valentin. Gestion du contact dans les simulations multi-domaines de crashs automobiles. 7e colloque national en calcul des structures, CSMA, May 2005, Giens, France. hal-01812892

HAL Id: hal-01812892

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

Submitted on 11 Jun 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
# Gestion du contact dans les simulations multi-domaines de crashs automobiles

Benjamin Bourel*,** — Alain Combescure* — Laurent Di Valentin**

*Laboratoire de Mécanique des Contacts et des Solides, INSA Lyon
20, avenue Albert Einstein, 69621 Villeurbanne
alain.combescure@insa-lyon.fr

**PSA Peugeot-Citroën
Route de Gizy, 78943 Vélizy-Villacoublay
{benjamin.bourel,laurent.divalentin}@mpsa.com

<footnote>RÉSUMÉ. Dans ce papier nous présentons une méthode de gestion du contact unilatéral, pour les calculs multi-domaines et multi-échelles temporelles en dynamique explicite. L'intérêt de la méthode, basée sur une modélisation du contact par pénalité et une représentation de type « maître-esclave », est de découpler le plus possible chaque sous-domaine : il peut ainsi évoluer avec son propre pas de temps de stabilité. L'algorithme proposé assure la continuité cinématique aux interfaces en échangeant un minimum d'informations entre les sous-domaines. Nous présentons un exemple montrant la validité et l'efficacité de la méthode par rapport à l'approche mono-domaine classique.</footnote>

<footnote>ABSTRACT. In this paper we present a method of management of the unilateral contact, for multi-domaines and multi-time-scales calculations in explicit dynamic. The interest of the method, based on a modeling of contact by penalty and a representation of the « master-slave » type, is to uncouple each domain as much as possible : it can thus evolve with its own time step of stability. The algorithm proposed ensures kinematic continuity throughout the interface by exchanging a minimum of information between the domains. We present an example showing the validity and the effectiveness of the method compared to the traditional mono-domain approach.</footnote>

<footnote>MOTS-CLÉS : Modélisation du contact, dynamique explicite, éléments finis, multi-domaines, multi-échelles.</footnote>

<footnote>KEYWORDS: Modeling of contact, explicit dynamic, finite elements, multi-domains, multi-scales.</footnote>
# 1. Introduction

En dynamique explicite, l'analyse des structures est souvent très coûteuse en temps de calcul à cause du grand nombre de pas de temps nécessaires à la stabilité des modèles. C'est pourquoi, il est crucial de trouver des algorithmes plus performants pour diminuer les temps de calcul. Tel est l'intérêt des méthodes de décomposition de domaine [CHA 00] qui consistent à diviser le modèle éléments finis en plusieurs sous-domaines ayant chacun leur propre schéma d'intégration. Cela permet d'appliquer à chaque sous-domaine la discrétisation la plus pertinente et par conséquent le pas de temps le plus économique. Ces travaux font suite à une première étude, traitant du couplage de sous-domaines incompatibles en espace et en temps [HER 02, FAU 03]. Le but est ici d'étendre la méthode mise en oeuvre par la gestion du contact unilatéral (avec ou sans frottement) entre les différents sous-domaines. Ce développement a une grande importance pour de nombreux cas pratiques tels que les simulations crashes automobiles où l'on ne peut à priori exclure des contacts de ce type.

# 2. Mise en place du problème

## 2.1. Modélisation du contact

L'interface de contact peut être définie par un ensemble de deux surfaces géométriques susceptibles de se rencontrer au cours d'une transformation. Le contact est ici modélisé par une technique classique de pénalité [BEL 00] et une représentation de type « maître-esclave » (cf. figure 1). La condition d'impénétrabilité des noeuds esclaves dans la surface maître se traduit alors par un effort normal pénalisant le long de la surface de contact.

## 2.2. Problème de référence

On considère un solide déformable occupant le domaine $\Omega$ de frontière $\Gamma$ susceptible de s'auto-impacter au cours de la déformation. Le problème d'équilibre discrétisé en espace s'écrit alors :

$$M^k \ddot{U} = F_{ext} - F_{int}(U) + F_{cont}(U, \dot{U}) \quad [1]$$

où $\ddot{U}$, $\dot{U}$, $U$ représentent les vecteurs accélération, vitesse et déplacement, $M$ représente la matrice masse, $F_{ext}$ les forces externes, $F_{int}$ les forces internes et $F_{cont}$ les forces de contact induites par la méthode de pénalité. Pour résoudre ce problème généralement non linéaire, on utilise ici un schéma explicite de type différence centrée :

$$U^{n+1} = U^n + \Delta t^n \dot{U}^{n+1/2} \quad \text{et} \quad \dot{U}^{n+1/2} = \dot{U}^{n-1/2} + \Delta t^{n-1/2} \ddot{U}^n \quad [2]$$

où $\dot{U}^{n+1/2}$ est la vitesse calculée au demi pas de temps et $\Delta t^{n-1/2}$ est le pas de temps moyen $(\Delta t^{n-1} + \Delta t^n)/2$. L'avantage d'un tel schéma est qu'à l'instant $t^{n+1}$ le second membre de l'équation d'équilibre [1] est connu explicitement. A chaque cycle,
la solution $\ddot{U}^{n+1}$ peut donc être calculée de manière directe. Il existe néanmoins un inconvénient majeur lié à la stabilité conditionnelle du schéma. En effet, le pas de temps utilisé au cours du calcul ne doit pas dépasser la valeur critique $\Delta t_{crit}$ imposée par la condition de Courant [BEL 00].

## 2.3. Décomposition de domaines

Il convient ici de distinguer « l'interface de contact » de « l'interface de liaison ». L'interface de liaison représente un lien permanent entre les deux sous-domaines alors que l'interface de contact n'agit que lorsque les deux sous-domaines s'impactent. Décomposons le domaine $\Omega$ en deux sous-domaines $\Omega^1$ et $\Omega^2$ et notons $\Gamma$ l'interface de liaison entre ces deux sous-domaines. Le problème discret [1] à résoudre peut alors désormais s'écrire :

$$ \forall k \in \{1, 2\} \qquad M^k \ddot{U}^k = F_{ext}^k - F_{int}^k + F_{cont}^k + F_{link}^k \qquad [3] $$

$$ \sum_{k=1,2} F_{link}^k = 0 \qquad [4] $$

Alors que le système [3] traduit l'équilibre de chaque sous-domaine $\Omega^k$, la relation [4] exprime l'équilibre de l'interface $\Gamma$. La résolution d'un tel problème peut être faite par une décomposition « free-link ». De manière générale le calcul multi-domaines se déroule alors en trois étapes :

$$ \forall k \in \{1, 2\} \qquad M^k \ddot{U}_{free}^k = F_{ext}^k - F_{int}^k + F_{cont}^k \qquad [5] $$

$$ \forall k, l \in \{1, 2\}, k \neq l \qquad M^k \ddot{U}_{link}^k = M^l \ddot{U}_{free}^l \qquad [6] $$

$$ \forall k \in \{1, 2\} \qquad \ddot{U}^k = \ddot{U}_{free}^k + \ddot{U}_{link}^k \qquad [7] $$

La première étape correspondant au problème free [5] est résolue par un code de dynamique explicite tel que le logiciel Radioss$^{①}$. Du point de vue multi-domaines, ces calculs constituent les programmes « clients ». La seconde étape correspondant au problème link [6] ainsi que la troisième correspondant au calcul des quantités globales [7] sont gérées par un programme externe dit programme « principal ».

# 3. Gestion du contact entre sous-domaines

## 3.1. Principe de la méthode

La méthode mise en oeuvre, peut être vue comme la juxtaposition d'une « interface de contact » avec une « interface de liaison » entre sous-domaines (cf. figure 1). On considère alors deux sous-domaines : le sous-domaine A sera choisi comme étant le sous-domaine maître et le sous-domaine B comme étant le sous-domaine esclave. Le principe consiste à copier une partie des noeuds du sous-domaine esclave B dans le sous-domaine A et de gérer le contact dans le sous-domaine A uniquement. L'interface de liaison entre les deux sous-domaines permet de lier rigidement les noeuds esclaves proprement dit et les noeuds esclaves copiés.
<figure><img src="image_2.png" /><figcaption>Figure 1. Principe de la juxtaposition « interface de contact - interface de liaison »</figcaption></figure>

## 3.2. Cas du calcul mono-échelle

En mono-échelle temporelle, tous les calculs clients évoluent avec le même pas de temps (correspondant au pas de temps minimum repéré sur l'ensemble des sous-domaines). Les instants d'équilibres de chaque sous-domaine étant synchronisés, la continuité cinématique à l'interface peut-être calculée exactement à chaque cycle. La condition de liaison imposée par la relation [6] revient à corriger les accélérations free des noeuds d'interface calculées par le programme client en y ajoutant le terme link. En notant $\ddot{U}_A(t^n)$ (resp. $\ddot{U}_B(t^n)$) l'accélération free des noeuds d'interface envoyée par le sous-domaine A (resp. B) au programme principal à l'instant $t^n$, et $\ddot{U}_A^*(t^n)$ (resp. $\ddot{U}_B^*(t^n)$) l'accélération corrigée qui lui est retournée, on peut écrire d'après [7] :

$$ \ddot{U}_A^*(t^n) = \ddot{U}_B^*(t^n) = \ddot{U}_A(t^n) + \ddot{U}_B(t^n) \quad [8] $$

## 3.3. Cas du calcul multi-échelles

Notre méthode multi-domaines prend néanmoins tout son intérêt avec l'utilisation du multi-échelles. On autorise alors chaque sous-domaine à évoluer avec son propre pas de temps de stabilité. Comme précédemment, à chaque cycle de calcul les accélérations renvoyées à chacun des sous-domaines sont corrigées en y ajoutant les accélérations free du sous-domaine voisin. Ainsi pour le sous-domaine A on a :

$$ \ddot{U}_A^*(t_A^n) = \ddot{U}_A(t_A^n) + \ddot{U}_B(t_A^n) \quad [9] $$

Contrairement au cas mono-échelle, les sous-domaines ne sont plus à l'équilibre aux mêmes instants ($t_A^n \neq t_B^n$). Par exemple, à l'instant $t_A^n$, la correction d'accélération n'est plus triviale puisque le terme $\ddot{U}_B(t_A^n)$ est à priori inconnu. On procède alors à une interpolation linéaire des accélérations entre les instants $t_B^n$ et $t_B^{n+1}$ [HER 02] :

$$ \ddot{U}_B(t_A^n) = \ddot{U}_B(t_B^n) \left[ \frac{t_B^{n+1} - t_A^n}{t_B^{n+1} - t_B^n} \right] + \ddot{U}_B(t_B^{n+1}) \left[ \frac{t_A^n - t_B^n}{t_B^{n+1} - t_B^n} \right] \quad [10] $$

A chaque cycle, la continuité d'accélération entre deux noeuds d'interface en vis à vis n'est plus imposée exactement mais à partir de grandeurs interpolées (cf. figure 2). L'erreur commise par l'interpolation se répercute et se cumule logiquement sur les grandeurs intégrées que sont les vitesses et les déplacements. Afin d'imposer un déplacement des noeuds copiés de A identique à celui de leur homologue coté esclave B,
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Evolution des calculs en multi-échelles</figcaption></figure>

on apporte une correction $\Delta \ddot{U}_A$ aux accélérations $\ddot{U}_A^*$ renvoyées à ces noeuds copiés. En notant $\Delta U_A$ le décallage ($U_B - U_A$) sur les noeuds d'interface, on a :

$$ \ddot{U}_A^*(t_A^n) = \ddot{U}_A^*(t_A^n) + \Delta \ddot{U}_A(t_A^n) \quad \text{avec} \quad \Delta \ddot{U}_A(t_A^n) = \frac{\Delta U_A(t_A^n)}{\Delta t_A^n \Delta t_A^{n-1/2}} \quad [11] $$

# 4. Exemple numérique

## 4.1. Description du cas test

On s'intéresse ici à la simulation d'un choc tête à 35km/h sur un véhicule projet PSA Peugeot-Citroën®. La tête est modélisée par une sphère venant impacter le haut du capot, au niveau du bras d'essuie vitre (cf. figure 3a). Dans ce modèle, le système d'essuie vitre maillé très finement dispose d'un pas de temps pénalisant fortement le reste du véhicule (cf. tableau 1). Nous cherchons donc à isoler par notre méthode de décomposition, ce sous-domaine essuie-vitre (PART1) du reste de la structure (PART2). Les deux sous-domaines sont ici liés par 5 interfaces dont 3 liaisons permanentes et 2 liaisons de contacts. Les calculs ont été réalisés sur une machine mono-processeur cadencée à 3.75GHz.

## 4.2. Validation des résultats

L'intérêt de telles simulations est d'étudier le comportement de la structure ainsi que certains critères bio-mécaniques. Un des critères les plus importants concerne la décélération calculée au centre de l'impacteur tête. Les courbes d'accélération (cf. figure 3b) permettent de calculer le critère HIC (Head Injury Criterion) très utilisé dans

<table><caption>**Tableau 1.** Comparaison des modèles mono-domaine et multi-domaines</caption><thead><tr><th rowspan="2"></th><th>Mono-domaine</th><th colspan="4">Multi-domaines</th></tr><tr><th>Modèle global</th><th>PART1</th><th>PART2</th><th>Rad2rad</th><th>Total</th></tr></thead><tbody><tr><th>Nombre d'éléments</th><td>317451</td><td>13803</td><td>303648</td><td>-</td><td>317451</td></tr><tr><th>Δt<sub>crit</sub></th><td>0,339.10's</td><td>0,339.10's</td><td>0,834.10's</td><td>-</td><td>-</td></tr><tr><th>Temps CPU</th><td>391390 s</td><td>19964 s</td><td>148567 s</td><td>9257 s</td><td>177788 s</td></tr></tbody></table>
<figure><img src="image_4.png" /><figcaption>**Figure 3.** (a) Modèle choc piéton - (b) Sous-domaine essuie-vitre maillé finement</figcaption></figure>

le domaine de la biomécanique. L’écart de ce critère calculé à partir des deux méthodes est inférieur à 1%. En ce qui concerne les temps CPU, ceux-ci sont résumés dans le tableau 1. On remarque alors que la méthode mis en oeuvre apporte un gain substantiel en temps de calcul de 2,2 par rapport au modèle mono-domaine de référence. Ce gain prend en compte le temps passé dans le programme principal (Rad2rad©) ce qui représente près de 5% du temps de calcul total.

**5. Conclusion**

Cette méthode s'avère être d'un intérêt très pratique dans les domaines d'appli-
cations tel que celui du crash automobile, où le phénomène de contact est inévitable.
Les résultats obtenus ont non seulement révélé la très bonne qualité de la méthode
proposée, mais ont également montré l'efficacité de l'approche multi-domaines sur
les temps de calcul. Cet intérêt sera bien évidemment accru par l'utilisation de calcu-
lateurs parallèles qui permettront le traitement simultané de chaque sous-domaine par
un ou plusieurs processeurs. Toutefois, le gain en temps de calcul dépend fortement de
la discrétisation des sous-domaines en contact et s'avère être d'autant plus important
que la différence de pas de temps entre les deux sous-domaines est grande. Le choix
de la décomposition est donc primordial vis à vis de l'efficacité de la méthode.

**6. Bibliographie**

[BEL 00] BELYTSCHKO T., LIU W., MORAN B., *Nonlinear Finite Elements for Continua and Structures*, Wiley, 2000.

[CHA 00] CHAMPANEY L., DUREISSEIX D., « Calcul de structures et parallélisme : un bilan et quelques développements récents. Mécanique et Industries », *Mécanique et Industries*, vol. 1, n° 1, 2000, p. 43-46.

[FAU 03] FAUCHER V., COMBESCURE A., « A time an space mortar method for coupling li-
near modal subdomains and non-linear subdomains in explicit structural dynamics », *Com-
put. Methods Appl. Mech. Engrg*, vol. 192, 2003, p. 509-533.

[HER 02] HERRY B., DI VALENTIN L., COMBESCURE A., « An approach to the connection between subdomains with non-matching meshes for transient mechanical analysis », *Int. J. Numer. Meth. Engng*, vol. 55, 2002, p. 973-1003.