# Modélisation du frittage par une approche dynamique des contacts

Sylvain Martin, Mohamed Guessasma, Jérôme Fortin, Khashayar Saleh, Jacques Lechelle, Frédéric Adenot

► To cite this version:

Sylvain Martin, Mohamed Guessasma, Jérôme Fortin, Khashayar Saleh, Jacques Lechelle, et al. Modélisation du frittage par une approche dynamique des contacts. 11e colloque national en calcul des structures, CSMA, May 2013, Giens, France. hal-01717079

HAL Id: hal-01717079

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

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
# Modélisation du frittage par une approche dynamique des contacts

Sylvain MARTIN<fnref n="1" />,<fnref n="3" />*, Mohamed GUESSASMA<fnref n="2" />, Jérôme FORTIN<fnref n="2" />, Khashayar SALEH<fnref n="3" />, Jacques LECHELLE<fnref n="1" />, Frédéric ADENOT<fnref n="1" />

**Résumé** — Le frittage est un procédé clé dans la fabrication des matériaux céramiques. Différents modèles ont été proposés dans la littérature, à l'échelle du grain ainsi que d'un empilement de particules, avec des hypothèses diverses sur l'évolution géométrique et les modes de diffusion prédominants. Depuis quelques années, plusieurs approches utilisant la Méthode des Éléments Discrets ont été proposées. L'avantage est de prendre en compte la nature granulaire des matériaux et de pouvoir accéder facilement aux paramètres caractéristiques de la microstructure. L'objectif de cette étude est d'intégrer le modèle discret de frittage à une méthode implicite de type Dynamique des Contacts, afin de s'affranchir de l'augmentation arbitraire de la masse volumique des particules qui est indispensable lorsque l'on utilise la Dynamique Moléculaire. Les résultats obtenus montrent une bonne cohérence avec les données de la littérature.

**Mots clés** — Frittage, Eléments Discrets, Dynamique des Contacts, Diffusion

## 1. Introduction

Le frittage est un procédé de solidification et de densification des poudres sous l'action de la chaleur qui intervient dans le procédé de fabrication des matériaux céramiques. Le phénomène physique principal est la diffusion à l'état solide. De nombreux modèles ont été développés à différentes échelles. Depuis quelques années, plusieurs simulations du frittage par la Méthode des Éléments Discrets ont été proposées dans la littérature [1-7]. Elles permettent de prendre en compte la nature granulaire du matériau en définissant les interactions entre les grains. Les calculs reposent sur un modèle géométrique simplifié de frittage qui décrit l'évolution sous forme de sphères tronquées [8]. La densification est représentée par l'interpénétration des particules sphériques. Ainsi, il est possible d'évaluer l'impact de paramètres comme la densité initiale ou la distribution granulométrique sur l'évolution de la microstructure. Cependant, la Méthode des Éléments Discrets généralement mise en œuvre dans ces simulations est de type Dynamique Moléculaire. Cette approche a l'avantage d'autoriser l'utilisation de différentes lois de contact. En revanche, le schéma explicite utilisé pour la résolution impose un pas de temps petit. Dans le cas du frittage, les forces mises en jeu sont très élevées devant l'échelle de la particule. Le pas de temps critique est donc très petit, typiquement de l'ordre de 1.10-15s pour des particules métallique d'une dizaine de micromètres [7]. Pour contourner ce problème, la masse volumique des particules est augmentée artificiellement de plusieurs ordres de grandeurs. Afin de s'affranchir de cette contrainte, nous avons proposé une méthode alternative basée sur la Dynamique des Contacts [9,10] qui utilise un schéma de résolution implicite. Il est alors possible de résoudre le modèle discret de frittage avec la véritable masse volumique des particules. Le modèle a été introduit dans le code de calcul MULTICOR [11-13] développé à l'UPJV.

<footnote n="1">1 CEA Cadarache, Sylvain.MARTIN3@cea.fr, jacques.lechelle@cea.fr, frederic.adenot@cea.fr</footnote>

<footnote n="2">2 UPJV, Laboratoire des Technologies Innovantes EA3899, mohamed.guessasma@u-picardie.fr, jerome fortin@u-picardie.fr</footnote>

<footnote n="3">3 UTC, Transformation Intégrée de la Matière Renouvelable EA4297, khashayar.saleh@utc.fr</footnote>

<footnote>* Auteur correspondant</footnote>
## 2. Modélisation du frittage par une approche discrète

### 2.1 Modélisation physique du frittage

Le modèle utilisé dans les modèles discrets de frittage repose sur le modèle géométrique de Coble qui décrit l'évolution des particules sous forme de sphères tronquées [8]. Le modèle physique est défini pour des grains monocristallins. Les diffusions aux joints de grains et en surface sont considérées comme les mécanismes prédominants. Les particules en contact sont soumises à une force, appliquée au point de contact, qui dépend des paramètres caractéristiques du contact et de l'avancement du frittage, Figure 1.

<figure><img src="image_2.png" /><figcaption>Fig. 1. Représentation du contact de deux particules avec formation d'un joint de grain</figcaption></figure>

L'angle dièdre $\Psi$ est supposé constant, ce qui revient à considérer que l'équilibre mécanique, qui dépend du rapport des énergies de surface et du joint de grain, est atteint à chaque instant. Considérons deux particules sphériques en contact, la force de contact normale s'écrit alors [3] :

$$f_n^{fritt} = \frac{\pi r^4}{2\Delta_b \beta} \times V_n - \frac{\pi \gamma_s R \alpha}{\beta}$$

Avec : : (1)

$$\Delta_b = \frac{\Omega}{k_b T} \delta_b D_b$$

R est le rayon des particules, $\gamma_s$ l'énergie de surface et $V_n$ la vitesse de rapprochement du centre des particules. $D_b$ est le coefficient de diffusion aux joints de grains, $\delta_b$ son épaisseur et $\Omega$ le volume atomique. Les coefficients $\alpha$ et $\beta$ dépendent du rapport des coefficients de diffusion aux joints de grains et surfacique [14]. r est le rayon du cou, il évolue, d'après le modèle de Coble [8], tel que:

$$\frac{dr}{dt} = \frac{R}{r} \times V_n \qquad (2)$$

Notons que l'angle dièdre intervenait dans la première expression de la force normale proposée par Parhami [1], ce qui n'est pas le cas dans notre équation (1). En effet, Martin [3], a simplifié l'équation initiale en s'appuyant sur les résultats de Bouvard [14] qui montraient une faible influence de l'angle dièdre sur la cinétique de densification.

Le premier terme de (1) est une résistance visqueuse qui s'oppose à l'interpénétration et représente la cinétique de diffusion au joint de grain. Le second terme représente la force attractive de frittage qui dérive de l'enthalpie libre.

La loi de frottement tangentiel s'écrit :
$$f_t = \eta_c \times \frac{\pi \times r^2 \times R^2}{8D_b} \times \dot{u}_t \quad (3)$$

Avec $\eta_c$ un paramètre inférieur à 0.3 qui dépend des propriétés du joint de grain mais pour lequel il n'existe pas de consensus dans la littérature [1,3,6].

Dans le cas du contact entre deux particules de rayons $R_1$ et $R_2$, il a été montré [15] que R peut simplement être remplacé dans (1,2,3) par :

$$R = 2 \times \frac{R_1 \times R_2}{R_1 + R_2} \quad (4)$$

## 2.2 Utilisation de la Dynamique des Contacts

La méthode des éléments discrets permet de modéliser des systèmes multi-corps qui interagissent entre eux par l'intermédiaire d'une loi de contact. Dans le cas de particules sphériques, chaque particule est repérée par les coordonnées de son centre ainsi que par les angles d'Euler définis dans le repère global de l'empilement. On associe, à chaque contact, un repère local qui permet de définir la vitesse de contact ainsi que les deux composantes normale et tangentielle de la force de contact, figure 2.

<figure><img src="image_3.png" /><figcaption>Fig. 2. Système Multi-contacts</figcaption></figure>

On distingue plusieurs familles d'Éléments Discrets. La méthode Smooth, plus généralement associée à la dynamique moléculaire, regroupe les méthodes explicites. Les méthodes dites Non Smooth comme la Dynamique des Contacts permettent l'utilisation d'une loi de contact non régulière avec un schéma de résolution implicite. Ainsi, le pas de temps utilisé peut être très supérieur à celui de la dynamique moléculaire. Les modèles existants de frittage sont basés sur la dynamique moléculaire. La loi de contact est définie par les forces normales et tangentielles de frittage (1) et (3). Pour pouvoir utiliser la Dynamique des Contacts, nous allons réécrire les forces normale et tangentielle de frittage afin de les intégrer au schéma de résolution implicite.
### 2.2.1 Loi de contact normale

Dans le cas le plus simple de la Dynamique des Contacts, les particules sont rigides et non pénétrables. Ces hypothèses, plus connues sous le nom de conditions de Signorini s'écrivent, dans le repère local du contact entre deux particules :

$$ \dot{u}_n \ge 0 \quad f_n \ge 0 \quad \dot{u}_n \times f_n = 0 \qquad (5) $$

Avec $\dot{u}_n$ la vitesse normale de contact dans le repère local du contact et $f_n$ la force normale unilatérale répulsive.

Dans le cas du frittage, la loi de contact normale s'écrit sous la forme d'une composante attractive et d'une résistance visqueuse à l'interprénétration. Il faut donc réécrire la loi de contact de frittage sous la forme de conditions unilatérales répulsives. En traitant séparément les deux termes de la force de frittage et en isolant la vitesse de rapprochement dans (2), on obtient :

$$ \begin{gathered} (\dot{u}_n + V_n) \ge 0 \quad f_n \ge 0 \\ (\dot{u}_n + V_n) \times f_n = 0 \end{gathered} \qquad \text{Avec} \qquad V_n = \frac{f_n}{\frac{\pi r^4}{2D_b \beta}} \qquad (6) $$

$$ f_{\text{cohe}} = - \frac{\pi \gamma_s R \alpha}{\beta} $$

L'expression de $V_n$ est assez similaire aux méthodes de résolution statiques [1] pour lesquelles l'équation (1) s'écrit sous la forme d'un système matriciel. On peut d'ailleurs retrouver l'équation quasi-statique, en considérant l'accélération nulle : on obtient la force normale répulsive égale et opposée à la force de cohésion, ce qui correspond bien au système utilisé par Parhami [1].

Cependant, ces techniques sont limitées à des empilements de quelques centaines de particules en trois dimensions [6] car les particules flottantes (sans contact) rendent la matrice singulière et donc plus longue à inverser, ce qui n'est pas le cas avec la Dynamique des Contacts.

La loi de contact s'écrit donc comme une translation des conditions usuelles de Signorini. Le graphe de cette relation est donné par la figure 3.

<figure><img src="image_4.png" /><figcaption>Fig. 3. Loi de contact de frittage représentée comme une translation des conditions unilatérales répulsives.</figcaption></figure>

Notons que, sur la figure 3, la translation selon l'axe des ordonnées, qui représente la force attractive, correspond au cas classique des matériaux granulaires cohésifs en Dynamique des Contacts. Cette approche a notamment été utilisée pour représenter les matériaux cimentaires [10,13].
Pour les systèmes multi-corps, la résolution se fait selon un schéma classique avec une itération sur tous les contacts jusqu'à ce que le critère de convergence fixé soit vérifié. La vitesse $V_n$ est calculée avec la force $f_n$ de l'itération précédente Ainsi, à l'itération i, le schéma de résolution local sur chaque contact s'écrit:

$$V_n^i = \frac{f_n^i}{\frac{\pi r^4}{2D_b\beta}} \qquad f_{cohe} = -\frac{\pi\gamma_s R\alpha}{\beta}$$
$$\tau_n^{i+1} = f_n^i + (\dot{u}_n^i + V_n^i) \times \rho \qquad (7)$$
$$Si \quad \tau_n^{i+1} \le 0 \quad \text{alors} \quad f_n^{i+1} = 0$$
$$\text{Sinon} \quad f_n^{i+1} = \tau^{i+1}$$

Le paramètre $\rho$ doit être strictement positif pour assurer la convergence. Dans le cas usuel, nous prendrons $\rho = \frac{m}{\Delta t}$ avec m la masse réduite des particules en contact.

Notons que c'est la somme de la force unilatérale répulsive $f_n^{i+1}$ et de la force attractive $f_{cohe}$ qui est utilisée pour corriger la vitesse des particules dans le repère global.

Ce schéma local de prédiction correction correspond à la formulation obtenue par le bi-potentiel [11,12] pour un coefficient de frottement de Coulomb nul.

### 2.2.2 Loi de contact tangentielle

La loi de contact tangentielle utilisée est généralement le frottement sec de Coulomb. Ici, la force de frottement s'écrit comme une résistance visqueuse au glissement.

L'algorithme de résolution dans le repère local s'écrit sous la forme d'une série $(\dot{u}_t^{(i)}, f_t^{(i)})$, qui représente les valeurs de la vitesse tangentielle et de la force visqueuse du contact à l'itération i. Il faut associer à $(\dot{u}_t^{(i)}, f_t^{(i)})$, les caractéristiques de la résistance visqueuse. Tout d'abord la définition nous permet d'écrire la convergence de la série pour chaque contact :

$$\lim_{i \to \infty} f_t^{(i)} = -\eta' \times \dot{u}_t^{(i)}$$
$$\text{Avec} \qquad \eta' = \eta_c \times \frac{\pi \times r^2 \times R^2}{8D_b} \qquad (8)$$

On peut exprimer la résolution implicite de la force tangentielle sous la forme d'un système de deux équations :

$$f_t^{(i+1)} = f_t^{(i)} + \Delta f_t = -\eta' \times \dot{u}_t^{(i+1)}$$
$$\dot{u}_t^{(i+1)} = \dot{u}_t^{(i)} + \Delta f_t \times \frac{\Delta t}{m} \qquad (9)$$

En combinant ces deux équations, on obtient la forme générale de $\Delta f_t$ :
$$ \Delta f_t = \frac{-\eta'(\dot{u}_t^{(i)} - f_t^{(i)} \times \frac{\Delta t}{m})}{1 + \frac{\eta \Delta t}{m}} - f_t^{(i)} \quad (10) $$

Il est facile de vérifier que les conditions de convergence sont respectées. Si l'on reprend l'expression de $f_t^{(i+1)}$, il est possible d'ajouter un paramètre de relaxation $\omega$ pour faciliter la convergence :

$$ f_t^{(i+1)} = f_t^{(i)} + \omega \times \Delta f_t \quad (11) $$

### 3. Application numérique : frittage de poudres céramiques

Ce modèle a été implémenté dans les versions 2D [11] et 3D [12] du logiciel MULTICOR. Les empilements aléatoires en 3D ont été générés à avec le logiciel open source LIGGGHTS [16]. Afin de comparer l'évolution du frittage pour différentes poudres, nous avons utilisé des variables normalisées [14]. Dans un premier temps, le frottement visqueux n'a pas été intégré aux simulations.

Nous avons simulé le frittage d'un empilement bimodal de 1500 particules en deux dimensions (Figure 4). L'évolution moyenne du rayon du cou normalisé a été comparée avec le modèle de Coblentz [17] qui prédit l'évolution pour le frittage de deux particules isolées:

$$ \frac{r}{R} = \left( \frac{54t}{\tau_g} \right)^{\frac{1}{6}} \quad \text{avec } \tau_g = \frac{R^4}{\Delta_b} \quad (12) $$

<figure><img src="image_5.png" /><figcaption>Fig. 4. Frittage d'un empilement bimodal de 1500 particules en 2D</figcaption></figure>

Le visu de l'empilement représente la vitesse associée à chaque particule dans le repère global de l'empilement. Les particules se déplacent en direction du centre de l'empilement conformément à la logique de densification. Le graphique représente la moyenne ainsi que les extremums de rayons du cou. Il apparaît que les contacts de l'empilement se comportent en moyenne comme un contact isolé. La courbe minimum stagne à une valeur très faible en raison de l'apparition régulière de nouveaux contacts. La différence entre les extrémums témoigne de l'hétérogénéité des contacts sur l'empilement. La résistance visqueuse à l'interpénétration, qui correspond au premier terme de l'équation (1), varie en fonction du rayon du cou à la puissance quatre. La coexistence de contacts
avec des degrés d'avancement très différents amène un écart de plusieurs ordres de grandeurs pour la résistance visqueuse à l'interprénétration, ce qui est à l'origine de difficultés numériques de résolution.

L'introduction du frottement visqueux ralentit la convergence et multiplie par deux environ le nombre d'itérations nécessaires pour obtenir une précision identique. Dans le cas d'un empilement poly-disperse en trois dimensions de 17000 particules (Figure 5), il faut environ 5000 itérations pour obtenir une précision de l'ordre de 1/10000 sur le critère de convergence. D'après les premiers résultats, le temps de calcul demeure raisonnable pour des empilements inférieurs à 20000 particules.

<figure><img src="image_6.png" /><figcaption>Fig. 5. Frittage d'un empilement polydisperse de 17000 particules en 3D</figcaption></figure>

La vitesse de frittage diminue avec la taille des particules. Le pas de temps devant être adapté à la vitesse de frittage des plus petites particules, il décroît avec l'amplitude de la distribution granulométrique. La modélisation d'empilement avec une distribution granulométrique étendue est, par conséquent, particulièrement consommatrice de temps de calcul.

## 4. Conclusion

L'algorithme proposé permet de modéliser le frittage des poudres à l'état solide, selon le modèle classique des sphères tronquées, en s'affranchissant de l'augmentation de la masse inhérente à l'utilisation des méthodes explicites. Cela devrait permettre d'étudier de façon plus précise l'évolution de la microstructure et notamment le réarrangement qui est très sensible à la densité des particules. L'intérêt est d'autant plus grand pour les particules submicroniques pour lesquelles l'échelle des forces de contact est très grande devant la taille des particules. D'autres études devront être menées pour mettre en évidence cet aspect.

La principale limite est l'augmentation du temps de calcul. En effet, contrairement aux méthodes explicites qui ne nécessitent pas d'itération pour calculer les forces de contact, dans le cas présent, il faut quelques milliers d'itérations pour obtenir une précision raisonnable, ce qui augmente d'autant la durée du calcul. Des empilements poly-disperses de 20000 particules maximum ont été utilisés avec la version 3D du code. Plusieurs pistes sont étudiées afin de repousser cette limite, comme l'optimisation de l'initialisation des forces de contact ou la parallélisation du code.

## Références

[1] F. Parhami, R. M. McMeeking, A network model for initial stage sintering, Mechanics of Materials, Elsevier, 1998

[2] C. L. Martin, L. C. R. Schneider, L. Olmos, D. Bouvard, Discrete element modelling of metallic powder sintering, Scripta Materialia, Elsevier, 2006
[3] C. L. Martin, R. K. Bordia, The effect of a substrate on the sintering of constrained films, Acta Materialia, Elsevier, 2009

[4] C. L. Martin, H. Camacho-Montes, L. Olmos, D. Bouvard, R. K. Bordia, Evolution of defects during sintering: Discrete Element Simulations, Journal of the American Ceramic Society, Wiley & Sons, 2009

[5] L. R. Olmos, C. L. Martin, F. Bouvard, D. Bellet, M. Di Michiel, Investigation of the sintering of heterogeneous Powders Systems by Synchrotron Microtomography and Discrete Element Simulation, Journal of the American Ceramic Society, 2009

[6] B. Henrich, A. Wonish, T. Kraft, M. Moseler, H. Riedel, Simulation of the influence of rearrangement during sintering, Acta Materialia, Elsevier, 2009

[7] L. R. Olmos Navarrete, Etude du frittage de poudres par microtomographie in situ et modélisation discrète, thèse de l'Institut Polytechnique de Grenoble, 2010

[8] R. L. Coble, Initial Sintering of Alumina and Hematite, Journal of the American Ceramic Society, John Wiley & Sons, 1958

[9] M. Jean, The non Smooth Contact Dynamics Methods, Computer methods in applied mechanics and engineering, Elsevier, 1999

[10] M. Jean, V. Acary, Y. Monerie, Non-smooth contact dynamics approach of cohesive materials, Philosophical transactions of the Royal Society A, Royal Society Publishing, 2001

[11] J. Fortin, O. Millet, G. de Sacxcé, Numerical simulation of granular materials by an improved discrete element method, International Journal of Numerical Methods in Engineering, John Wiley & Sons, 2005

[12] I. Sanni, E. Bellenger, J. Fortin, P. Corevits, A reliable algorithm to solve 3D frictional multi-contact problems: Application to granular media, Journal of Computational and Applied Mathematics, Elsevier, 2010

[13] M. Guessasma, J. Fortin, Modélisation d'un milieu discret cohésif par l'approche dynamique des contacts : Application à la rupture d'un matériau cimentaire, CSMA, 2011

[14] D. Bouvard, R. M. McMeeking, Deformation of Interparticle Necks by Diffusion-Controlled Creep, Journal of the American Ceramic Society, John Wiley & Sons, 1996

[15] F. Parhami, R. M. McMeeking, A. C. F. Cocks, Z. Suo, A model for the sintering and coarsening of rows of spherical particles, Mechanics of Materials, Elsevier, 1998

[16] C. Kloss, C. Goniva, A. Hager, S. Amberger, S. Pirker, Models, algorithms and validation for opensource DEM and CFD-DEM, Progress in Computational Fluid Dynamics, 2012

[17] W. S. Coblenz, J. M. Dynys, R. M. Cannon, R. L. Coble, Initial Stage Sintering Models: A Critical Analysis and Assessment, Proceedings of the Fifth International Conference on Sintering and Related Phenomena, G. C. Kuzynski Plenum Press, 1980