# Modélisation par la méthode des éléments discrets du comportement élastique de milieux continus hétérogènes

Hamza Haddad, Willy Leclerc, Mohamed Guessasma, Nabil Ferguen, Christine Pelegris, Emmanuel Bellenger

## ► To cite this version:

Hamza Haddad, Willy Leclerc, Mohamed Guessasma, Nabil Ferguen, Christine Pelegris, et al. Modélisation par la méthode des éléments discrets du comportement élastique de milieux continus hétérogènes. 12e Colloque national en calcul des structures, CSMA, May 2015, Giens, France. hal-01517277

HAL Id: hal-01517277

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

Submitted on 3 May 2017

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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
# Modélisation par la méthode des éléments discrets du comportement élastique de milieux continus hétérogènes

H. Haddad<fnref n="1" />, W. Leclerc<fnref n="2" />, M. Guessasma<fnref n="2" />, N. Ferguen<fnref n="2" />, C. Pélegris<fnref n="2" />, E. Bellenger<fnref n="2" />

**Résumé** — La méthode des éléments discrets s'avère une approche adaptée pour la simulation du comportement mécanique de milieux continus homogènes pour lesquels des phénomènes locaux tels que l'endommagement et la fissuration interviennent. Dans cette communication, nous nous intéressons aux milieux continus hétérogènes pour lequel un travail important de validation reste encore à réaliser. Nous nous limitons au domaine de l'élasticité linéaire en 2D et considérons le cadre des éléments discrets cohésifs qui couple le modèle particulaire à un modèle lattice de type poutre. Des tests de validation sont réalisés et des comparaisons effectuées avec des méthodes classiques d'homogénéisation pour différentes configurations et contrastes de propriétés.

**Mots clés** — Méthode des éléments discrets, matériaux composites, propriétés élastiques, homogénéisation.

## 1 Introduction

Le comportement thermomécanique des matériaux homogènes ou hétérogènes est usuellement simulé à l'aide de méthodes basées sur la mécanique des milieux continus. Les phénomènes intervenant au sein du matériau comme la fissuration ou l'endommagement nécessitent une prise en compte naturelle des discontinuités induises par des interfaces de forme géométrique souvent complexe. Les modèles discrets se révèlent donc mieux adaptés pour simuler ce type de phénomène. Afin de modéliser un domaine continu à l'aide d'éléments discrets, une méthodologie consiste à associer au modèle particulaire un modèle lattice [5] permettant de simuler la cohésion entre les éléments discrets, généralement de forme sphérique en 3D, ou cylindrique en 2D. La liaison peut être réalisée à l'aide d'éléments de type ressort mais il a été montré que le modèle poutre permet une meilleure simulation des géométries de fissuration [7]. Ainsi, André et al. [1] ont récemment utilisé cette approche afin de modéliser le comportement mécanique au sein d'un milieu continu homogène, et Haddad et al. [4] l'ont couplé avec la méthode des éléments finis pour simuler le transfert thermique par conduction. Notre objectif est d'adapter et valider le modèle hybride particulière-lattice pour les matériaux continus hétérogènes. Dans le cadre de cet article, nous nous limitons au domaine de l'élasticité linéaire avec des liaisons parfaites aux interfaces. Deux configurations 2D sont étudiées, le cas d'une inclusion circulaire noyée au sein d'une matrice, et celui d'un composite particulaire. Des tests de validation sont réalisés et des comparaisons effectuées en termes de propriétés et de contraintes avec d'autres approches numériques et analytiques.

## 2 Cas d'un milieu continu homogène

### 2.1 Modèle de poutre 2D

Nous considérons le modèle hybride particulière-lattice pour lequel les particules sont représentées par des cylindres droits parallèles et les intéractions entre particules par des poutres de longueur $L_μ$, de module de Young $E_μ$, de section droite $A_μ$ et de moment d'inertie $I_μ$. Nous supposons un contact entre particules cohésif sans critère de rupture. Le modèle de poutre considéré ainsi que les intéractions associées entre deux particules i et j sont décrits ci-dessous :

<footnote n="1">¹ LAMFA, Université de Picardie Jules Verne, hamza.haddad@u-picardie.fr</footnote>

<footnote n="2">² LTI, Université de Picardie Jules Verne, {willy.leclerc,mohamed.guessasma,nabil.ferguen,christine.pelegris, emmanuel.bellenger}@u-picardie.fr</footnote>
<figure><img src="image_2.png" /><figcaption>FIGURE 1 – Modèle de poutre 2D</figcaption></figure>

où $F_n^{i \to j}$, $F_t^{i \to j}$ et $M_{i \to j}^{int}$ sont respectivement les efforts d'intéraction normal et tangentiel, et le moment d'intéraction de la particule $i$ sur $j$. $\theta_{i,j}$, $u_n^{i,j}$ et $u_t^{i,j}$ représentent respectivement la rotation et les déplacements normal et tangentiel. La section droite $A_\mu$ est rectangulaire de côté $e$ et $h$, où $h$ dépend des rayons des particules en intéraction $R_i$ et $R_j$ et d'un paramètre strictement positif $r_\mu = \frac{2h}{R_i+R_j}$.

$$ \begin{bmatrix} F_n^{i \to j} \\ F_t^{i \to j} \\ M_{i \to j}^{int} \end{bmatrix} = \begin{bmatrix} \frac{E_\mu A_\mu}{L_\mu} & 0 & 0 \\ 0 & \frac{12E_\mu I_\mu}{L_\mu^2} & \frac{6E_\mu I_\mu}{L_\mu^2} \\ 0 & \frac{6E_\mu I_\mu}{L_\mu^2} & \frac{4E_\mu I_\mu}{L_\mu} \end{bmatrix} \begin{bmatrix} u_n^i - u_n^j \\ u_t^i - u_t^j \\ \theta_i \\ \theta_j \end{bmatrix} $$

## 2.2 Calibration et validation

Les paramètres macroscopiques du milieu continu homogène sont déduits des paramètres microscopiques $E_\mu$ et $r_\mu$ à l'aide d'un test de traction quasi-statique lequel mène à la détermination du module de Young macroscopique $E_M$ et du coefficient de Poisson $v_M$. Des abaques ont été construits afin de déterminer le comportement local à partir des paramètres macroscopiques. Les tests effectués ont été réalisés à l'aide d'empilements granulaires aléatoires générés par la méthode de compression isotropique des parois au sein d'un domaine rectangulaire d'élancement 4,78. Nous avons opté pour des systèmes discrets polydisperses lesquels sont moins sujets aux effets d'anisotropie. La fraction surfacique occupée par les cylindres est proche de 87% avec un nombre de coordination proche de 4,5 ce qui correspond à un random close packing en 2D [3]. Afin d'assurer l'isotropie du système, le nombre de particules a été fixé à 15000 après étude de la distribution des angles de contact. Les Figures 2 et 3 décrivent une évolution

<figure><img src="image_3.png" /><figcaption>FIGURE 2 – Evolution de $E_M$ en fonction de $r_\mu$</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>FIGURE 3 – Evolution de $v_M$ en fonction de $r_\mu$</figcaption></figure>

quadratique des paramètres macroscopiques $E_M$ et $v_M$ en fonction de $r_\mu$. Les résultats sont similaires à ceux de André et al. [1] obtenus dans le cadre de particules sphériques en 3D. Afin de valider la méthode de calibration, des essais quasi-statiques de traction ont été menés à la fois sur les modèles discret et continu de type éléments finis. Les Figures 4 et 5 illustrent l'évolution au cours du temps de la réaction sur le bord encastré et du déplacement sur le bord sollicité pour un milieu homogène de paramètres $E_M = 65\,\text{GPa}$ et $v_M = 0,247$.

<figure><img src="image_5.png" /><figcaption>FIGURE 4 – Réaction en fonction du temps</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIGURE 5 – Déplacement en fonction du temps</figcaption></figure>
# 3 Modèle hétérogène à une inclusion

## 3.1 Hypothèses et conservation de masse

Nous nous intéressons au cas d'un milieu continu hétérogène constitué d'une inclusion unique noyée dans une matrice. Nous considérons un motif représentatif carré pour lequel l'inclusion circulaire est centrée et a un rayon égal à 0,25 fois la longueur du motif. $E_μ$ à l'interface entre une particule i de masse $m_i$ au sein de la matrice et une particule j de masse $m_j$ au sein de l'inclusion est déterminé comme suit :

$$E_{\mu} = \frac{m_i E_{\mu}^{i} + m_j E_{\mu}^{j}}{m_i + m_j} \quad (1)$$

Le nombre de particules à générer doit être au minimum de 3000 afin de respecter la densité de particules considérée dans le cas homogène. Nous avons vérifié que ce nombre suffit pour respecter la masse et le volume de l'inclusion après discrétisation. Cependant, à des fins de précision, le calcul de contraintes nécessite un nombre total plus important estimé à 7000 particules.

## 3.2 Etudes et comparaisons

Les contraintes sont estimées à l'aide de la formule de Love-Weber sans terme inertiel via un pavage du domaine discret en 332 volumes élémentaires triangulaires [2,6]. Nous avons estimé qu'un minimum de 20 particules doit être considéré dans chaque volume élémentaire, soit donc le total avancé de 7000 particules. Des tests quasi-statiques de traction et de cisaillement simple avec conditions aux limites de symétrie et d'antisymétrie respectivement sont réalisés sur le motif représentatif. Le module de Young $E^m$ de la matrice est fixé à 65 GPa et différents modules de Young $E^i$ sont considérés pour l'inclusion. Ainsi, le contraste de propriétés $c_r = E^i/E^m$ est varié de 0,02 à 50. Les Figures 6 et 7 décrivent le

<figure><img src="image_7.png" /><figcaption>FIGURE 6 – Contraintes de Von-Mises (MED - cas de l'inclusion seule)</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>FIGURE 7 – Contraintes de Von-Mises (MEF - cas de l'inclusion seule)</figcaption></figure>

champ de contraintes de Von-Mises obtenu dans le cas du test de traction avec la MED et la MEF. Nous considérons ici les propriétés d'un composite $BaTiO_3$/époxy avec $E^m = 65GPa$ et $E^i = 2,9GPa$. Les résultats sont qualitativement très proches avec un maximum cependant supérieur (15,5 MPa contre 11,1 MPa) dans le cas de la MEF du fait d'une discrétisation plus fine (720000 éléments triangulaires à trois noeuds). Les modules de Young $E$ et de cisaillement $G$ effectifs du motif élémentaire sont estimés à partir des tests quasi-statiques puis adimensionnés par rapport aux propriétés de la matrice ($E^a = E/E^m$, $G^a = G/G^m$). Des comparaisons sont réalisées avec la MEF dans la même configuration : la MEF avec conditions aux limites périodiques (2EMEF), la méthode de Transformée de Fourier Rapide (TFR) et le modèle de Mori-Tanaka (MT). Les Figures 8 et 9 décrivent l'influence de $c_r$ sur $E^a$. Les résultats montrent une bonne adéquation entre les différents modèles quel que soit le contraste considéré. Des résultats similaires ont été obtenus pour $G^a$.

# 4 Modèle de composite particulaire

Nous considérons un modèle de composite particulaire composé d'inclusions monodisperses en contact ponctuel avec une fraction surfacique d'inclusions fixée à 85% ± 0,5. Un nombre de 40 inclusions est retenu afin d'assurer la représentativité du milieu. Le nombre d'éléments discrets au sein du système granulaire est fixé à 65000 particules selon l'étude réalisée sur l'inclusion seule. Les contraintes de Von-Mises sont estimées selon la même méthodologie que celle décrite pour l'inclusion seule. Nous
<figure><img src="image_9.png" /><figcaption>FIGURE 8 – $E^a$ vs $c_r$ ($E^a > 1$)</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>FIGURE 9 – $E^a$ vs $c_r$ ($E^a < 1$)</figcaption></figure>

considérons un pavage triangulaire du domaine composé de 3300 éléments de sorte que chaque élément contienne en moyenne 20 particules. Des tests quasi-statiques de traction et de cisaillement sont réalisés de façon identique au cas de l'inclusion seule mais nous nous limitons ici aux propriétés du composite $BaTiO_3$/époxy. Les Figures 10 et 11 illustrent les contraintes de Von-Mises obtenues au sein du milieu hétérogène pour le test de traction avec la MED et la MEF. La contrainte maximale est de 6,37 MPa dans le cas de la MEF contre 8,04 MPa dans le cas de la MED. Cette différence s'explique par la présence aux interfaces des inclusions d'éléments triangulaires contenant parfois moins de 5 particules. Les modules

<figure><img src="image_11.png" /><figcaption>FIGURE 10 – Contraintes de Von-Mises (MED - cas du composite)</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>FIGURE 11 – Contraintes de Von-Mises (MEF - cas du composite)</figcaption></figure>

de Young $E$ et de cisaillement $G$ effectifs sont estimés à l'aide de la MED et les prédictions sont comparées avec la MEF, la 2EMEF et la TFR. La Table 1 montre que les coefficients élastiques obtenus avec la MED sont en accord avec ceux obtenus avec ces différentes approches.

<table><caption>TABLE 1 – Prédicitions numériques de *E* et *G* pour le modèle de composite</caption><thead><tr><th></th><th>MEF</th><th>2EMEF</th><th>TFR</th><th>MED</th></tr></thead><tbody><tr><td><i>E</i> (GPa)</td><td>4,30</td><td>4,47</td><td>4,20</td><td>4,31</td></tr><tr><td><i>G</i> (GPa)</td><td>1,64</td><td>1,56</td><td>1,63</td><td>1,66</td></tr></tbody></table>

Références

[1] D. André, I. Iordanoff, JC. Charles, J. Neauport *Discrete element method to simulate continuous material by using the cohesive beam model*, Computational Methods in Applied Mechanics and Engineering, 213, 113-128, 2012.

[2] G. De Saxcé, J. Fortin, O. Millet *About the numerical simulation of the dynamics of granular media and the definition of the mean stress tensor*, Mechanics of Materials, 36(12), 1175-1184, 2004.

[3] A. Donev, I. Cissé, D. Sachs, EA. Variano, FH. Stillinger, R. Connelly et al. *Improving the density of jammed disordered packings using ellipsoids*, Science, 303, 990-993, 2004.

[4] H. Haddad, M. Guessasma, J. Fortin *Heat transfer by conduction using DEM-FEM coupling method*, Computational Materials Science, 81, 339-347, 2013.

[5] A. Ibrahimbegovic, A. Delaplace *Microscale and mesoscale discrete models for dynamic fractures of structures built of brittle material*, Computers and Structures, 81(21), 1225-1265, 2003.

[6] F. Nicot, N. Hadda, M. Guessasma, J. Fortin, O. Millet *On the definition of the stress tensor in granular media*, International Journal of Solids and Structures, 50(14-15), 2508-2517, 2013.

[7] E. Schlangen, E.J. Garboczi *New method for simulating fracture using an elasticity uniform random geometry lattice*, International Journal of Engineering Science, 34(10), 1131-1144, 1996.