# Modélisation et simulation du transport de matière par diffusion surfacique à l'aide d'une approche Level-Set

Julien Bruchon, Guillaume Pacquaut, Sylvain Drapier, François Valdivieso

► To cite this version:

Julien Bruchon, Guillaume Pacquaut, Sylvain Drapier, François Valdivieso. Modélisation et simulation du transport de matière par diffusion surfacique à l'aide d'une approche Level-Set. 9e Colloque national en calcul des structures, CSMA, May 2009, Giens, France. hal-01413789

HAL Id: hal-01413789

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

Submitted on 11 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_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Modélisation et simulation du transport de matière par diffusion surfacique à l'aide d'une approche Level-Set

J. Bruchon¹, G. Pacquaut¹, S. Drapier¹, F. Valdivieso²

Ecole Nationale Supérieure des Mines de Saint-Etienne
Centre Sciences des Matériaux et des Structures

¹ LTDS CNRS UMR 5513

² PECM CNRS UMR 5146

158 cours Fauriel

42024 Saint-Etienne - Cedex 2

bruchon@emse.fr

**Résumé** — Dans le contexte général du frittage entre grains de céramique, cet article présente une stratégie numérique pour la simulation directe du transport de matière par diffusion surfacique. Une méthode Level-Set est utilisée pour décrire l'évolution de la surface libre des grains. La vitesse d'interface est alors fonction du Laplacien surfacique de la courbure qu'il faut donc calculer rigoureusement. Plusieurs exemples de simulations sont présentés.

**Mots clés** — diffusion surfacique, méthode Level-Set, frittage.

## 1 Introduction

Le cadre général de ce travail est la simulation directe de la formation et de l'évolution de la microstructure d'une céramique lors de son frittage. Ce procédé d'élaboration permet aux grains d'une poudre compactée de densifier, à une température bien inférieure à celle de fusion, éliminant ainsi la porosité. Classiquement [1], le transport de matière durant le frittage s'effectue de la surface du grain vers le cou (Figure 1), et du centre du joint de grain vers le cou, via des chemins de diffusion surfaciques, volumiques, de joint de grain, ou par le milieu environnant (gaz).

Dans ce premier travail, nous traitons le cas du transport de matière de la surface vers le cou par diffusion surfacique. Ceci nous permet de nous affranchir dans un premier temps de la réponse mécanique du matériau, puisque la diffusion est uniquement guidée par un critère géométrique : le Laplacien de la courbure. Dans la suite, nous présentons la modélisation physique du phénomène de diffusion surfacique, puis la stratégie numérique mise en oeuvre. Nous conclurons par deux applications d'évolutions d'interfaces par diffusion surfacique de la matière.

<figure><img src="image_2.png" /><figcaption>Figure 1 – Transfert de matière par diffusion surfacique entre deux grains identiques.</figcaption></figure>
## 2 Modèle physique

Considérons un ensemble de grains $\Omega_g$, et notons par $S_g$ la surface libre de cet ensemble (voir Figure 1). Nous supposons, comme dans [1] et [2], qu'il existe le long de $S_g$ un flux de matière $\mathbf{j}_s$ induit par des gradients de potentiel chimique, lesquels dépendent à leur tour du gradient de la courbure moyenne $K$ de la surface :

$$ \mathbf{j}_s = \frac{\delta_s D_s \Omega \gamma_s}{kT} \nabla_s K \quad (1) $$

où $T$ est la température absolue, $k$ la constante de Boltzmann, $D_s$ le coefficient de diffusion surfacique, $\delta_s$ l'épaisseur dans laquelle a lieu la diffusion, $\gamma_s$ l'énergie de surface libre, et $\Omega$ le volume molaire du matériau. Enfin, $\nabla_s$ désigne l'opérateur de gradient surfacique.

Pour déterminer la vitesse d'interface engendrée par ce flux, prenons une surface élémentaire $S \subset S_g$. Nous faisons alors l'hypothèse que toute matière accumulée sur cette surface par unité de temps entraîne un déplacement normal de la surface. Ceci se traduit par une vitesse de surface $\mathbf{v}_s$ normale : $\mathbf{v}_s = v_n \mathbf{n}$, avec $\mathbf{n}$ la normale unitaire extérieure à $S$. L'égalité entre le volume de matière traversant la frontière $\partial S$ et le déplacement induit s'écrit :

$$ \int_{\partial S} \mathbf{j}_s \cdot \mathbf{n}_l dl = \int_S \mathbf{v}_s \cdot \mathbf{n} dS \quad (2) $$

où $\mathbf{n}_l$ est le vecteur normal à $\partial S$ (situé dans le plan tangent de $S$). D'où, en appliquant le théorème de la divergence :

$$ \int_S [\nabla_s \cdot \mathbf{j}_s - (v_n \mathbf{n}) \cdot \mathbf{n}] dS = 0 \quad (3) $$

où $\nabla_s \cdot$ est la divergence surfacique. Ceci étant valide pour tout $S \subset S_g$, la vitesse de l'interface induite par le flux de matière $\mathbf{j}_s$ s'écrit :

$$ \mathbf{v}_s = (\nabla_s \cdot \mathbf{j}_s)\mathbf{n} = C_0(\Delta_s K)\mathbf{n} \quad (4) $$

avec $C_0 = \frac{\delta_s D_s \Omega \gamma_s}{kT}$ que l'on suppose ici constant.

Remarquons pour conclure que le volume des grains reste inchangé par le transport par la vitesse $\mathbf{v}_s$. En effet, la variation de volume des grains vaut : $d|\Omega_g|/dt = \int_{S_g} \mathbf{v}_s \cdot \mathbf{n} dS = \int_{S_g} \nabla_s \cdot \mathbf{j}_s dS$. Donc, comme $S_g$ est une surface fermée, que $\mathbf{j}_s \cdot \mathbf{n} = 0$ et que $\nabla \mathbf{n} \cdot \mathbf{n} = 0$ :

$$ \frac{d|\Omega_g|}{dt} = - \int_{S_g} \mathbf{j}_s \cdot \nabla 1 dS + \int_{S_g} (\mathbf{j}_s \otimes \mathbf{n}) : \nabla \mathbf{n} = 0 $$

## 3 Modèle numérique

Nous discrétisons l'Equation (4) à l'aide d'une méthode éléments finis. Appelons $\Omega$ le domaine de calcul, le carré unité en 2D ou le cube unité en 3D. $\Omega$ est composé d'une partie occupée par les grains, $\Omega_g$, et de la partie complémentaire, occupée par l'air, $\Omega_a : \Omega = \Omega_g \cup \Omega_a$. Le domaine $\Omega$ est discrétisé par un maillage $\mathcal{T}_h(\Omega)$ de simplexes (des triangles en 2D, des tétraèdres en 3D). Le choix de discrétiser l'ensemble du domaine de calcul, et non uniquement la surface des grains, vient de ce que notre objectif à terme est de simuler le procédé de frittage à l'échelle des grains, incluant une diffusion surfacique de la matière, mais aussi une diffusion en volume induite par le gradient des contraintes dans le matériau. D'où, la nécessité de représenter les grains globalement et non uniquement d'un point de vue surfacique.

De plus, nous choisissons une approche Eulérienne de la description des grains, à savoir que ceux-ci sont décrits par une fonction $\phi : \Omega \to \mathbb{R}$ qui évolue au travers du maillage fixe $\mathcal{T}_h(\Omega)$.
Ainsi, les changements topologiques pouvant intervenir dans la surface des grains ("fusion" de deux surfaces) seront pris en compte implicitement.

## 3.1 Description des grains

Nous utilisons la méthode Level-Set développée dans [3] pour décrire et convecter l'ensemble des grains. Nous cherchons donc, pour tout temps $t$ et tout $x \in \Omega$, la valeur de $\phi(x,t)$ telle que $\phi(x,t) > 0$ si $x \in \Omega_a$ et $\phi(x,t) < 0$ si $x \in \Omega_g$. La surface libre des grains est alors définie par

$$S_g(t) = \{x \in \Omega; \phi(x,t) = 0\}$$

A l'instant initial $t=0$, nous choisissons $\phi(x,0) = \phi_0(x)$ avec :

$$\phi_0(x) = \begin{cases} \frac{2E}{\pi} \sin\left(\frac{\pi d(x, S_g(0))}{2E}\right) & \text{si } |d(x, S_g(0))| \le E \\ +\frac{2E}{\pi} & \text{si } d(x, S_g(0)) \ge E \\ -\frac{2E}{\pi} & \text{si } d(x, S_g(0)) \le -E \end{cases} \quad (5)$$

où $d(x, S_g(0))$ est la distance signée du point $x$ à la surface libre initiale $S_g(0)$, positive si $x \in \Omega_a$, négative sinon. Le paramètre numérique $E$ est une épaisseur, choisie dans la pratique égale à $20 \times h$, $h$ étant la taille de maille. Remarquons que dans le voisinage de l'interface $S_g(0)$, $\phi_0(x) \approx d(x, S_g(0))$ car $\sin b \approx b$ lorsque $|b|$ est petit. De plus, $\phi_0$ est constante au-delà de la distance $E$ de l'interface, tout en étant dérivable sur tout $\Omega$. Pour tout $t > 0$, nous construisons la fonction $\phi(.,t)$ qui hérite de ces trois propriétés, ce qui nous permettra d'une part d'obtenir la courbure correcte dans un voisinage de l'interface, et d'autre part de ne pas trop nous soucier de ce qui se passe au-delà de la distance $E$ de l'interface, puisque $\phi$ sera alors constante.

Pour ce faire, $\phi$ est la solution de l'équation de transport "classique" dans un voisinage de l'interface

$$\frac{\partial \phi}{\partial t} + \mathbf{v}_s \cdot \nabla \phi = 0 \quad (6)$$

tandis qu'en dehors de ce voisinage, $\phi$ est cherchée telle que

$$|\nabla \phi| = \sqrt{1 - \left(\frac{\pi \phi}{2E}\right)^2}$$

Ceci est l'analogue de la réinitiation classique $|\nabla \phi| = 1$, mais pour une fonction $\phi$ du type (5) (en effet, $\nabla \sin(x) = \cos(x) = \sqrt{1 - \sin^2(x)}$). Ainsi, les étapes de convection et de réinitialisation sont réunies en une seule équation, comme décrit dans [3].

Nos simulations sont effectuées en prenant $\phi_h$, fonction approchant $\phi$, comme étant continue sur tout $\Omega$ et linéaire par élément.

## 3.2 Calcul de la courbure

La courbure moyenne est définie par (voir [4]):

$$K = \nabla \cdot \frac{\nabla \phi}{\| \nabla \phi \|}$$

Cependant, la fonction approchée $\phi_h$ étant linéaire par élément, son gradient est constant par élément, et son Laplacien est donc nul. Nous considérons donc la forme faible de l'équation suivante :
$$ K + \varepsilon_K \Delta K = \nabla \cdot \frac{\nabla \phi}{\| \nabla \phi \|} \qquad (7) $$

où le terme de diffusion $\varepsilon_K \Delta K$ a été ajouté afin de régulariser la courbure, avec $\varepsilon_K > 0$, $\varepsilon_K \ll 1$. La formulation faible discrète que nous résolvons s'écrit alors :

$$ \int_{\Omega_h} K_h q_h d\Omega - \varepsilon_K \int_{\Omega_h} \nabla K_h \cdot \nabla q_h d\Omega = - \int_{\Omega_h} \frac{\nabla \phi_h}{\| \nabla \phi_h \|} \cdot \nabla q_h d\Omega \quad (8) $$

pour toute fonction test $q_h \in H^1(\Omega)$. L'équation (8) a été obtenue en supposant que $S_g$ est à une distance supérieure à $E$ de $\partial\Omega_h$. En effet, nous avons alors $\phi_h|_{\partial\Omega_h} = 2E/\pi$, et donc $\nabla\phi_h|_{\partial\Omega_h} = 0$. Les intégrales surfaciques s'annulent alors. De plus, nous devons souligner que $K_h$, solution de (8), est définie sur tout le maillage, et correspond dans un voisinage de l'interface à la courbure réelle de l'interface.

<figure><img src="image_3.png" /><figcaption>Figure 2 – Isovaleurs {$K_h = 5$} pour $\varepsilon_K = 0$ et $\varepsilon_K = 10^{-4}$, superposées à l'interface $S_g$, ici le cercle de rayon 0,2.</figcaption></figure>

Dans nos simulations, $K_h$ est choisie continue et linéaire par élément. La figure 2 montre l'isovaleur {$\phi_h = 0$} correspondant à un cercle de rayon 0,2, ainsi que les isvaleurs {$K_h = 5$} de la courbure (la courbure "théorique" du cercle), obtenues respectivement avec $\varepsilon_K = 0$ et $\varepsilon_K = 10^{-4}$. Nous voyons que dans le premier cas l'isovaleur 5 de la courbure est très irrégulière, tandis que dans le second cas elle est bien localisée à l'interface $S_g$. D'autre part, notons la présence d'une seconde isovaleur 5 de la courbure, située près du centre du grain, liée au fait que la fonction level-set $\phi$ définie plus haut n'est pas deux fois dérivable à la distance $E$ de l'interface, engendrant ainsi un saut dans la courbure. Cependant ceci se passe suffisamment loin de l'interface pour ne pas en perturber l'évolution. Ainsi, cet exemple simple montre que la courbure de l'interface est correctement calculée lorsque l'on résout (8).

### 3.3 Calcul de la vitesse de diffusion surfacique

La vitesse de diffusion surfacique $v_s$ est donnée par l'équation (4). Cependant, puisque nous considérons une description Eulérienne de $S_g$, (4) ne peut pas être directement appliquée. Nous devons prolonger la définition de $v_s$ dans tout le domaine de calcul $\Omega$. Nous avons déjà vu que la définition de la courbure $K_h$ donnée par (8), est valide sur tout $\Omega$. Sachant que l'interface est définie par $S_g = \{\phi = 0\}$, nous exprimons le Laplacien surfacique $\Delta_s$ en fonction des opérateurs classiques $\nabla \cdot$ et $\nabla$ (voir [5]) :

$$ \Delta_s K = \frac{1}{\| \nabla \phi \|} \nabla \cdot (\| \nabla \phi \| P \nabla K) $$
avec P la matrice de projection sur le plan tangent :

$$P = I - \frac{\nabla \phi}{\|\nabla \phi\|} \otimes \frac{\nabla \phi}{\|\nabla \phi\|}$$

où I est le tenseur identité. Au final, en effectuant la même régularisation que pour la courbure, nous obtenons, à partir de (4), l'équation scalaire suivante sur $v_n$ ($v_s = v_n \mathbf{n}$) :

$$v_n \|\nabla \phi\| + \varepsilon_v \Delta v_n = C_0 \nabla \cdot (\|\nabla \phi\| P \nabla K) \quad (9)$$

Nous résolvons alors la formulation faible discrète suivante :

$$\int_{\Omega_h} \|\nabla \phi_h\| v_{nh} q_h d\Omega - \varepsilon_v \int_{\Omega_h} \nabla v_{nh} \cdot \nabla q_h d\Omega = - \int_{\Omega_h} C_0 \|\nabla \phi\| (P_h \nabla K_h) \cdot \nabla q_h d\Omega \quad (10)$$

pour toute fonction test $q_h \in H^1(\Omega)$. Dans nos simulation, $v_{nh}$ est prise continue sur $\Omega_h$, et linéaire par élément. Nous avons alors :

$$\mathbf{v}_{sh} = v_{nh} \mathbf{n}_h \quad (11)$$

où $\mathbf{n}_h$ vaut, au nœud $i$, la moyenne des $\nabla \phi_h / \|\nabla \phi_h\|$ sur les éléments partageant $i$ comme sommet commun.

## 4 Simulation numérique de la diffusion surfacique

Dans les simulations présentées ici, les calculs de $\phi_h$, $K_h$ et $v_{nh}$ sont découplés, à savoir que, connaissant $\phi_h$ au temps $t$, nous en déduisons $K_h$ puis $v_{nh}$ par (8) et (10), avant de transporter $\phi_h$ par la vitesse $\mathbf{v}_{sh}$ obtenue et de passer au temps suivant $t + \Delta t$. Nous avons de plus pris $C_0 = 10^{-4}$, $\Delta t = 10^{-2}$, $\varepsilon_k = 10^{-4}$, $\varepsilon_v = 5 \cdot 10^{-3}$. Le maillage est isotrope et adapté près de l'interface $S_g$ avec une taille de maille minimale variant selon les cas, et une taille maximale de 0,03 (voir [6]).

Tous nos développements ont été effectués dans la librairie éléments finis C++ CIMLIB développée au CEMEF (Centre de Mise en Forme des Matériaux - Ecole des Mines de Paris).

<figure><img src="image_4.png" /><figcaption>Figure 3 – Diffusion surfacique entre deux grains de mêmes tailles : isovaleur $\{\phi_h = 0\}$.</figcaption></figure>

La première simulation met en jeu deux grains initialement de même rayon $r = 0,2$, comme représenté sur la Figure 3(a). Le centre du premier grain est en (0,301;0,5), celui du second grain en (0,699;0,5), si bien que les deux grains "s'interpénetrent" légèrement. La taille de maille minimale est de 0,008. Tout d'abord, remarquons que la diffusion surfacique intervient bien d'une zone de $S_g$ à faible courbure, vers la zone à plus forte courbure, i.e., vers le cou. Ensuite, les zones à courbure constante (les deux extrémités extérieures des grains) ne bougent pas durant
la simulation, comme attendu. Enfin, le volume des grains est conservé, et leurs centres ne se rapprochent pas, comme souhaité (les figures sont à la même échelle).

La deuxième simulation présentée ici, est réalisée avec cinq grains de tailles différentes : quatre de rayon 0,18, et un de rayon 0,1, positionnés comme sur la Figure 4. Les remarques faites précédemment restent valables. De plus, la symétrie du système est préservée. Enfin, le maillage adapté est présenté sur la Figure 4(c) : la taille de maille minimale est ici de 0,003 pour un total de 20 000 noeuds. Le temps $t = 0,05$ a été atteint en 50 incréments, pour un temps de calcul de 5 minutes.

<figure><img src="image_5.png" /><figcaption>Figure 4 – Diffusion surfacique entre 4 grains de rayons 0,18 et 1 grain de rayon 0,1 : isovaleur $\{\phi_h = 0\}$.</figcaption></figure>

# 5 Conclusion et perspectives

Nous avons présenté une stratégie numérique basée sur une approche de type Level-Set permettant de simuler le transport de matière par diffusion surfacique, en combinant adaptation de maillage et régularisation de la courbure et de la vitesse de diffusion. Les premières simulations, effectuées en 2 dimensions, sont encourageantes et montrent un comportement cohérent de l'interface. Les prochaines étapes consistent à prendre en compte la diffusion aux joints de grains, dans le cadre de la simulation directe du frittage, puis, naturellement, la diffusion en volume.

# Références

[1] M.F. Ashby. A first report on sintering diagrams. *Acta Metallurgica*, Elsevier, 275-289, 1974.

[2] D. Bouvard, R.M. McMeeking. Deformation of interparticle necks by Diffusion-Controlled Creep. *Journal of the American Ceramic Society*, Wiley-BlackWell, 665-672, 1996.

[3] T. Coupez. Réinitialisation convective et locale des fonctions Level Set pour le mouvement de surfaces et d'interfaces. *Journées Activités Universitaires de Mécanique, La Rochelle*, ISBN 2-9526-8123-8, 33-41, 2006.

[4] S. Osher, F. Fedkiw. Level Set Methods : An Overview and Some Recent Results. *Journal of Computational Physics*, Elsevier, 463-502, 2001.

[5] M. Burger, F. Hausser, C. Stöcker, A. Voigt. A level-set approach to anisotropic flows with curvature regularization. *Journal of computational physics*, Elsevier, 183-205, 2007.

[6] J. Bruchon, H. Digonnet, T. Coupez. Using a signed distance function for the simulation of metal forming processes : Formulation of the contact condition and mesh adaptation. From a Lagrangian approach to an Eulerian approach. *International Journal for Numerical Methods in Engineering*, Wiley InterScience, pré-publié en ligne, 2009.