# Simulation DEM des matériaux élasto-plastiques endommageables

Iñigo Terreros, Ivan Iordanoff, Jean-Luc Charles

► To cite this version:

Iñigo Terreros, Ivan Iordanoff, Jean-Luc Charles. Simulation DEM des matériaux élasto-plastiques endommageables. 11e colloque national en calcul des structures, CSMA, May 2013, Giens, France. hal-01717844

HAL Id: hal-01717844

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

Submitted on 26 Feb 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
# Simulation DEM des matériaux élasto-plastiques endommageables

Iñigo TERREROS<fnref n="1" />, Ivan IORDANOFF<fnref n="2" />, Jean-Luc CHARLES<fnref n="3" />

**Résumé** — De nos jours, la majorité des procédés de fabrication sont simulés à l'aide des méthodes basées sur la mécanique des milieux continus comme, par exemple, la méthode FEM. Ces méthodes de simulation trouvent souvent des difficultés pour simuler certains phénomènes complexes comme le frottement, l'écoulement de matière... typiquement liés aux procédés de mise en forme.

La méthode des éléments discrets (DEM) a désormais montré son adéquation à certains de ces problèmes, notamment ceux liés à la fracture du matériau et au contact sec en présence d'un troisième corps. Néanmoins, elle doit être adaptée pour simuler le comportement mécanique d'un matériau continu élasto-plastique. Ce travail décrit une méthode pour simuler des matériaux élasto-plastiques avec la DEM. Le matériau est modélisé par des sphères de différentes tailles afin d'obtenir un domaine compact et iso-trope. La loi d'interaction entre les éléments discrets est basée sur un potentiel de Lennard – Jones. Un algorithme de détection de surfaces est implémenté afin de suivre la propagation des fissures créées lors des chargements en traction. Finalement, l'influence des différents paramètres de la loi d'interaction a été étudiée.

**Mots clés** — DEM mécanique, élasto-plasticité, matériaux continus, propagation des fissures

## 1 Introduction

Comme mentionné dans le résumé précédent, la majorité des procédés de fabrication sont actuellement simulés avec des approches de type mécanique des milieux continus. Ces simulations doivent faire face à des problèmes liés à leur approche continue. La méthode des éléments discrets (DEM) [1] a désormais montré son adéquation à certains de ces problèmes, notamment à ceux liés à la fracture [2] et au contact sec en présence d'un troisième corps [3]. Néanmoins, elle doit être adaptée pour simuler le comportement mécanique du matériau continu élasto-plastique.

L'objectif à terme de ce travail est d'obtenir une méthode de simulation des matériaux élasto-plastiques avec la DEM. Cette méthode peut en plus être couplée avec des méthodes continues comme la FEM ou la NEM [4] dans les zones où ces approches sont plus rapides en temps de calcul (figure 1).

<figure><img src="image_2.png" /><figcaption>Fig. 1 – Exemple simplifié d'un procédé de fraisage. La zone en rouge et vert autour du contact outil-matière est modélisée avec la DEM et la zone hachurée avec une méthode continue.</figcaption></figure>

Plusieurs travaux ont été réalisés afin de simuler des matériaux continus avec la méthode DEM mais soit :

- Ils étaient peu adaptés aux problèmes d'écoulement [6].

- Ils étaient plus orientés vers des matériaux fragiles [2], [7], [8].

- Ils étaient très coûteux en temps de calcul et, donc, peu adaptés à la simulation des procédés de fabrication [5].

<footnote n="1">¹ I2M UMR 5295 CNRS, MPI, inigo.terreros@u-bordeaux1.fr</footnote>

<footnote n="2">² I2M UMR 5295 CNRS, MPI-DuMAS, ivan.iordanoff@ensam.eu</footnote>

<footnote n="3">³ I2M UMR 5295 CNRS, DuMAS, jl.charles@i2m.u-bordeaux1.fr</footnote>
Cela explique la décision de développer une nouvelle démarche qui puisse simuler la mécanique d'un continuum en utilisant la méthode DEM étant à la fois rapide, efficace et capable de gérer l'élasto-plasticité et les grands déplacements.

## 2 Description de la méthode

Le type de méthode des éléments discrets utilisée pour la simulation des matériaux élasto-plastiques endommageables est celle connue comme « Smooth DEM » ou « DEM explicite », dont les caractéristiques principales sont :

- Le comportement du matériau est défini par la loi d'interaction entre les éléments discrets.

- La force résultante exercée sur un élément discret est égale à la somme des forces exercées par ses voisins.

- Les interactions entre éléments discrets sont indépendantes.

- Un algorithme d'intégration explicite est utilisé pour faire des calculs dynamiques.

Le fait que les interactions entre les éléments discrets soient indépendantes permet de les calculer par paires pour après faire l'addition à chaque itération. Par exemple, la force résultante $F_{ri}$ d'un élément discret $i$ avec $N_v$ voisins est :

$$F_{ri} = \sum_{j=1}^{N_v} F_{ij}$$

$F_{ij}$ étant la force exercée par chacun des voisins $j$ sur l'élément discret regardé $i$.

### 2.1 Loi d'interaction

La loi d'interaction entre éléments discrets utilisée dans un premier temps est une force issue d'un potentiel de Lennard-Jones :

$$F_{ij}(d) = A \left( \frac{d_{eq}^{12}}{(Bd-C)^{13}} - \frac{d_{eq}^{6}}{(Bd-C)^{7}} \right) \quad (1)$$

Où $d$ est la distance entre les éléments discrets, $d_{eq}$ est la distance d'équilibre et $A$, $B$ et $C$ sont des paramètres qui vont définir la forme de la courbe. Ce type de forces est fréquemment utilisé dans les problèmes de dynamique moléculaire et elle a été aussi utilisée par [5] entre autres pour la simulation DEM mécanique des milieux continus.

Une des caractéristiques de la DEM explicite c'est la diminution du pas de temps critique selon la relation :

$$\Delta t_{\text{critique}} \propto \sqrt{\frac{M_{\text{min}}}{K_{\text{max}}}}$$

Où $K_{max}$ est la raideur maximale et $M_{min}$ est la masse équivalente minimale entre 2 éléments discrets en contact. Une des faiblesses de l'équation 1 est la forte augmentation de la raideur quand les éléments discrets en contact sont proches.

Afin d'éviter cette faiblesse, une simplification de la loi d'interaction a été testée.

<figure><img src="image_3.png" /><figcaption>Fig. 2 – Loi d'interaction entre deux éléments discrets.</figcaption></figure>
<table><caption>Tableau 1 – Paramètres liés à la loi d'interaction.</caption><thead><tr><td>PARAMÈTRE</td><td>DESCRIPTION</td><td>UNITÉ</td></tr></thead><tbody><tr><td>E<sub>μ</sub></td><td>Module d'Young du contact entre deux éléments discrets avant d'arriver au plateau de force maximale d'attraction. E<sub>μ</sub> = K/R<sub>moy</sub>.</td><td>N/m<sup>2</sup></td></tr><tr><td>γ</td><td>Force maximale d'attraction entre deux éléments discrets.</td><td>N</td></tr><tr><td>d<sub>max</sub></td><td>Distance maximale d'interaction entre les deux éléments discrets. Si la distance entre les éléments discrets est plus grande ils n'interagissent pas.</td><td>m</td></tr></tbody></table>

La loi d'interaction qui va définir notre modèle mécanique est schématisée dans la figure 2. Étant donné que la somme des rayons des deux éléments discrets $r_i + r_j$ ne dépend que de la géométrie des deux éléments discrets en contact, la loi proposée dépend des paramètres de la table 1.

À part ces 3 paramètres liés à la courbe de la force d'interaction, un paramètre sans dimension $C_s$ est utilisé pour sélectionner les éléments discrets appartenant aux surfaces (paragraphe 2.2). Un paramètre $c_d$ est également utilisé pour calculer la force de dissipation au contact et introduire de l'amortissement numérique.

## 2.2 Détermination des surfaces

Afin de déterminer si un élément discret $i$ appartient à une surface on va tracer une sphère imaginaire autour de lui et regarder l'ensemble des éléments discrets dont leur centre de gravité est placé dedans (figure 3). Un indicatif de la convenance de définir l'élément discret comme étant « de surface » va être la distance entre le barycentre de cet ensemble et la position de l'élément discret regardé (figure 4). En effet, plus cette distance est grand, plus les éléments discrets qui entourent l'élément discret $i$ vont être d'un côté d'un plan imaginaire qui peut définir une surface. L'élément discret concerné est alors un élément discret de surface. Cette méthode permet de détecter la création de fissures.

L'algorithme est divisé donc en deux parties :

- Dans une première partie l'ensemble des éléments discrets des alentours de $i$ va être identifié et en suite le barycentre de cet ensemble va être calculé.

- Finalement la distance entre le barycentre et le centre de gravité de $i$ va définir la appartenance à une surface selon un critère qui sera expliqué par la suite.

<figure><img src="image_4.png" /><figcaption>Fig. 3 – Détail du domaine discret avec l'élément discret $i$ en bleu, le rayon d'interaction mécanique en noir et le rayon de recherche en rouge.</figcaption></figure>

Il faut remarquer que chaque élément discret du domaine a un rayon d'action mécanique $R_a$ qui est plus grand que son rayon géométrique ($R_a \approx 1,2R_g$ dans [2]). Ce fait va limiter le choix du rayon de recherche des éléments discrets « indicatifs » car il n'est pas désirable qu'un élément discret qui
appartient à une surface opposée joue aucun rôle sur le critère utilisé. Cela limite alors la recherche de l'ensemble à ces éléments discrets qui se trouvent dans un rayon 2 fois plus grand que le rayon d'action $R_a$ de l'élément discret $i$ (figure 3).

Une fois que les éléments discrets de référence ont été choisies et leur barycentre a été calculé, il faut introduire dans le modèle un paramètre $C_s$ qui va permettre de définir l'appartenance de l'élément discret à une surface. Ce paramètre va définir une sphère autour de l'élément discret et si le barycentre est à l'extérieur de la sphère il sera considéré comme appartenant à une surface. La figure 4 montre un exemple en 2D pour un rayon d'action égal au rayon de l'élément discret.

<figure><img src="image_5.png" /><figcaption>Fig. 4 – Barycentre des éléments discrets en interaction avec l'élément discret $i$. Le critère marqué en jaune considère $i$ ($C_s = 0,3$) comme étant un élément discret de surface au contraire du critère rouge ($C_s = 0,8$).</figcaption></figure>

## 2.3 Modification de la loi de comportement

Une fois que les éléments discrets qui appartiennent aux surfaces sont détectés, la loi d'interaction mécanique peut être modifiée selon ce critère. Dans la figure 5 on peut voir un élément discret $i$ (en bleu) et la loi de comportement qu'il aurait s'il serait en contact avec un élément discret qui est placé du même coté d'une fissure (en orange) et un élément discret qui est placé dans le lèvre opposé de la fissure (en rouge). Dans ce dernier cas, le contact est unilatéral, ce qui caractérise la formation d'une fissure dans le matériau.

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

<figure><img src="image_7.png" /><figcaption>Fig. 5 – Éléments discrets dans une fissure et modification de la loi de comportement.</figcaption></figure>

# 3 Adimensionnement

Afin que deux matériaux aient le même comportement mécanique indépendemment de la taille moyenne des ses éléments discrets, les valeurs adimensionnés des paramètres qui définissent leur comportement mécanique doivent être identiques.

**Longueur** Le rayon moyen des éléments discrets a été choisi comme longueur de référence.
$$L_{\mu} = R_{moy}$$

**Masse** La masse de référence a été calculée à partir de la densité du matériau continu.

$$M_{\mu} = \rho_{0} \cdot R_{moy}^{3}$$

**Temps** Le pas de temps de référence est calculé à partir du module d'Young $E_{\mu}$. Si on fait l'analyse dimensionnel :

$$ 
\begin{align*}
[E_{\mu}] &= N \, m^{-2} = \text{kg} \, \text{m} \, s^{-2} \, \text{m}^{-2} \\
&= \text{kg} \, \text{m}^{-3} \cdot \text{m}^{2} \cdot \text{s}^{-2} \\
&= [\rho_{0}] \cdot [R_{moy}]^{2} \cdot [T_{\mu}]^{-2}
\end{align*} $$

Le pas de temps de référence retenu est alors :

$$T_{\mu} = R_{moy} \sqrt{\frac{\rho_{0}}{E_{\mu}}}$$

## 4 Étude de sensibilité

Des éprouvettes parallélépipèdiques de 30000 éléments discrets avec un rayon moyen de 53,180 µm ont été créées afin de mener les essais numériques de traction uniaxiale. Les dimensions nominales des éprouvettes sont 5 mm × 3 mm × 2 mm.

Les propriétés du matériau de base sont les suivantes :

<table><thead><tr><th>K (kN/m)</th><th>γ (N)</th><th>d<sub>max</sub> (R<sub>moy</sub>)</th><th>C<sub>s</sub> (-)</th><th>c<sub>d</sub> (N s<sup>2</sup>/m<sup>2</sup>)</th><th>ρ<sub>0</sub> (kg/m<sup>3</sup>)</th></tr></thead><tbody><tr><td>1500</td><td>25</td><td>2,4</td><td>0,4</td><td>0,08</td><td>2790</td></tr></tbody></table>

Différents essais numériques de traction vont être menés afin d'étudier l'influence de chacun des paramètres dans le comportement mécanique du matériau. Les différents valeurs adoptées sont montrées dans la table suivante :

<table><thead><tr><th>K (kN/m)</th><th>γ (N)</th><th>d<sub>max</sub> (R<sub>moy</sub>)</th><th>C<sub>s</sub> (-)</th><th>c<sub>d</sub> (N s<sup>2</sup>/m<sup>2</sup>)</th></tr></thead><tbody><tr><td>1000~8000</td><td>15~70</td><td>2,2~2,6</td><td>0,2~∞</td><td>0,04~0,06</td></tr></tbody></table>

Après avoir testé 4 éprouvettes différentes en matériau de base, les résultats de l'étude de sensibilité aux variations des paramètres sont montrés dans les figures suivantes :

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

<figure><img src="image_9.png" /></figure>
<figure><img src="image_10.png" /></figure>

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

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

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

- La variation du paramètre K fait notamment varier le module d'Young du matériau ainsi que la déformation à rupture.

- La variation du paramètre $\gamma$ augmente la contrainte maximale, mais fait aussi varier le module d'Young et la déformation à rupture.

- La principale influence de la variation de $d_{max}$ est le changement de la déformation maximale du matériau, mais la contrainte maximale et le module d'Young sont aussi affectés.

- Le paramètre $C_s$ joue un rôle principal dans la fragilité du matériau. Si ce paramètre est petit le matériau a un comportement fragile tandis que sa valeur es élevée il sera ductile et sa rupture progressive.

- Le paramètre $c_d$ a peu d'importance dans le comportement macroscopique du matériau.

# 5 Conclusions et perspectives

La campagne d'essais menée montre un bon comportement qualitatif du comportement élasto-plastique endommageable pour un essai de traction. Néanmoins, les paramètres de la loi de comportement mécanique utilisée dans ce travail sont fortement couplés. Cela rend difficile une calibration directe des paramètres microscopiques pour obtenir les propriétés macroscopiques désirées.

# Références

[1] P. Cundall. *A computer model for simulating progressive large scale movements of blocky rock systems*, In Proceedings of the symposium of the International Society of Rock Mechanics, 1971.

[2] F.A. Tavarez, M.E. Plesha. *Discrete element method for modelling solid and particulate materials* International Journal for Numerical Methods in Engineering, 70, 379–404, 2006.

[3] D. Richard, I. Iordanoff, M. Renouf, Y. Berthier. *Thermal study of the dry sliding contact with third body presence*, ASME Journal of Tribology, 130, Issue 3(031404):10, 2008.

[4] F. Chinesta, S. Cescotto, E. Cueto, P Lorong. *Natural element method for the simulation of structures and processes*, ISTE-Wiley, London, 2010.

[5] Z.P. Tang. *Three-dimensional DEM theory and its application to impact mechanics*, Science in China, 44(6), 2001.
[6] F. Fleissner, T. Gaugele, P. Eberhard. *Application of the discrete element method in mechanical engineering*, Multibody System Dynamics, 18:81–94, 2007.

[7] D. André, I. Iordanoff, J. L. Charles, J. Néauport. *Discrete element method to simulate continuous material by using the cohesive beam model*, Computer Methods in Applied Mechanics and Engineering, 213–216:113–125, 2012.

[8] F. Camborde, C. Mariotti, F. V. Donzé. *Numerical study of rock and concrete behaviour by discrete element modelling*, Computers and Geotechnics, 27:225-247, 2000.