# Simulation d'un problème de contact par des méthodes de réduction de modèle pour une application temps réel

Bastien Lefèvre, Frédéric Druesne, Jean-Luc Dulong, Pierre Villon

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

## To cite this version:

Bastien Lefèvre, Frédéric Druesne, Jean-Luc Dulong, Pierre Villon. Simulation d'un problème de contact par des méthodes de réduction de modèle pour une application temps réel. 9e Colloque national en calcul des structures, CSMA, May 2009, Giens, France. hal-01422169

HAL Id: hal-01422169

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

Submitted on 23 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_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Simulation d'un problème de contact par des méthodes de réduction de modèle pour une application temps réel

B. Lefèvre, F. Druesne, J.-L. Dulong, P. Villon

Laboratoire Roberval - Université de Technologie de Compiègne - UTC-CNRS, UMR 6253
Centre de Recherches de Royallieu
BP 20529 - 60205 COMPIEGNE cedex FRANCE
{bastien.lefevre, frederic.druesne, jean-luc.dulong, pierre.villon}@utc.fr

**Résumé** – L'interaction en temps-réel entre un opérateur et un prototype virtuel, permettant de simuler par exemple un assemblage, est une aide précieuse à la conception. Nous proposons ici une approche en deux phases : apprentissage et temps-réel. L'objectif de la première phase est d'obtenir le plus rapidement possible une surface de réponse réduite. Des méthodes de réduction de modèle sont utilisées avec une formulation en minimisation du résidu.

**Mots clés** – contact, minimisation du résidu, réduction de modèle.

## 1 Introduction

La possibilité de réaliser une simulation temps-réel d'un assemblage présente un intérêt certain lors de la phase de conception. Ce type de simulateurs existe mais se limite à des hypothèses d'élasticité linéaire [2]. La méthode des éléments finis est appropriée pour tenir compte des non linéarités liées à ce genre de problème : géométrique, matérielle et de conditions aux limites (associés au contact). Or le calcul en temps-réel s'avère impossible dès lors que la taille du problème devient importante.

Une approche en deux phases a été développée : une phase d'apprentissage et une phase temps-réel [5]. Au cours de la première, une campagne de calcul préliminaire est réalisée [1]. Cette campagne de calcul représente une grille de différentes manipulations possibles du modèle, également appelées cas de charge. L'objectif de cette première phase est d'obtenir une surface de réponse réduite qui sera utilisée en temps-réel au cours de la seconde phase par interpolation. Nous nous focalisons ici sur la première phase. Afin d'obtenir une surface de réponse réduite, nous présentons trois niveaux de réductions de modèle : la méthode adaptative [4], la Transformée de Karhunen-Loeve [3] [6], et l'hyperreduction [8]. Ces techniques de réduction de modèle sont ensuite employées pour étudier l'écrasement d'un cylindre en caoutchouc sur un plan rigide.

## 2 Méthode adaptative

### 2.1 Présentation du problème

Le problème à résoudre conduit à l'équation (1), avec $r(u)$ le résidu numérique :
$$K(u)u - f(u) = r(u) \qquad (1)$$

Afin de résoudre (1) par l'intermédiaire d'un schéma de Newton-Raphson, le résidu est linéarisé par un développement de Taylor du premier ordre, à l'itération $i$ :

$$r(u^i) = r(u^{i-1}) + K_T(u^{i-1})du \qquad (2)$$

L'équation (2) sera noté $R(du)$. L'équation (1) présente $n$ degrés de liberté, potentiellement grand. Dans le but de réduire le coût de résolution de ce problème, une méthode adaptative est utilisée.

## 2.2 Formulation par minimisation du résidu

Le déplacement $u_k$, correspondant au $k$<sup>ième</sup> cas de charge, est projeté sur une base $\Phi=\{\phi_1, ..., \phi_m\}$ de taille réduite $m$ avec $m \ll n$. Cette base est associée à un vecteur de coefficients $a$ :

$$u_k = \Phi a_k \qquad (3)$$

L'équation (2) s'écrit alors :

$$R(da) = r(a_k^{i-1}) + K_T(a_k^{i-1})\Phi da \qquad (4)$$

Le problème présente $n$ équations pour $m$ inconnus, il est donc surdéterminé. Afin de résoudre ce système, l'équation (4) est minimisée au sens d'une norme $Q$ :

$$J_Q(da) = \frac{1}{2} \|R\|_Q^2 = \frac{1}{2} R^T(da)QR(da) \qquad (5)$$

Le choix retenu pour la norme est $Q = I$. La dérivée de (5) dans la direction $da$, pour tout $\delta da$, est :

$$J'.\langle \delta da \rangle = R^T K_T \Phi.\langle \delta da \rangle \qquad (6)$$

L'équation (6) est annulée :

$$[\Phi^T K_T^T K_T \Phi] da = -\Phi^T K_T^T r \qquad (7)$$

La matrice de rigidité tangente réduite (8) est de dimension ($m \times m$) et le résidu réduit (9) de dimension ($m \times 1$).

$$\hat{K}_T = \Phi^T K_T^T K_T \Phi \qquad (8)$$

$$\hat{r} = \Phi^T K_T^T r \qquad (9)$$

Le nombre $m$ de modes est très petit par rapport au nombre de ddl $n$, l'inversion de la matrice de rigidité tangente réduite pleine (8) est donc moins coûteux que l'inversion de la matrice de rigidité tangente non réduite bande.

Durant les itérations de Newton-Raphson permettant de résoudre le système ayant pour inconnu le vecteur de coefficients $a$, la base $\Phi$ est supposée connue. Cette base est initialisée par un premier cas de charge $u_l$ obtenue par des itérations de Newton-Raphson sur $u$, d'où les premiers mode (10) et coefficient (11) :

$$\phi_l = \frac{u_l}{\|u_l\|} \qquad (10)$$

$$a_l = \|u_l\| \qquad (11)$$

Au cours de la campagne de calculs, correspondant à aux différents cas de charge, la base peut s'avérer trop pauvre pour décrire correctement les nouveaux cas de charge et donc induit un déplacement calculé incorrect. Il est donc nécessaire d'enrichir la base avec de nouveaux modes. Plusieurs choix sont envisageables pour calculer ces nouveaux modes, l'enrichissement par le résidu [7] par exemple. Le choix retenu, permettant une bonne convergence de l'algorithme est d'ajouter un nouveau mode calculé par :
$$J'.\langle \delta da \rangle = R^T PPK_T \Phi.\langle \delta da \rangle \quad (17)$$

L'équation (17) est annulée :

$$[\Phi^T K_T^T P K_T \Phi] da = -\Phi^T K_T^T P r \quad (18)$$

La matrice de rigidité tangente hyperréduite est symétrique et pleine.

# 5 Modèle et résultats

## 5.1 Modèle

Le problème étudié est l'écrasement d'un cylindre en caoutchouc sur un plan indéformable. Ce modèle est maillé en éléments de type Q4, avec $n = 10000$ degrés de liberté. Les déplacements sont imposés sur la ligne supérieure. Pour la gestion du contact, l'algorithme du Lagrangien augmenté est utilisé. Les cas de charge de la campagne de calculs correspondent à des niveaux successifs d'écrasement. Ce modèle présente des non linéarités matérielle (matériau hyper élastique), géométrique et de conditions aux limites.

## 5.2 Résultats et discussion

L'objectif étant d'obtenir une surface de réponse réduite le plus rapidement possible, il est donc important d'observer l'évolution de la taille de la base ainsi que le comportement en temps de nos méthodes.

A travers ce papier, trois niveaux de réduction de modèle sont étudiés : la méthode adaptative, la méthode adaptative combinée à la TKL et pour finir l'hyperreduction combinée à la TKL.

La solution de référence est obtenue par un schéma de résolution Newton-Raphson classique sur le déplacement $u$. Il est important de noter que, quelque soit la méthode de réduction utilisée, l'erreur relative par rapport à la solution de référence est toujours inférieure à $10^{-4}$.

<figure><img src="image_3.png" /><figcaption>Figure 1 – taille de la base au cours de la campagne</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 2 – enrichissements cumulés au cours de la campagne</figcaption></figure>

La méthode adaptative a une évolution régulière en terme d'ajout de nouveaux modes, jusqu'à atteindre 101 à la fin de la campagne (Fig. 1). L'utilisation de la TKL pour les deux autres méthodes limite la taille de la base à 35 modes, correspondant à la valeur $m_{limite}$ fixée. Après une TKL, le déplacement obtenu est approché, on observe ainsi une augmentation locale du nombre d'enrichissements (Fig. 2). Puis le taux d'enrichissement redevient identique à celui d'avant la TKL.
<figure><img src="image_5.png" /><figcaption>Figure 3 – temps CPU des différentes méthodes au cours de la campagne</figcaption></figure>

Le second aspect important à étudier est le temps CPU nécessaire pour réaliser la campagne de calculs. Par rapport à la méthode classique de référence, on observe un gain de 55% pour la méthode adaptative (Fig. 3). Le gain lorsque cette méthode est combinée à la TKL est minime, en effet on observe un gain de 57% par rapport à la méthode classique. Ce faible écart est du au faible nombre de modes nécessaires à la méthode adaptative pour finaliser la campagne. Concernant l'hyperreduction combinée à la TKL, le gain est beaucoup plus important, de l'ordre de 76% par rapport à la méthode classique. Ce gain s'explique par le fait que le coût numérique de (18) est très faible en raison du nombre de ddl retenus.

<figure><img src="image_6.png" /><figcaption>Figure 4 – Maillages utilisés par l’hyperreduction</figcaption></figure>

Il est intéressant de noter que le maillage initial (Fig. 4.a) retenu correspond aux éléments associés aux ddl à conditions aux limites imposées. Au cours de la campagne, des ddl sont ajoutés, ces ddl sont localisés sur la zone de contact comme le montre le maillage final (Fig. 4.b). Pour conclure, on peut constater que le nombre ddl nécessaire pour mener à bien la campagne est $q = 523$, ce nombre étant très petit par rapport aux $n = 10000$ ddl du modèle global.

# 6 Conclusion

Les méthodes de réduction de modèle que nous avons utilisées, s’appuyant sur une formulation en minimisation du résidu, permettent d’obtenir une surface de réponse réduite pour un problème mettant en jeu des non linéarités géométrique, matérielle et de conditions aux limites. En effet, la quantité de données stockées à l’aide notre méthodologie s’avère considérablement réduite par rapport à une surface de réponse obtenue en utilisant un schéma
classique Newton-Raphson. De plus, notre méthodologie permet également de réaliser la campagne de calcul préliminaire de manière significativement plus rapide pour les deux méthodes de réduction de modèle utilisées.

La prochaine étape de nos travaux consistera à optimiser une campagne de calculs 2D. Pour ce faire, un modèle 3D tel que la mise en place d'une durit en caoutchouc sur un embout rigide sera utilisé.

Références

[1] F. Druesne, J.L. Dulong, P. Villon. Model reduction applied to simulation of mechanical behaviour for flexible parts. *The 8th International Conference on Computational Structures Technology*, 2006.

[2] C. Duriez. *Contact frottant entre objets déformables dans des simulations temps-réel avec retour haptique*, Thèse de Doctorat de l'Université d'Evry, 2004.

[3] K. Karhunen. *Uber lineare methoden für wahrscheinigkeitsrechnung*, Annales of Academic Science Fennicae Series A1 Mathematical Physics, Vol. 37, page 3 – page 79, 1943.

[4] P. Krysl, S. Lall, J.E. Marsden. Dimensional Model Reduction in Non-linear Finite Element Dynamics of Solids and Structures. *International Journal for Numerical Methods in Engineering*, 2000.

[5] B. Lefevre, F. Druesne, J.L. Dulong, P. Villon. Simulation of a mechanical assembly using model reduction. *The 8th. World Congress on Computational Mechanics (WCCM8)*, 2008.

[6] Loève M. *Probability Theory*, Van Nostrand, New Jersey, 1955.

[7] F. Risler et C. Rey. Iterative accelerating algorithms with Krylov subspaces for the solution to large-scale nonlinear problems. *Numerical Algorithms*, Vol. 23, page 1 – page 30, 2000.

[8] D. Ryckelynck. A priori hyperreduction method: an adaptive approach. *Journal of computational physics*, Vol. 202, page 346 – page 366, 2005.