# Gestion des contacts 3D pour une modélisation Éléments Discrets

Issa Sanni, Jérôme Fortin, Patrice Coorevits

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

To cite this version:

Issa Sanni, Jérôme Fortin, Patrice Coorevits. Gestion des contacts 3D pour une modélisation Éléments Discrets. 7e colloque national en calcul des structures, CSMA, May 2005, Giens, France. hal-01813008

HAL Id: hal-01813008

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

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

Distributed under a Creative Commons CC0 1.0 - Universal - International License
# Gestion des contacts 3D pour une modélisation Eléments Discrets

## Issa Sanni - Jérôme Fortin - Patrice Coorevits

Laboratoire de Technologies Innovantes (EA 3899)
Université de Picardie Jules Verne
48 rue d'Ostende,
Saint-Quentin 02 100, France
jerome.fortin@insset.u-picardie.fr

**RÉSUMÉ.** L'objectif de ce papier est de proposer une stratégie pour optimiser la gestion des contacts 3D dans un code de calcul par Eléments Discrets (ED). La modélisation ED utilisée s'appuie sur la méthode de Dynamique des Contacts (DC) et sur la méthode du bipotentiel. On aboutit à un algorithme comportant une phase globale de résolution de l'équation de la dynamique et une phase locale de résolution du schéma de prédiction-correction qui fournit une nouvelle valeur de l'impulsion. On constate alors que la gestion des contacts consomme une grande partie du temps machine (90 % du temps CPU). Afin de restreindre le nombre de cas à examiner, on propose dans ce papier un algorithme performant de détection des contacts 3D. En terme de temps CPU, la stratégie mise en place pour gérer les contacts 3D permet de diminuer considérablement le temps machine.

**ABSTRACT.** The aim of this paper is to propose a strategy to optimize the 3D contacts management in a computer code by Discrete Elements Method (DEM). The DEM modelling used is based on the Contact Dynamics method and on the bipotential method. This leads to an algorithm comprising a global stage of resolution of the dynamic equations and a local stage of resolution of the prediction- correction equations which provides a new value of the impulse. We noted whereas that the contacts management consumes most of the machine time (90 percent of CPU time). In order to restrict the number of cases to be examined, we propose in this paper a reliable algorithm of 3D contacts detection. In term of CPU time, this strategy enable to decrease considerably the machine time (CPU time).

**MOTS-CLÉS :** contact frottant non linéaire ; milieu granulaire ; formulation variationnelle

**KEYWORDS:** non-smooth frictional contact; Granular media; variational formulation
# 1. Introduction

Le matériau granulaire concerne des secteurs d'activité aussi divers que le Génie Civil, les industries minières, agro-alimentaires, pharmaceutiques, etc... À la différence avec d'autres mélanges fluide-solide, l'interaction entre les particules joue un rôle très important dans la manipulation de ces systèmes complexes. L'objet des travaux sur ce type de matériau est donc de mieux comprendre le comportement d'ensemble d'un grand nombre de grains ou de corps rigides à partir de leur cinématique et des efforts qu'ils transmettent. La diversité des modélisations Éléments Discrets (ED) provient pour une bonne part de la difficulté à représenter correctement l'interaction entre les corps.

L'objectif de ce papier est de proposer une stratégie pour optimiser la gestion des contacts 3D en fonction des non-linéarités rencontrées (ici la loi de contact unilatéral avec frottement de Coulomb). La modélisation ED utilisée s'appuie sur la méthode de Dynamique des Contacts (DC) développée par Jean [JEA 99] et Moreau [MOR 94] et sur la méthode du bipotentiel développée par de Saxcé [SAX 98] et appliquée aux systèmes granulaires 2D par Fortin [FOR 00]. Les systèmes mécaniques relevant de la DC exigent un traitement implicite de l'équation de la dynamique et des relations interparticulaires. On aboutit alors à un algorithme comportant une phase globale de résolution de l'équation de la dynamique qui calcule une nouvelle approximation de la vitesse et une phase locale de résolution du schéma de prédiction-correction qui fournit une nouvelle valeur de l'impulsion. On constate alors que la gestion des contacts consomme une grande partie du temps machine (90 % du temps CPU). Pour minimiser cette consommation, il nous faut satisfaire à deux conditions : trier rapidement les candidats potentiels aux contacts et calculer une solution "possible" des impulsions avec un nombre limité d'itérations locales.

# 2. Détection des contacts potentiels en 3D

Concernant la sélection des candidats aux contacts [FOR 04], classiquement, pour chaque candidat $\omega_i$ et pour tous les antagonistes $\omega_j$ avec $i > j$, on doit effectuer $\frac{N(N+1)}{2}$ vérifications ce qui devient très rapidement coûteux dès que N devient élevé. Afin de restreindre le nombre de contacts potentiels à examiner, on propose un algorithme performant de détection des contacts. Celui-ci se décompose en quatre étapes : maillage du système à étudier, identification des tables de connectivités, appartenance des particules au patch et vérification du critère de sélection.

L'espace $\sum$ étudié est subdivisé en plusieurs sous-espaces (ou mailles), figure1 $\Omega_i, i = 1, ..., p, p$ représentant le nombre maximun de maille qu'il est possible de construire dans le système étudié,

$$p = \frac{(\max_x - \min_x)(\max_y - \min_y)(\max_z - \min_z)}{d_x d_y d_z}, \quad [1]$$
avec $d_j = k_j \emptyset_{max}$, où $\emptyset_{max}$ représente le diamètre maximal de la particule à simuler et $k_j$ un paramètre qui permet d'ajuster la taille de la maille en fonction du nombre de particules désiré par maille et par axe.

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Exemple de maillage de l'espace de simulation.</figcaption></figure>

Maintenant pour chaque maille $\Omega_i$, nous allons calculer sa position puis identifier ses voisins au moyen de celle-ci pour ainsi constituer sa table de connectivité. La position d'une maille $\Omega_i$ est identifier par $ix$, $iy$ et $iz$. Pour les calculer on utilise la procédure suivante

$$k = 1, nm_z \quad \left[ \begin{array}{l} l = (k-1)nm_{xy} + 1 \\ i = l, knm_{xy} \\ \left[ \begin{array}{l} iy(i) = I_p \left( \frac{(i-(k-1)nm_{xy})-1}{nm_x} \right) + 1 \\ ix(i) = (i - (k-1)nm_{xy}) - (iy(i) - 1)nm_x \\ iz(i) = k \end{array} \right. \end{array} \right. \quad [2]$$

où $I_p$ représente la partie entière alternative, $nm_z$ le nombre de mailles suivant l'axe z et $nm_{xy} = nm_x \times nm_y$. Une fois les mailles positionnées, il nous reste à identifier pour chaque maille ses voisins et établir les tables de connectivité. L'objectif étant de limiter notre champ d'exploration des contacts potentiels à ces tables. En fonction de sa position, la maille peut avoir entre 1 et 14 voisins. Ce calcul inclut la maille elle même. Le processus utilisé ici pour identifier le nombre de voisins et le nom des mailles voisines utilise les dimensions du maillage $nm_x$, $nm_y$ et $nm_z$. A titre d'exemple, on montre à la figure2 une maille et ses voisins.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Table de connectivité : la maille étudiée est $\Omega_5$. En gris foncé sont représentés ses voisins, ici au nombre de 14 le maximum</figcaption></figure>

Il nous reste alors à identifier l'appartenance, pour chaque particule, à tel ou tel patch représentant l'ensemble des particules d'une même table de connectivité. Pour ce faire on utilise les coordonnées des particules, $qc$.

$$i = I_p \left( \frac{qc_x - min_x}{d_x} \right) + I_p \left( \frac{qc_y - min_y}{d_y} \right) nm_x + I_p \left( \frac{qc_z - min_z}{d_z} \right) nm_{xy} + 1 \quad [3]$$

Enfin, dans chaque patch de particules, chaque couple de candidats au contact potentiel est testé sur le critère de sélection utilisé. La condition de contact s'appuie sur l'analyse de l'interstice $u_n$.

$$u_n = \|qc_i - qc_j\| - (a_i + a_j) \le \beta \quad [4]$$

où $qc_i$, $qc_j$, $a_i$, $a_j$ sont respectivement les coordonnées des centres de gravité et les rayons des corps $\Omega_i$ et $\Omega_j$. $\beta$ étant la distance à considérer, dans notre cas elle est de l'ordre de $10^{-8}$. Ce qui nous donnera alors une liste de candidats au contact effectif sur lesquels on calculera les impulsions $s_c$.

### 3. Calcul des impulsions et qualité de la solution

La loi de contact utilisée ici est la loi de contact unilatéral avec frottement de Coulomb. On a montré dans [FOR 00] qu'elle pouvait se mettre sous la forme d'un seul schéma de prédiction-correction :

$$ \begin{aligned} \text{prédiction :} \quad & \tau^{i+1} = s^i - \rho[\dot{u}_T^i + (\dot{u}_n^i + \mu|| - \dot{u}_T^i||)n], \\ \text{correction :} \quad & s^{i+1} = \text{proj}(\tau^{i+1}, K_\mu). \end{aligned} \quad [5] $$

où $K_\mu$ représente le cône de Coulomb, $s$ l'impulsion de contact, $\dot{u}_n$ la vitesse normale, $\dot{u}_T$ la vitesse de glissement et $\rho$ un paramètre numérique qui est fonction du pas de temps, du coefficient de frottement $\mu$, des coefficients de restitution normale et tangentielle et de la masse réduite. La projection de $\tau$ sur le cône de Coulomb conduit à la possibilité de trois statuts :
<table><tr><td>si $\tau^{i+1} \in K_{\mu}^*$</td><td>non contact,</td><td></td></tr><tr><td>si $\tau^{i+1} \in K_{\mu}$</td><td>contact avec frottement,</td><td>[6]</td></tr><tr><td>si $\tau^{i+1} \in R^3 - (K_{\mu} \cup K_{\mu}^*)$</td><td>contact avec glissement.</td><td></td></tr></table>

Pour la qualité de la solution, on utilise comme critère d’arrêt l’indicateur d’erreur en loi de comportement [LAD 96] au sens strict basé sur la “violation” du bipotentiel de contact [FOR 00].

$$ \epsilon = \sum b_c(-\dot{\mathbf{u}}^{i+1}, \mathbf{s}^{i+1}) + \dot{\mathbf{u}}^{i+1} \cdot \mathbf{s}^{i+1} \quad [7] $$

### 4. Application

A titre d’application, nous proposons la vidange d’un système granulaire constitué de 1000 particules sphériques soumis à la seule pesanteur et aux efforts interparticulaires, figure 3. Le frottement est de l’ordre de 0.2. Les coefficients de restitution normales et tangentielles sont égales respectivement à 0. et 1. Le pas de temps est égal à $10^{-4}$ seconde. A l’initialisation, les particules de taille homogène sont positionnées de la même façon qu’un assemblage rectangulaire.

<figure><img src="image_5.png" /><figcaption>**Figure 3.** simulation 3D d'un matériau granulaire ensilé.</figcaption></figure>

Pour information, il faut 15 heures pour simuler la vidange sur une station de travail HP B2000 avec un processeur de 400 MHz et 256 Mb de RAM. En terme de temps CPU, la détection des contacts potentiels prend beaucoup plus de temps que le calcul des impulsions, voir tableau 4.
<table><tr><td>temps de la simulation (s)</td><td>2.10<sup>-4</sup></td><td>2,87</td><td>4,5</td><td>6,16</td></tr><tr><td>sélection des contacts potentiels (CPU)</td><td>0,94</td><td>0,82</td><td>0,82</td><td>0,84</td></tr><tr><td>calcul des impulsions (CPU)</td><td>0,48</td><td>0,1</td><td>0,16</td><td>0,3</td></tr><tr><td>Nombre de contacts</td><td>2900</td><td>667</td><td>860</td><td>1194</td></tr><tr><td>Nombre d'itérations locales</td><td>4</td><td>3</td><td>4</td><td>5</td></tr></table>

En moyenne, sur ce type de simulation, il nous faut moins de 5 itérations locales pour converger vers une solution acceptable. Il est à noter que dans la détection des contacts est inclut, pour les simples contacts effectifs, le calcul des matrices de passage entre les variables locales et globales [FOR 00].

## 5. Conclusion

En terme de temps CPU, la stratégie mise en place pour gérer les contacts 3D (triage des candidats potentiels et le calcul des impulsions) permet de diminuer considérablement le temps machine. En terme de qualité, l'utilisation d'une seule formulation variationnelle nous permet de converger très rapidement (entre 2 et 5 itérations) avec des résultats de simulations très correcte. L'objectif à terme est d'utiliser un maillage adaptatif qui doit nous permettre de diminuer le temps CPU dans cette procédure.

## 6. Bibliographie

[FOR 00] FORTIN J., « Simulation numérique de la dynamique des systèmes multicorps appliquée aux milieux granulaires », PhD thesis, Université de Lille 1, january 2000.

[FOR 03] FORTIN J., COOREVITS P., Stratégie de calcul du pas optimal pour une modélisation Éléments Discrets, 20-23 mai 2003.

[FOR 04] FORTIN J., COOREVITS P., « Selecting contact particles in dynamics granular mechanics systems », Journal of Computational and Applied Mathematics, vol. 168, 2004, p. 207-213.

[JEA 99] JEAN M., « The non-smooth contact dynamics method », Computer Methods in Applied Mechanics and Engineering, vol. 177, 1999, p. 235-257.

[LAD 96] LADEVÉZE P., Mécanique non linéaire des structures : nouvelle approche et méthodes de calcul non incrémentales, Hermés, 1996.

[MOR 94] MOREAU J.-J., « Some numerical methods in multibody dynamics : application to granular materials », Eur.J. Mech., A/Solids, vol. 13, 1994, p. 93-114.

[SAX 98] DE SAXCÉ G., FENG Z.-Q., « The bipotential method : a constructive approach to design the complete contact law with friction and improved numerical algorithms », Math. Comput. Modelling, vol. 28, 1998, p. 225-245.