# Simulation du comportement thermique d'un milieu continu équivalent par méthode des éléments discrets

Jérôme Rivière, Mathieu Renouf, Yves Berthier

► To cite this version:

Jérôme Rivière, Mathieu Renouf, Yves Berthier. Simulation du comportement thermique d'un milieu continu équivalent par méthode des éléments discrets. 11e colloque national en calcul des structures, CSMA, May 2013, Giens, France. hal-01717076

HAL Id: hal-01717076

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

Submitted on 25 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 du comportement thermique d'un milieu continu équivalent par méthode des éléments discrets

Jérôme RIVIERE<fnref n="1" />, Mathieu RENOUF<fnref n="2" />, Yves BERTHIER<fnref n="1" />

## Résumé —

La méthode des éléments discrets est largement utilisée pour représenter non seulement le comportement mécanique de milieux hétérogènes et discontinus tels que les milieux divisés (grains, sols, maçonneries) mais aussi leur comportement thermique. Toutefois, lorsqu'elle est utilisée pour représenter un milieu continu équivalent, les approches thermo-mécaniques au sens discret ne sont plus utilisables en tant que telles. Ce travail a pour objectif d'étendre la méthode des éléments discrets afin de modéliser le transfert de chaleur par conduction dans un milieu continu équivalent. L'influence de l'arrangement spatial des éléments est étudiée afin que l'évolution de la température du système ne dépende pas de cette discrétisation. Les résultats sont comparés à une approche continue de type éléments finis. Pour conclure, le cas d'un milieu endommagé est traité afin de pouvoir effectuer une transition entre milieux continu et discret.

**Mots clés** — méthode des éléments discrets, milieu continu équivalent, thermique.

# 1 Introduction

D'un point de vue physique, en fonction de l'échelle d'observation, un milieu microscopiquement divisé peut être vu macroscopiquement comme un milieu continu. A l'inverse chaque grain, issu de la division, peut lui-même être vu comme un milieu continu. En tribologie, ces deux points de vus entrent dans la description du triplet tribologique où les premiers corps (corps en contact) et troisième corps (interface) coexistent [1] et où des phénomènes de cohésion [2] peuvent rendre continus des corps initialement discontinus.

D'un point de vue numérique, l'une des possibilités offertes par la méthode des éléments discrets peut être de modéliser des phénomènes ayant lieu à des échelles discrètes ou continues à l'aide d'une collection d'objets. La difficulté réside en la détermination des lois d'interaction appropriées entre éléments permettant de retrouver le comportement global du milieu équivalent. On trouve de nombreux travaux traitant de ce type de modélisation utilisant par exemple des lois de contact interparticulaires dérivées de techniques d'homogénéisation [3], des interactions de type poutre [4] ou encore des lois bilatérales endommageables [5, 6]. C'est sur ce type d'approche que nous appuierons notre travail nous permettant ainsi, grâce à une variable $β$, de caractériser l'endommagement de chaque liaison ($β = 0$, le lien est cassé et $β = 1$, le lien est intègre).

En se focalisant sur les phénomènes de transfert de chaleur dans les solides, il existe également différents travaux traitant du comportement thermique équivalent d'un milieu par approche discrète. Cependant, certaines de ces méthodes utilisent un coefficient "artificiel" sans réelle signification physique pour faire coïncider les résultats issus de la méthode des éléments discrets avec des méthodes continues [7]. Plus récemment, d'autres travaux [8] tentent d'introduire un effet de microstructure au travers de nombreux paramètres qu'il est difficile d'appréhender. C'est dans cette lignée que s'articule notre travail, en proposant une conductance au contact entre éléments discrets reliée directement à la microstructure du milieu. Nous cherchons donc à déterminer son comportement thermique équivalent et le traiterons sous trois types de représentations : discrète ($β = 0$), continue ($β = 1$) et endommagée ($β = 0, 1$).

<footnote n="1">¹ Université de Lyon, CNRS, LaMCoS, INSA Lyon, UMR 5259, F-69621 Villeurbanne Cedex, France {Jerome.Riviere,Yves.Berthier}@insa-lyon.fr</footnote>

<footnote n="2">² Université Montpellier 2, CNRS, LMGC, UMR 5508, F-34095 Montpellier Cedex, France, Mathieu.Renouf@univ-montp2.fr</footnote>
## 2 Étude thermique incluse dans une représentation discrète d'un milieu

La conduction thermique dans un milieu granulaire nécessite de pouvoir déterminer quelle quantité de chaleur va être échangée au contact. Suite aux travaux de Vargas et al. [9, 10], le flux de chaleur $Q_{ij}$ traversant deux sphères rigides indépendantes i et j est conditionné par la conductance au contact H. Cette dernière correspond à la constriction des lignes de flux ainsi qu'à leur température respectives $T_i$ et $T_j$ :

$$Q_{ij} = H(T_i - T_j). \qquad (1)$$

En approximant les solutions analytiques proposées par Yovanovitch [11], il est possible de relier la conductance H aux forces de contact en utilisant la théorie élastique de Hertz :

$$H = 2k_{th}\sqrt{2aL} = 2k_{th}\left(\frac{8F_n r^*}{\pi L E^*}\right)^{\frac{1}{4}}, \qquad (2)$$

où $k_{th}$ représente la conductivité thermique des corps rigides en contact, L, la longueur du contact (unitaire dans un modèle granulaire bi-dimensionnel) et a, le rayon de contact. La théorie élastique de Hertz permet d'écrire a en fonction de la force normale $F_n$ entre les deux particules i et j, du rayon effectif $r^*$ et du module de Young effectif $E^*$ du couple de particules.

En considérant les multiples interactions qu'une particule i peut subir avec n particules voisines, le flux de chaleur $Q_i$ reçu vaut donc :

$$Q_i = \sum_{j=1}^{n} Q_{ij}. \qquad (3)$$

L'évolution temporelle de la température $T_i$ s'exprime alors par :

$$\frac{dT_i}{dt} = \frac{Q_i(t)}{\alpha_i}, \qquad (4)$$

où $\alpha_i = \rho_i c_i V_i$ représente la capacité thermique massique de la particule i. Cette approche considère que les gradients thermiques sont confinés près du point de contact entre les deux sphères de rayon r (nombre de Biot : $Bi = 2a/\pi r \ll 1$) et que le flux thermique est quasi-statique ($dT_i/(T_j - T_i) \ll 1$).

## 3 Étude thermique incluse dans une représentation continue d'un milieu

### 3.1 Formulation

De façon naïve, il est facile de faire une analogie entre un empilement de disques sur un réseau carré dans une boîte rectangulaire et une discrétisation par différences finies du même milieu (cf. figure 1). On peut ainsi décrire la température d'un noeud (c'est à dire d'une particule) de deux façons différentes : soit en utilisant l'approche présentée dans la section 2, soit par une approche continue type différences finies.

<figure><img src="image_2.png" /><figcaption>Fig. 1 – Milieu continu discrétisé par (a) des disques empilés, (b) un maillage régulier et (c) leur analogie</figcaption></figure>
L'équation (2) n'a de sens que si le milieu est naturellement divisé. Si la division devient numérique, au sens méthode sans maillage, elle ne doit pas influencer les comportements mécanique et thermique. Nous proposons ici d'utiliser la microstructure de l'empilement afin d'avoir le comportement thermique équivalent d'un milieu continu. Au travers de l'équation (1), nous déduisons que c'est dans la définition de H que va résider cette équivalence.

Dans ce cadre, on considère une plaque solide de longueur L et d'épaisseur $2R$ où R est le rayon d'une particule servant à discrétiser le domaine. D'après la section 2, la température $T_i$ de la particule $P_i$ est égale à :

$$T_i^+ = T_i^- + \frac{dt}{\alpha_i} [Q_{i,i-1}^- + Q_{i,i+1}^-], \qquad (5)$$

$$\Leftrightarrow T_i^+ = T_i^- + \frac{dt}{\alpha_i} \tilde{H} [(T_{i-1}^- - T_i^-) + (T_{i+1}^- - T_i^-)], \qquad (6)$$

en supposant que $\tilde{H}$, la conductance au contact recherchée est le même pour $Q_{i,i-1}$ et $Q_{i,i+1}$. On peut alors réécrire l'équation (6) comme :

$$T_i^+ = \left[ 1 - \frac{2dt}{\alpha_i} \tilde{H} \right] T_i^- + \frac{dt}{\alpha_i} \tilde{H} [T_{i-1}^- + T_{i+1}^-]. \qquad (7)$$

L'évolution de la température en un point de la même plaque peut être également donnée par l'équation de la chaleur unidimensionnelle :

$$k_{th} \frac{d^2 T}{dx^2} = \rho c \frac{dT}{dt}, \qquad (8)$$

où $k_{th}$ représente la conductivité thermique de la matière, $\rho$, sa masse volumique et $c$, sa capacité calorifique. Une approximation de l'équation de la chaleur (8) nous donne :

$$k_{th} \frac{(T_{i+1}^- + T_{i-1}^- - 2T_i^-)}{dx^2} = \rho c \frac{T_i^+ - T_i^-}{dt}. \qquad (9)$$

Ici, $dx$ qui représente la discrétisation spatiale est donc égale à $2R$. On peut alors écrire que :

$$T_i^+ = \frac{dtk_{th}}{\rho c} \left[ \frac{(T_{i+1}^- + T_{i-1}^- - 2T_i^-)}{(2R)^2} \right] + T_i^-, \qquad (10)$$

$$\Leftrightarrow T_i^+ = \left[ 1 - \frac{dtk_{th}\pi}{2\alpha_i} \right] T_i^- + \frac{dtk_{th}\pi}{4\alpha_i} [T_{i+1}^- + T_{i-1}^-]. \qquad (11)$$

Par analogie entre les deux modèles représentés par les équations (7) et (11), on en déduit que :

$$\tilde{H} = \frac{k_{th}\pi}{4}. \qquad (12)$$

Il est intéressant de noter que la conductance au contact $\tilde{H}$ ne fait intervenir apparemment aucun paramètre concernant la microstructure du milieu.

## 3.2 Comparaison de modèles

Afin d'évaluer la pertinence du modèle donné par l'équation (12), on construit une plaque solide de 100 mm par 40 mm par empilement granulaire, en utilisant 3 discrétisations différentes : un empilement carré, un autre hexagonal et un dernier aléatoire où la polydispersité est de 20 % autour du rayon moyen (cf. figure 2). Nous considérons une collection de disques rigides contenus dans une boîte rectangulaire et nous faisons l'hypothèse de transferts thermiques évoluant strictement par conduction entre éléments et sans gradients thermiques dans leur volume. Le milieu étudié est statique et purement homogène.

Initialement, ce dernier est à une température de 300 K et des conditions adiabatiques sont fixées aux bords. On le soumet ensuite à une source de chaleur de 1000 K, constante dans le temps, sur la paroi
inférieure. L'évolution de la température est mesurée à différentes distances de la source de chaleur. Les paramètres des simulations et les propriétés volumiques des particules, indépendantes de la température, sont répertoriés dans les tableaux 1 et 2.

<figure><img src="image_3.png" /><figcaption>Fig. 2 – Empilements discrets après 60 s de chauffe : (a) carré, (b) hexagonal, (c) aléatoire</figcaption></figure>

<table><caption>Tableau 1 – Paramètres de simulation</caption><thead><tr><th>Empilement</th><th>Carré</th><th>Hexagonal</th><th>Aléatoire</th></tr></thead><tbody><tr><td>Nombre de particules</td><td>1000</td><td>1131</td><td>1000</td></tr><tr><td>Polydispersité (%)</td><td>0</td><td>0</td><td>20</td></tr><tr><td>Longueur de la plaque (mm)</td><td></td><td>100</td><td></td></tr><tr><td>Largeur de la plaque (mm)</td><td></td><td>40</td><td></td></tr><tr><td>Rayon des particules (mm)</td><td></td><td>1</td><td></td></tr><tr><td>Température initiale (K)</td><td></td><td>300</td><td></td></tr><tr><td>Température de chauffe (K)</td><td></td><td>1000</td><td></td></tr><tr><td>Temps de simulation (s)</td><td></td><td>60</td><td></td></tr><tr><td>Pas de temps (ms)</td><td></td><td>0,1</td><td></td></tr></tbody></table>

<table><caption>Tableau 2 – Propriétés volumiques des particules</caption><tr><td>Conductivité thermique k<sub>th</sub> (W.m<sup>-1</sup>.K<sup>-1</sup>)</td><td>110</td></tr><tr><td>Chaleur spécifique c (J.kg<sup>-1</sup>.K<sup>-1</sup>)</td><td>1500</td></tr><tr><td>Masse volumique ρ (kg/m<sup>3</sup>)</td><td>1750</td></tr><tr><td>Module de Young effectif E* (MPa)</td><td>30000</td></tr><tr><td>Coefficient de Poisson effectif ν*</td><td>0,25</td></tr></table>
Afin de valider notre modèle thermique, nous proposons de le confronter avec une simulation par éléments finis (FEM). Sous le logiciel Ansys™ 12.1, nous avons conçu un milieu ayant les mêmes propriétés géométriques et thermiques que celles utilisées dans la méthode des éléments discrets. Nous avons maillé le milieu de telle sorte que la longueur d'une maille fasse deux fois le rayon d'une particule numérique. Il y a donc 1000 éléments mesurant 2 mm chacun.

La figure 3 illustre l'évolution de la température au cours du temps mesurée à 25 mm et 75 mm de la source de chaleur. Nous notons un très bon accord entre le résultat de la simulation par éléments finis et celui de la simulation discrète avec un empilement carré et aléatoire. L'écart est en revanche beaucoup plus important avec un empilement hexagonal.

<figure><img src="image_4.png" /><figcaption>Fig. 3 – Évolution de la température à (a) 25 mm et (b) 75 mm de la source de chaleur pour trois types d'empilements avec $\tilde{H} = k_{th}\pi/4$ et comparaison avec la méthode des éléments finis (FEM)</figcaption></figure>

Cette différence est cependant naturelle car l'équation (12) est seulement valable pour un empilement (discrétisation) carré. De plus, comme indiqué dans le tableau 1, l'empilement hexagonal se caractérise par une plus grande compacité. De fait, un nombre plus important de particules est nécessaire pour discrétiser une même surface de la plaque. Ainsi, le nombre de contacts augmentant, les échanges de chaleur dans le milieu sont alors plus nombreux. Afin de généraliser $\tilde{H}$ pour tous types d'empilements, réécrivons $V_i$ et $dx$ de manière à faire apparaître des paramètres qui vont représenter la microstructure du milieu. Dans ce cadre, introduisons $N_b^{tot}$, le nombre total d'éléments discrets et $N_c^{tot}$, le nombre de contacts total entre les éléments discrets tels que :

$$ \left\{ \begin{array}{l} dx = r_i, \\ V_i = \frac{N_b^{tot}}{N_c^{tot}} \pi r_i^2. \end{array} \right. \qquad (13) $$

En observant que $N_c^{tot}/N_b^{tot} = <N_c>$ représente le nombre moyen de contacts dans le milieu, nous pouvons exprimer $\tilde{H}$ comme :

$$ \tilde{H} = \frac{k_{th}\pi}{<N_c>} \qquad (14) $$

De ce résultat, nous en concluons que la conductance au contact $\tilde{H}$ dépend exclusivement du nombre moyen de contacts possédé par chaque particule numérique du milieu. Ce dernier est relié à la compacité et donc à la manière dont elles sont arrangées dans le milieu. Vérifions cette expression précédente dans le cas de l'échauffement d'une plaque solide.

La figure 4 prend en compte l'expression (14) de $\tilde{H}$ dans le calcul de l'évolution de la température. Elle montre que cette température ne dépend maintenant que très peu de l'empilement granulaire et ceci, à différentes positions du milieu équivalent. Les légers écarts entre les différents empilements sont des effets de bords. En effet, dans le cas de l'empilement hexagonal, les disques se trouvant au milieu
du système ont chacun 6 voisins alors que les disques en contact avec les parois ont 4 voisins. Par conséquent, le nombre moyen de contacts dans le milieu hexagonal est de 5,81 et ceux des empilements carré et aléatoire sont respectivement de 4 et 4,07. Par rapport à une simulation par éléments finis, la figure 4 indique une très bonne corrélation avec les simulations par éléments discrets.

<figure><img src="image_5.png" /><figcaption>Fig. 4 – Évolution de la température à (a) 25 mm et (b) 75 mm de la source de chaleur pour trois types d'empilements avec $\tilde{H} = k_{th}\pi / <N_c>$ et comparaison avec la méthode des éléments finis (FEM)</figcaption></figure>

Enfin, dans le but de vérifier ce résultat sur le cas aléatoire, nous avons réalisé 5 empilements différents avec une polydispersité identique et égale à 20 % autour du rayon moyen (cf. figure 5).

<figure><img src="image_6.png" /><figcaption>Fig. 5 – Évolution de la température à 75 mm de la source de chaleur dans 5 empilements aléatoires polydisperses à 20 %</figcaption></figure>

Nous remarquons une superposition des évolutions de température. Ce constat souligne donc que l'équation (12) exprimant $\tilde{H}$ est bien appropriée dans l'objectif de s'abstraire de la distribution des éléments discrets dans la représentation du comportement thermique d'un milieu continu équivalent et bidimensionnel.

## 4 Etude thermique d'un milieu endommagé

Modéliser l'endommagement d'un milieu peut s'effectuer par exemple par des modèles de zone cohésive qui consistent à décrire des phénomènes locaux de fragmentation en conservant les approches
globales de la mécanique de la rupture [6]. Soit $\beta$, le paramètre global d'endommagement qui traduit le degré d'endommagement des liaisons entre les particules numériques. Dans le cas d'un tel milieu endommagé, $\beta \in ]0, 1[$ puisque des liaisons sont rompues entre la particule concernée et ses voisines. Nous reprenons alors l'équation (14) de la conductance au contact $\tilde{H}$ d'un milieu continu et plaçons ce paramètre $\beta$ en facteur :

$$ \tilde{H} = \frac{\beta k_{th} \pi}{< N_c >}. \qquad (15) $$

Nous illustrons le caractère plus isolant du milieu d'autant plus que des phénomènes d'endommagement tels que de la fissuration ou la présence de pores sont importants. Plus la liaison bilatérale entre deux particules est endommagée, plus la conduction de chaleur est difficile. Nous avons représenté l'influence du paramètre d'endommagement $\beta$ pour l'empilement carré sur la figure 6. A 30 s de simulation et 75 mm de la source de chaleur, dans un milieu endommagé à 75 % tel que $\beta = 0,25$, la température ne représente plus que 75 % de la température calculée pour un milieu continu équivalent. A 60 s, celle-ci ne vaut plus que 58 % (cf. tableau 3).

<figure><img src="image_7.png" /><figcaption>Fig. 6 – Évolution de la température de l'empilement carré pour un point situé à 75 mm de la source de chaleur en fonction du paramètre d'endommagement $\beta$</figcaption></figure>

<table><caption>Tableau 3 – Variations de température en fonction de l'endommagement du milieu à une distance de 75 mm de la source de chaleur</caption><thead><tr><th>β</th><th>Température à 30 s (K)</th><th>Part relative de la température calculée avec β = 1</th></tr></thead><tbody><tr><td>0,25</td><td>301</td><td>1 %</td></tr><tr><td>0,5</td><td>322</td><td>22 %</td></tr><tr><td>0,75</td><td>358</td><td>58 %</td></tr></tbody><thead><tr><th>β</th><th>Température à 60 s (K)</th><th>Part relative de la température calculée avec β = 1</th></tr></thead><tbody><tr><td>0,25</td><td>322</td><td>9 %</td></tr><tr><td>0,5</td><td>400</td><td>40 %</td></tr><tr><td>0,75</td><td>479</td><td>71 %</td></tr></tbody></table>

## 5 Conclusion

A travers cette étude, la méthode des éléments discrets a été utilisée pour décrire le comportement thermique d'un milieu en fonction de son endommagement. En partant de la méthode des différences
finies et d'expressions de la conductance au contact habituellement prises en compte pour modéliser des milieux granulaires rigides, nous avons déduit une solution correspondante à une représentation continue de la matière. Nous avons notamment montré son indépendance vis à vis de l'arrangement spatial des éléments discrets. Par ailleurs, en choisissant des modèles de zone cohésive, la solution a été modifiée pour évaluer l'évolution de la température dans un milieu endommagé. L'ensemble permettra de mieux décrire l'évolution du comportement thermique d'un contact sous sollicitations tribologiques (pression et cisaillement simultanées) où la matière et en particulier le troisième corps sont plus ou moins continus. Il s'agit alors de tendre vers la représentation des phénomènes multiphysiques [13] (aspects thermomécaniques et physico-chimiques) pouvant être causes de variations de continuité de la matière.

## Remerciements

Les modélisations par éléments discrets ont pu être réalisées par l'intermédiaire de la plateforme logicielle LMGC90 [14].

## Références

[1] M. Godet, *The third-body approach : a mechanical view of wear*, Wear, 100, 437-452, 1984.

[2] S. Descartes, Y. Berthier, *Rheology and flows of solid third bodies : background and application to an MoS<sub>2</sub> coating*, Wear, 252, 546-556, 2002.

[3] J. F. Jerier, J. F. Molinari, *Normal contact between rough surfaces by the discrete element method*, Tribology International, 47, 1-8, 2012.

[4] 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.

[5] Y. Berthier, H.-P. Cao, V.-H. Nhu, M. Renouf, A. Saulot, *Couplage et dialogue entre méthodes discrète et continue pour la tribologie numérique*, Actes du 9<sup>ème</sup> colloque national en calcul des structures, 2009.

[6] F. Perales, *Fissuration des matériaux à gradient de propriétés. Application au zircaloy hydruré*, Thèse de doctorat, Université Montpellier II, 2005.

[7] Q. Zhou, H. Zhang, Y. Zheng, *A homogenization technique for heat transfer in periodic granular materials*, Advanced Powder Technology, 23, 104-114, 2012.

[8] I. Terreros, I. Iordanoff, J. L. Charles, *Simulation of continuum heat conduction using DEM domains*, Computational Materials Science, 69, 46-52, 2013.

[9] W. L. Vargas, J. J. McCarthy, *Heat conduction in granular materials*, AIChE Journal, Wiley Subscription Services, Inc., A Wiley Company, 47, 1052-1059, 2001.

[10] W. L. Vargas, J. J. McCarthy, *Thermal expansion effects and heat conduction in granular materials*, Phys. Rev. E, American Physical Society, 76, 041301-1-041301-8, 2007.

[11] M. M. Yovanovitch, *Thermal contact resistance across elastically deformed spheres*, Journal of Spacecraft and Rockets, 4, 119-122, 1967.

[12] G. K. Batchelor, R. W. O'Brien, *Thermal or electrical conduction through a granular material*, Proceedings of the Royal Society of London, A 355, 313-333, 1977.

[13] M. Renouf, H.-P. Cao, V.-H. Nhu, *Multiphysical modeling of third-body rheology*, Tribology International, 44, 417-425, 2011.

[14] F. Dubois, M. Jean, *LMGC90 : une plateforme de développement dédiée à la modélisation des problèmes d'interaction*, Actes du 6<sup>ème</sup> colloque national en calcul de structures, 1, 111-118, 2003.